#include #include #include int mmax = 30; //globale Variable, Maximale Anzahl der Schritte double f1(double x, double a, double z) // Zu integrierende Funktion { return (exp(-(x * x) / (a * a)) / sqrt(x * x + z * z)); } double f2(double x, double a, double z) { return((exp(-(x*x)/(a*a))*z)/pow((x*x)+(z*z),1.5)); } double trapez(int n, double a, double z, double (*func)(double, double, double), double x1, double x2) { //Trapeze double sum; double h = (x2 - x1) / (n); //Trapezsumme für Grenzen berechnen sum = (0.5 * func(x1, a, z)) + (0.5 * func(x2, a, z)); //in der Summe Funtion aufrufen um Werte an Grenzen x1 und x2 zu bekommen for (int i = 1; i < n; i++) { sum += func(x1 + (h * i), a, z); } return sum * h; } double IntRomb(double x1, double x2, int n, double a, double z, double (*func)(double, double, double), int eps) { int i, k; //Variablen für Schleifen double ergebnis; //zum speichern des Ergebnis double **T; //Trapezsummen (später als 2-dim Array) double *h; //Schrittweiten Array h = (double*) malloc((mmax + 1) * sizeof(double)); //Speicherplatz für Array allozieren T = (double**) malloc(((mmax + 1)) * sizeof(double*)); for (i = 0; i <= mmax; i++) { T[i] = (double*) malloc(((mmax + 1)) * sizeof(double)); //Speicherplatz für 2-dim Array allozieren } h[0] = (x2 - x1) / (n); //Erste Schrittweite berechnen T[0][0] = trapez(n, a, z, func, x1, x2); //erstes Trapez berechnen for (i = 1; i <= mmax; i++) { n = n * 2; //Trapeze für k=0 berechnen, indem man n vergrößert T[i][0] = trapez(n, a, z, func, x1, x2); h[i] = h[i - 1] / 2.0; //h berechnen für i>1 indem man durch 2 teilt for (k = 1; k <= i; k++) { //mit Formel aus VL T für k>1 berechnen T[i][k] = (-h[i - k] * h[i - k] / ((h[i] * h[i]) - (h[i - k] * h[i - k])) * T[i][k - 1]) + (h[i] * h[i] / (h[i] * h[i] - h[i - k] * h[i - k]) * T[i - 1][k - 1]); } if (abs(T[i][i] - T[i - 1][i - 1]) <= eps) { //Abschätzen ob das Ergebnis genau genug ist wenn ja aus Schleife raus ergebnis = T[i][i]; break; } } if (i==mmax){ printf("Fehler, Integral ist ungenau"); // Fehler falls Abbruchbedingung der vorherigen Schleife nicht greift } //Speicherplatz befreien for (i = 0; i <= mmax; i++) { //Speicherplatz in Schleife freigeben free(T[i]); } free(T); free(h); return ergebnis; //Rückgabewert ist Ergebnis } double InfInt(int n, double a, double z, double eps,double (*func)(double, double, double)) { double x1 = -1; //Grenzen definieren double x2 = 1; double R=IntRomb(x1,x2,n,a,z,func,eps); //Variable zum speichern des Ergebnis, Integral für erste Grenzen ausrechnen double lastR; //Vergleichsvariable definieren do { lastR = R;//letztes Ergebnis speichern x1 = x1*2; //Grenzen vergrößern x2 = x2*2; R = IntRomb(x1, x2, n, a, z, func, eps); //neues Ergebnis berechnen } while ((lastR/R<(1-eps)) || (lastR/R>(1+eps))); //Integralwerte vergleichen ,sodass deren Änderung kleiner als eps ist return R; } double potenzial(int n, double a, double z, double Q, double (*func)(double, double, double), double eps) { double V; //Variable Potential double alpha = 1/sqrt(3.141); //normierung alpha V = alpha * (Q / a) * InfInt(n, a, z, eps,func); //Potential berechnen return V; } double ges_potenzial(int n, double a1, double a2, double z, double Q1, double Q2, double d, double (*func)(double, double, double), double eps) { /*Berechnung des gesamten Potenzials der beiden Drähte durch Superposition, wobei für die z-Koordinate, bzw den Abstand der Ladung zum zweiten Drahtes z-d gilt.*/ double V_ges; V_ges = potenzial(n,a1,z,Q1,f1,eps) + potenzial(n,a2,z-d,Q2,f1,eps); return V_ges; } void efeld(double a, double z, int n, double Q, double eps, double *dfp,double (*func)(double, double, double)) { double f[2]; //Array zum speichern der Funktionswerte double h =0.5* z; //h relativ groß wird in Schleife verkleinert double lastdf; f[0] = potenzial( n, a, z +(h), Q, func, eps); //f(z+h) berechnen f[1] = potenzial(n, a, z - h, Q, func, eps); //f(z-h) berechnen *dfp = -1* (f[0] - f[1]) / (double)(2 * h); do { h=h/2; //h kleiner machen lastdf=*dfp; //Alte Ableitung auf neuer Variable speichern f[0] = potenzial( n, a, z +(h), Q, f1, eps); //f(z+h) berechnen f[1] = potenzial(n, a, z - h, Q, f1, eps); //f(z-h) berechnen *dfp = -1* (f[0] - f[1]) / (2 * h); //Differenzenquotient } while((lastdf/(*dfp)<(1-eps)) || (lastdf/(*dfp)>(1+eps))); //Werte vergleichen ,sodass deren Änderung kleiner als eps ist } double ges_efeld(int n, double a1, double a2, double z, double Q1, double Q2, double d, double *df1, double *df2, double (*func)(double, double, double), double eps) { //verwendet numerisch berechnetes E-Feld efeld(a1, z, n, Q1, eps, df1, func); //Berechnung E-Feld des 1. Drahtes efeld(a2, z-d , n, Q2, eps, df2, func); //Berechnung E-Feld des 2. Drahtes double E_ges = *df1 + *df2; //Superposition der E-Felder der einzelnen Drähte return E_ges; } double ges_efeld2(int n, double a1, double a2, double z, double Q1, double Q2, double d, double (*func)(double, double, double), double eps) { //verwendet analytisch berechnetes E-Feld double E_ges2; //Superposition der E-Felder der einzelnen Drähte E_ges2 = potenzial(n, a1, z, Q1, func, eps) + potenzial(n, a2, z-d, Q2, func, eps); return E_ges2; } double secant(int n, double z1, double z2, double a1, double a2, double Q1, double Q2, double d, double *df1, double *df2, double (*func)(double, double, double), double eps, int *schritt) { //verwendet numerisch berechnetes E-Feld //Bestimmung des Glechgewichtspuntes bzw.der Nullstelle des addierten E-Feldes mit den Startwerten z1, z2 double zn; //Variable zum Speichern des neu berechneten z-Wertes *schritt = 0; while (fabs(z2-z1) > eps) //solange bis gewünschte Genauigkeit erreicht wurde { //Berechnung des nächsten Schätzwertes und speichern dieses Wertes als neuen wert zn = z2 - ges_efeld(n, a1, a2, z2, Q1, Q2, d, df1, df2, func, eps) * (z2-z1)/ (ges_efeld(n, a1, a2, z2, Q1, Q2, d, df1, df2, func, eps) - ges_efeld(n, a1, a2, z1, Q1, Q2, d, df1, df2, func, eps)); z1 = z2; z2 = zn; (*schritt)++; //Schrittezähler } return zn; } double secant2(int n, double z1, double z2, double a1, double a2, double Q1, double Q2, double d, double (*func)(double, double, double), double eps, int *schritt) { //verwendet analytisch berechnetes E-Feld //Bestimmung des Glechgewichtspuntes bzw.der Nullstelle des addierten E-Feldes mit den Startwerten z1, z2 double zn; //Variable zum Speichern des neu berechneten z-Wertes *schritt = 0; while (fabs(z2-z1) > eps) //solange bis gewünschte Genauigkeit erreicht wurde { //Berechnung des nächsten Schätzwertes und speichern dieses Wertes als neuen wert zn = z2 - ges_efeld2(n, a1, a2, z2, Q1, Q2, d, func, eps) * (z2-z1)/ (ges_efeld2(n, a1, a2, z2, Q1, Q2, d, func, eps) - ges_efeld2(n, a1, a2, z1, Q1, Q2, d, func, eps)); z1 = z2; z2 = zn; (*schritt)++; //Schrittezähler } return zn; } int main(void) { double z, z1, z2, a1, a2; int n = 30; const double eps = 1e-8; //genauigkeit epsilon double d = 12; //Abstand 2.Draht zum 1. double Q1 = 1; //Gesamtladung der Drähte 1 und 2 double Q2 = 4; double df1; //Variable zum speichern des Differenzenquotienten double df2; //Eingabe der Variablen: printf("Bitte geben Sie den Abstand der Punktladung von der x-Achse z ein\n "); scanf("%lf", &z); printf("Bitte geben Sie a1 ein (0 < a1 <= 4):\n"); //Eingabe der Breiten der Drähte scanf("%lf", &a1); printf("Bitte geben Sie a2 ein (a1/10 <= a2 <= a1):\n"); scanf("%lf", &a2); printf("Gesamtes Potenzial der beiden Drähte: %f\n", ges_potenzial(n, a1, a2, z, Q1, Q2, d, f1, eps)); printf("Gesamtes E-Feld (numerisch berechnet) durch Addition ist: %f\n", ges_efeld(n, a1, a2, z, Q1, Q2, d, &df1, &df2, f1, eps)); printf("Gesamtes E-Feld (analytisch berechnet) durch Addition ist: %f\n", ges_efeld2(n, a1, a2, z, Q1, Q2, d, f2, eps)); //H3.2 Nullstellenberechnung mit Sekantenverfahren int schritt; //Variable zum speichern der Schrittanzahl des Sekantenverfahrens printf("Bitte geben Sie die Startwerte z1, z2 für die Bestimmung der Gleichgewichtslage an:\n"); scanf("%lf %lf", &z1, &z2); /*Numerische Methode dauert aber sehr lange printf("Gleichgewichtspunkt des Teilchens: %1.8f\n", secant(n, z1, z2, a1, a2, Q1, Q2, d, &df1, &df2, f1, eps, &schritt)); */ //analytische Methode /*printf("Gleichgewichtspunkt des Teilchens: %1.8f\n", secant2(n, z1, z2, a1, a2, Q1, Q2, d, f2, eps, &schritt)); printf("Schrittanzahl:%d\n",schritt);*/ return 0; }