Dateien hochladen nach „“

hausaufgabe2.c

@@ -0,0 +1,247 @@
+//gcc hausaufgabe2.c -o hausaufgabe2 -lm -lh2_cip -L.
+//./hausaufgabe2
+//Pauline Maas, Ariane Ufer
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#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);
+
+}