Browse Source

Mit secant

master
maasp 5 years ago
parent
commit
2b9e4d7567
1 changed files with 115 additions and 69 deletions
  1. 184
      Hausaufgaben/ComputerphysikAbgabe/src/ComputerphysikAbgabe.c

184
Hausaufgaben/ComputerphysikAbgabe/src/ComputerphysikAbgabe.c

@ -1,15 +1,18 @@
#include <stdio.h>
#include <stdlib.h>
#include<math.h>
#include <math.h>
int mmax = 30; //globale Variable, Maximale Anzahl der Schritte
double InfInt(int n, double a, double z,double eps,double (*func)(double, double, double)); //InfInt wird bei Potenzial gebraucht deshalb ist es hier definiert
double f1(double x, double a, double z) { // Zu integrierende Funktion
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 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;
@ -34,11 +37,8 @@ double IntRomb(double x1, double x2, int n, double a, double z,
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++) {
@ -52,17 +52,14 @@ double IntRomb(double x1, double x2, int n, double a, double z,
/ ((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]);
}
//printf("%d",i);
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){
if (i == mmax) {
printf("Fehler, Integral ist ungenau"); // Fehler falls Abbruchbedingung der vorherigen Schleife nicht greift
}
//Speicherplatz befreien
@ -74,85 +71,134 @@ double IntRomb(double x1, double x2, int n, double a, double z,
return ergebnis; //Rückgabewert ist Ergebnis
}
double potenzial(int n, double a, double z,
double (*func)(double, double, double), double eps) {
double V; //Variable Potential
double alpha = 1/sqrt(M_PI); //normierung alpha
double Q = 4; //Gesamtladung Q
V = alpha * (Q / a) * InfInt(n, a, z, eps,func); //Potential berechnen
return V;
}
double InfInt(int n, double a, double z, double eps,double (*func)(double, double, double)) {
double InfInt(int n, double a, double z, double eps,
double (*func)(double, double, double)) {
double x1 = -1; //Grenzen definieren
double x2 = 0;
double R=IntRomb(x1,x2,n,a,z,func,eps); //Variable zum speichern des Ergebnis, Integral für erste Grenzen ausrechnen
double x2 = 0; //hier ist x2 null, da das Integral achsensymmetrisch ist
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;
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 2*R;
} while ((lastR / R < (1 - eps)) || (lastR / R > (1 + eps))); //Integralwerte vergleichen ,sodass deren Änderung kleiner als eps ist
return 2 * R; //auf Grund von achsensymmetrie 2 mal
}
void efeld(double a, double z, int n, double eps,
double *dfp,double (*func)(double, double, double)) {
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 h = 0.5 * z; //h relativ groß wird in Schleife verkleinert
double lastdf;
f[0] = potenzial( n, a, z +(h), func, eps); //f(z+h) berechnen
f[1] = potenzial(n, a, z - h, func, eps); //f(z-h) berechnen
*dfp = -1* (f[0] - f[1]) / (double)(2 * h);
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), func, eps); //f(z+h) berechnen
f[1] = potenzial(n, a, z - h, func, eps); //f(z-h) berechnen
h = h / 2; //h kleiner machen
lastdf = *dfp; //Alte Ableitung auf neuer Variable speichern
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]) / (2 * h); //Differenzenquotient
} while ((lastdf / (*dfp) < (1 - eps)) || (lastdf / (*dfp) > (1 + eps))); //Werte vergleichen ,sodass deren Änderung kleiner als eps ist
}
*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, d - z, 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 secant_pauline(int n, double a1, double a2, double Q1, double Q2,
double d, double *df1, double *df2,
double (*func)(double, double, double), double eps) {
//Setze die startgrenzen ein bisschen von den Drähten entfernt, da das Feld auf den Drähten null ist
double g1 = eps;
double g2 = d - eps;
double zero;
//Calculating secants until the zero changes less than eps
while ((g1 / g2) > (1 + eps) || (g1 / g2) < (1 - eps)) {
//Berechnen von mx + c aus Grenzen p1 und p2
double m = (ges_efeld(n, a1, a2, g2, Q1, Q2, d, df1, df2, func, eps)
- ges_efeld(n, a1, a2, g1, Q1, Q2, d, df1, df2, func, eps))
/ (g2 - g1);
double c = ges_efeld(n, a1, a2, g1, Q1, Q2, d, df1, df2, func, eps)
- (g1 * m);
//Mit dem Sekantenverfahren die nullstelle berechnen
zero = (-1) * (c / m);
g2 = g1;
g1 = zero;
}
//returning the result
return zero;
}
int main(void) {
double z, a;
int n=30;
double eps = 1e-8; //genauigkeit epsilon
double z, a1, a2;
int n = 30;
const double eps = 1e-8; //genauigkeit epsilon
const double eps2 = 1e-4;
double d = 12; //Abstand 2.Draht zum 1.
double Q1 = 1; //Gesamtladung der Drähte 1 und 2
double Q2 = 4;
double df;
double df1; //Variable zum speichern des Differenzenquotienten
double df2;
double df; //Variable zum speichern des Differenzenquotienten
//Eingabe der Variablen:
printf("Bitte geben Sie den Abstand der Punktladung von der x-Achse z ein\n ");
printf(
"Bitte geben Sie den Abstand der Punktladung von der x-Achse z ein\n ");
scanf("%lf", &z);
printf("Bitte geben Sie a ein\n");
scanf("%lf", &a);
printf("Unendliches Inegral %lf\n", InfInt(n, a, z, eps,f1)); //Integral berechnen mit InfInt
printf("Das Potenzial ist %f\n", potenzial(n,a,z,f1,eps)); //Potenzial aus Integral berechnen
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);
//H1 berechnung des Potentials
printf("Das Potenzial ist %f\n", potenzial(n, a1, z, Q1, f1, eps)); //Potenzial aus Integral berechnen
//H2.1 Berechnung des Efeldes mit Hilfe numerischer Ableitung
efeld(a,z,n,eps,&df,f1); //E-feld berechnen
efeld(a1, z, n, Q1, eps, &df, f1); //E-feld berechnen
printf("Das E-Feld an der Stelle %3.2lf ist:%f\n", z, df);
//H2.2 Berechnung des Efeldes mit Hilfe analytischer Ableitung
printf("Das E-Feld an der Stelle %3.2lf ist:%f\n", z, potenzial(n,a,z,f2,eps));
double diff=fabs(df - potenzial(n,a,z,f2,eps));
printf("Differenz %1.10f", diff);
printf("Achtung schleife");
for(int i=0; i<20; i++){
efeld(a,z,n,eps,&df,f1);
printf("Das E-Feld an der Stelle %3.2lf ist:%f\n", z, df);
printf("Das E-Feld an der Stelle %3.2lf ist:%f\n", z, potenzial(n,a,z,f2,eps));
double diff=fabs(df - potenzial(n,a,z,f2,eps));
printf("ZZ %3.2lf\n", z);
printf("Differenz %1.10f \n", diff);
z=z*2;
}
printf("Das E-Feld an der Stelle %3.2lf ist:%f\n", z,
potenzial(n, a1, z, Q1, f2, eps));
double diff = fabs(df - potenzial(n, a1, z, Q1, f2, eps));
printf("Differenz %1.10f", diff);
//H3 zwei Drähte: Setzte a1 als 4
a1 = 4;
printf("Gesamtes Potenzial der beiden Drähte: %f\n",
ges_potenzial(n, a1, a2, z, Q1, Q2, d, f1, eps2));
printf("Gesamtes E-Feld (numerisch berechnet) durch Addition ist: %f\n",
ges_efeld(n, a1, a2, z, Q1, Q2, d, &df1, &df2, f1, eps2));
printf("Nullstelle %f\n",
secant_pauline(n, a1, a2, Q1, Q2, d, &df1, &df2, f1, eps2));
z = 0.1;
}
Loading…
Cancel
Save