Hausaufgaben/H3.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

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;
}