From e11bbe80e8d67ac48ccc768b5f814b394880583c Mon Sep 17 00:00:00 2001 From: Ariane Date: Fri, 22 May 2020 20:35:07 +0200 Subject: [PATCH] =?UTF-8?q?Dateien=20hochladen=20nach=20=E2=80=9E=E2=80=9C?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- hausaufgabe2.c | 247 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 247 insertions(+) create mode 100644 hausaufgabe2.c diff --git a/hausaufgabe2.c b/hausaufgabe2.c new file mode 100644 index 0000000..e3d61d3 --- /dev/null +++ b/hausaufgabe2.c @@ -0,0 +1,247 @@ +//gcc hausaufgabe2.c -o hausaufgabe2 -lm -lh2_cip -L. +//./hausaufgabe2 +//Pauline Maas, Ariane Ufer + +#include +#include +#include +#include "h2.h" + +//Globale Variablen definieren, damit nicht so viele Werte übergeben werden müssen +//Hier: Variablen haben Werte aus Aufgabenstellung + +double g = 0.1; //gamma +double a = 0.5; //alpha +double w = 1.5; //omega +double x0 = 1.5; //Startwerte +double v0 = 2.75; +int neq = 2; //Anzahl Gleichungen + +//Definiere die externe Kraft für H5 +double ExterneKraft(double t) { + + return a * cos(w * t); //Rückgabe der Kraft +} + +//Damit man eine RK4 Funktion für alle Aufgaben verwenden kann hier für H4 eine Funktion die als Kraft Null zurückgibt +double KeineKraft(double t) { + return 0; +} + +//rechte Seite der DGL +void dgl(double w0, double *y, double *f, double t, double (*fext)(double)) +{ + f[0] = y[1]; //Definition der Hilfsfunktionen + f[1] = -w0 * w0 * y[0] - 2 * g * y[1] - (fext(t)); +} + +void euler( double h, double w0, double *y, double *f, + void (*dgl)(double, double*, double*, double, double(double)), double t, + double (*fext)(double)) { + int i; + + dgl(w0, y, f, t, fext); // Informationen über rechte Seite holen + + //Euler-Verfahren anwenden und x_n+1 berechnen + for (i = 0; i < neq; i++) { + y[i] += h * f[i]; + } +} + +void RK4(double h, double w0, double *y, double *f, + void (*dgl)(double, double*, double*, double, double(double)), + double t, double (*fext)(double)) { + int i; + double yTemp[neq]; //Temporäre Variable für y und t, die beim RK4 zur Berechnung von k verändert werden + double tTemp; + double *k[4]; //Pointer Array um die 4 zweidimensionalen k's abzuspeichern + for(int i=0;i<4;i++){ //Speicherplatz für die k-Werte zur Verfügung stellen (2 dimensionaler Array) + k[i]=(double*)malloc(sizeof(double)*neq); + } + //k's berechnen: Hier k[0]=k1, k[1]=k2 usw + dgl(w0, y, k[0], t, fext); // dgl berechnet k1 + //for Schleife um mit den k1 die nach RK4 veränderten y und t zu berechnen + for (i = 0; i < neq; i++) { + yTemp[i] = y[i] + h / 2 * k[0][i]; + tTemp = t + h / 2; + } + dgl(w0, yTemp, k[1], tTemp, fext); //dgl berechnet mit verändertem y und t k2 + //for Schleife berechnet neue y und t für k3 nach RK4 + for (i = 0; i < neq; i++) { + yTemp[i] = y[i] + h / 2 * k[1][i]; + tTemp = t + h / 2; + } + dgl(w0, yTemp, k[2], tTemp, fext); //k3 mit Hilfe von Dgl und neuen y und t berechnen + //y und t für k4 berechnen + for (i = 0; i < neq; i++) { + yTemp[i] = y[i] + h * k[2][i]; + tTemp = t + h; + } + dgl(w0, yTemp, k[3], tTemp, fext); //k4 berechnen + //mit allen k's y nach der RK4 vorschrift berechnen + for (i = 0; i < neq; i++) { + y[i] += (h / 2 * (k[0][i] + k[3][i]) + h * (k[1][i] + k[2][i])) / 3; + } + for(int i=0;i<4;i++){ // Speicherplatz wieder freigeben + free(k[i]); + } +} + +//berechnet analytische lösung der Differentialgleichung für H4 +double exactH4(double w0, double t) { + return (exp(-g * t) + * (x0 * cos(sqrt(w0 * w0 - g * g) * t) + + v0 * sin(sqrt(w0 * w0 - g * g) * t))); +} + +//berechnet analytische Lösung der DGL für H5 +double exactH5(double w0, double t) { + //berechnen der Koeffizienten siehe Pdf + double A = a / (sqrt(pow((w * w - w0 * w0), 2) + 4 * g * g * w * w)); + double D = atan(2 * g * w / (w * w - w0 * w0)); //D ist äquivalent zu delta + double C1 = x0 - A * cos(D); + double C2 = (v0 - w * A * cos(D) + g * C1) / sqrt(w0 * w0 - g * g); + + return (exp(-g * t) + * (C1 * cos(sqrt(w0 * w0 - g * g) * t) + + C2 * sin(sqrt(w0 * w0 - g * g) * t)) + A * cos(w * t + D)); + +} + +//berechnet Leistung P(t) für H.6 +double Leistung(double *y, double t,double (*fext)(double)) +{ + double P = fext(t) * y[1]; + return P; +} + +//berechnet den Leistungsmittelwert über eine Anzahl von Perioden für H.6 +int Leistung_Mittelwert(int n, double h, double w0, double *y, double *f, void (*dgl) (double,double*, + double*,double, double(double)), double t, double (*fext)(double),double (*P)(double*, double, double(double))) +{ + double x1,x2; + int j; //Definition Schrittezähler innerhalb einer Periode + double Pges; // Variable zum speichern der addierten Leistungswerte + double P_Mittelwert[n]; //Array zum speichern des Ergebnisses + double P_M_add = 0; //Hilfsvariablen zur Berechnung des Mittelwerts der einzelnen Mittelwerte + double Ges_Mittelwert; + double s_i; // Hilfsvariablen zur Berechnung der Standardabweichung + double s_add = 0; + double Ges_s; + + x1 = y[0]; //entspricht am Anfang x0 + RK4(h,w0,y,f,dgl,t,fext); // mit t = t0 berechnet nächsten x-Wert + x2 = y[0]; + printf("Omega_0 = %20.5le\n",w0); + while (x1 * x2 > 0) //solange kein Nustellenübergang (bzw.Vorzeichenwechsel) zwischen x1 und x2 vorliegt + { + RK4(h,w0,y,f,dgl,t,fext); + x1 = x2; //Wert übergeben + x2 = y[0]; // neuen Wert übernehmen + t += h; + } + RK4(h,w0,y,f,dgl,t,fext); // einen weiteren Schritt berechnen, so dass x1 den ersten Wert nach der NST annimmt (wenn nicht sogar x1 = 0 wird) + x1 = x2; + x2 = y[1]; + t += h; + for (int l = 0; l < n; l++) //wird so oft durchlaufen wie die gewünschte Periodenanzahl n + { + Pges = 0; + j = 0; + for (int i = 0; i < 2; i++)// Zwei Durchläufe, da eine Periode durch jede zweite NST begrenzt wird + { + while (x1 * x2 > 0) // überprüft auf Vorzeichenwechsel (NST) + { + Pges += P(y,t,fext); //Leistung am Punkt x1 wird berechnet und zu Pges hinzugefügt + t += h; + j ++; // zählt die Anzahl der berechneten Leistungswerte + x1 = x2; //Bestimmung der neuen x-Werte + RK4(h,w0,y,f,dgl,t,fext); + x2 = y[0]; + } + Pges += P(y,t,fext); //den letzten Wert vor der NST berechnen + j ++; + RK4(h,w0,y,f,dgl,t,fext); //und einen Schritt weiter gehen, um aus dem Nullstellenübergang herauszukommen + x1 = x2; + x2 = y[1]; + t += h; + } + P_Mittelwert[l] = Pges/ j; + printf("Mittelwert der Leistung (Periode %d): %20.5le\n", l+1 , P_Mittelwert[l]); + P_M_add += P_Mittelwert[l]; + } + Ges_Mittelwert = P_M_add/n; + for (int i = 0; i < n; i++) // Berechnung Standardfehler + { + s_i = pow (P_Mittelwert[i]-Ges_Mittelwert,2); + s_add += s_i; + } + Ges_s = sqrt(s_add/n*(n-1)); + printf("Mittelwert über alle %d Perioden:%20.5le\n",n,Ges_Mittelwert); + printf("Standardfehler:%20.5le\n",Ges_s); + return 0; +} + + +int main() { + double h = 0.001; //Schrittweite h + double t0 = 0; //Startzeit + + double tend = 1; //Endzeit einstellbar + + double *y, *f; //Pointer auf denen die Verktoren y und f gespeichert werden + double t; //Variable ,die momentane Zeit beschreibt + double w0 = 1; //omega null + int n; // Dafiniere Variable für die Periodenanzahl + + /*printf("Bitte geben Sie h, x0, v0 und tend, sowie omega0 und gamma ein:\n"); + scanf("%le %le %le %le %le %le", &h, &t0, &x0, &v0, &tend, &w0, &g);*/ + + y = (double*) malloc(sizeof(double) * neq); //Speicherplatz für y allocieren(y hat dim neq (hier 2)) + f = (double*) malloc(sizeof(double) * neq); //Speicherplatz für f allocieren (dim 2 siehe y) + + //startwerte zuweisen + y[0] = x0; + y[1] = v0; + + //H4: + + printf("t\tLösungen von H.4: x(t)\tdx(t)/dt\n"); + //Für jeden Zeitpunkt RK4 ausführen + for (t = t0; t <= tend; t += h) { + printf(" %20.5le\t%20.5le\t%20.5le\n", t, y[0], y[1]);//Ergebnisse als Tabelle ausgeben (y[0] ist x, y[1] ist xPunkt) + //Oder: + //printf(" %20.5le\t%20.5le\n", t, abs(y[0]-exactH4(w0,t))); //Vergleich von numerischer und analytischer Lösung für x + + RK4(h, w0, y, f, dgl, t, KeineKraft); //RK4 ohne anregende Kraft + } + + //H5: + + printf("t\tLösung von H.5: x(t)\tdx(t)/dt\n"); + for (t = t0; t <= tend; t += h) { + printf(" %20.5le\t%20.5le\t%20.5le\n", t, y[0], y[1]); //Ergebnisse in Tabelle ausgeben + //Oder: + //printf(" %20.5le\t%20.5le\n", t, y[0]- exactH5(w0,t)); //Vgl numerische und analytische Lösung + + + + RK4(h, w0, y, f, dgl, t, ExterneKraft); //Rk4 mit Anregung, die in ExterneKraft gespeichert ist + } + + //H.6: + + printf("Bitte geben sie die Periodenanzahl an:\n"); + scanf("%d", &n); + + for (double w0_scan = 1; w0_scan < 3.2; w0_scan += 0.2) //stellt den scan für w0 dar + { + Leistung_Mittelwert(n, h, w0_scan, y, f, dgl, t0, fext, Leistung); + } + + + //Speicherplatz frei machen + free(y); + free(f); + +}