4. Brauche dringend Hilfe bei implizitem Eulerverfahren

Brauche dringend Hilfe bei implizitem Eulerverfahren

  • ?

    Hi, hoffe Ihr könnt mir weiter helfen. Ich muss ein implizites Eulerverfahren in c-Programmieren und komme nicht ganz weiter.

    Aufganbenstellung ist:
    Es soll ein Programm erstellt werden, das eine gewöhnliche DGL numerisch integriert.

    Zum Testen: y´=t/(t+y)^2 ; Anfanfgsbedingung: y(t0=1,0)=y0=0,0
    Gesucht ist: y(tE=5,0)

    Es soll das implizite Eulerverfahren programmiert werden:
    yi+1=yi+h*f(ti+1,yi+1) ; Anfangsbedingung: y(t0)=y0

    Zur Lösung der nichtlinearen algebrarischen Gleischung ist das Newton- Verfahren zu verwenden.

    Das Eulerverfahren haben ich schon geschieben, aber ich weiss nicht ganz wie ich es impliziet gestalte:

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

    int main()
    {
    /* f=ti/(ti+yi)^2
    Anfangsbedingung:y(t0=1.0)=y0=0.0
    Schrittweite:h=2
    */

    double y, t, K1, h;
    int i, schrittzahl, endwert, schluss;

    printf("Bitte Anzahl Schritte eingeben\n");
    scanf("%d", &schrittzahl);

    printf("Bitte Endwert eingeben\n");
    scanf("%d", &endwert);

    i = 0;
    t = 1;
    y = 0;
    h = (endwert - t) / schrittzahl;
    printf("i = %d, t = %f, y = %f, h = %f", i, t, y, h);

    while (i != schrittzahl) {
    K1 = h * t/((t+y)*(t+y));
    y = y + K1;
    ++i;
    t = t + h;
    }

    printf("\ny = %f", y);
    printf("\nSchluss?\n");
    scanf("%d", &schluss);
    return 0;
    }

    und das Newton-Verfahren dazu kapiere ich leider auch nicht. Würde mich echt freunen wenn Ihr mir weiter helfen könnt. 🙂

    0
  • Mod

    Hallo Terob.,

    wie du siehst ist die Resonanz auf deine Frage nicht gerade groß. Dies liegt vermutlich an drei Dingen, die du nicht beachtet hast, weil du dich weder mit allgemeiner Netiquette auskennst, noch die als wichtig markierten Threads gelesen hast (die dir die Netiquette erklärt hätten).

    -Unleserlicher Beitrag, da du keine Formatierungstags benutzt hast
    -Keine klare Frage. Bloß ein bisschen Code und "geht nicht"/"komm nicht weiter".
    -Kleinigkeit, die eventuell das Fass zum Überlaufen bringt: Auch noch die Behauptung, dein Problem wäre "dringend". Es ist nicht dringend. Es ist dir vielleicht wichtig, das ist aber ganz was anderes. Dringend wäre es, wenn du in einem Atomkraftwerk arbeitest, der letzte Überlebende bist und ganz schnell eine Lösung brauchst, sonst fliegt der ganze Laden in die Luft.

    Siehe dieser Thread (und den dortigen Link!):
    http://www.c-plusplus.net/forum/136013
    Dieser auch:
    http://www.c-plusplus.net/forum/310201

    Wenn du den Beitrag nochmal vernünftig verfasst, dann findet sich wahrscheinlich auch Hilfe. Und ein netter Moderator (ich) ändert vielleicht auch deinen Eingangsbeitrag und die Überschrift, da du das als Unregistrierter nicht kannst, so dass Leser nicht vom ersten Eindruck abgeschreckt werden.

    0
  • W

    Terob. schrieb:

    Aufganbenstellung ist:
    Es soll ein Programm erstellt werden, das eine gewöhnliche DGL numerisch integriert.

    Zum Testen: y´=t/(t+y)^2 ; Anfanfgsbedingung: y(t0=1,0)=y0=0,0
    Gesucht ist: y(tE=5,0)

    Es soll das implizite Eulerverfahren programmiert werden:
    yi+1=yi+h*f(ti+1,yi+1) ; Anfangsbedingung: y(t0)=y0

    Zur Lösung der nichtlinearen algebrarischen Gleischung ist das Newton- Verfahren zu verwenden.

    Du beginnst damit in die Gleichung für das implizite Eulerverfahren 'Dein' f(t_{i+1},y_{i+1}) einzusetzen.
    Es gilt y'=f(t_{i+1},y_{i+1})=t_{i+1}/(t_{i+1}+y_{i+1})^2 also
    y_{i+1} = y_i+h*f(t_{i+1},y_{i+1}) = y_i+h*t_{i+1}/(t_{i+1}+y_{i+1})^2
    und das ist eine Gleichung in der alles bis auf y_{i+1} bekannt ist (t_{i+1}=t_i+h). Man kann versuchen das nach y_{i+1} aufzulösen:
    0 = y_i - y_{i+1} + ht_{i+1}/(t_{i+1}+y_{i+1})^2
    0 = (y_i-y_{i+1}) * (t_{i+1}+y_{i+1})^2 + h    t_{i+1}
    das wird ein Polynom 3.Grades

    Wenn man dies mit Newton-Verfahren lösen möchte, setzt man
    g(y_{i+1}) = (y_i-y_{i+1}) * (t_{i+1}+y_{i+1})^2 + h*t_{i+1}
    und sucht die Nullstelle, für die gilt
    g(y_{i+1}) = 0
    dazu bestimmt man g'(y_{i+1})=dg/d(y_{i+1})
    und iteriert sich mit
    y_{i+1} := y_{i+1} - g(y_{i+1})/g'(y_{i+1})
    an die Nullstelle heran. Und diese ist dann der Wert für y_{i+1}.

    Die Lösung ist y(t=5)=1,09694 (mit Runge-Kutta berechnet). Mit einem h von 0,1 müsstest Du mit dem obigen Verfahren ca. bei 1,085 raus kommen.

    SeppJ hat Recht. Besser wäre es gewesen im Mathe-Forum zu posten. Dass ich Deine Frage gesehen habe, ist purer Zufall.

    Gruß
    Werner

    0
  • ?

    Also erst mal Danke für die schnelle Antwort. Ich versuche die Raatschäge in Zukunft zu berücksichtigen, da es mir schon wichtig ist ;). Ich hoffe für mein ersten Post habe ich eure Gemühter nicht zu sehr beansprucht. Was wäre denn eine passende Überschrift, bin ja für jede Hilfe dankbar.

    Handschriftlich haben wir das alles schon berechnet. Hatte es auch gezielt hier gepostet, da ich Probleme mit dem Quelltext habe und nicht weiß wie ich es implizit zu gestalten habe, dass das Programm es umsetzt.

    0
  • ?

    "Implizit" ist keine Eigenschaft des Codes, sondern des Verfahrens. Es geht darum, das implizite Eulerverfahren zu verwenden, im Gegensatz zum expliziten Eulerverfahren, das dein Code derzeit verwendet.

    Im impliziten Verfahren wird statt

    yk+1 = yk + h * f(tk, yk)

    die Gleichung

    yk+1 = yk + h * f(tk+1, yk+1)

    als Näherung verwendet. tk+1 ist bekannt, aber yk+1 ist es noch nicht -- du versuchst es gerade auszurechnen. Aus diesem Grund brauchst du ein numerisches Suchverfahren (wie das Newton-Verfahren), um yk+1 zu finden.

    0
  • W

    Terob. schrieb:

    Handschriftlich haben wir das alles schon berechnet. Hatte es auch gezielt hier gepostet, da ich Probleme mit dem Quelltext habe und nicht weiß wie ich es implizit zu gestalten habe, dass das Programm es umsetzt.

    Es ist unklar, wo Dein eigentliches Problem liegt. Wie es 'implizit zu gestalten' ist meine ich Dir oben beschrieben zu haben. Wie man es in Code gießt ist was anderes, aber in diesem Fall nicht wirklich schwierig, wenn Du das, was oben steht, allein geschafft hast, wo ist denn das Problem?

    Ich mache mal eine Vorlage:

    #include <stdio.h>
#include <stdlib.h>
#include <math.h> /* fabs */

double newton( double y, double y_i, double t_ip1, double h, double eps );

int main()
{
    /* f=ti/(ti+yi)^2
    Anfangsbedingung:y(t0=1.0)=y0=0.0
    Schrittweite:h=2
    */

    double y, t, K1, h;
    int i, schrittzahl, endwert, schluss;

    printf("Bitte Anzahl Schritte eingeben\n");
    scanf("%d", &schrittzahl);

    printf("Bitte Endwert eingeben\n");
    scanf("%d", &endwert);

    i = 0;
    t = 1;
    y = 0;
    h = (endwert - t) / schrittzahl;
    printf("i = %d, t = %f, y = %f, h = %f\n", i, t, y, h);

    while (i != schrittzahl) {
        /* K1 = h * t/((t+y)*(t+y));
        y = y + K1; */
        y = newton( y, y, t+h, h, 1.E-9 ); /* y alias y_i wird als Startwert für die Näherung gewählt */
        ++i;
        t = t + h;
    }
    printf("i = %d, t = %f, y = %f, h = %f\n", i, t, y, h);

    printf("\nSchluss?\n");
    scanf("%d", &schluss);
    return 0;
}

    wie Du siehst habe ich in Deinem Code lediglich die beiden Zeilen in denen das nächste y berechnet wird, durch einen Aufruf einer Funktion newton() ersetzt.
    In der newton-Funktion benötigst Du noch die Funktion g(y_{i+1}) und deren Ableitung gs(y_{i+1}).
    Erstere sieht so aus, genau wie ich es oben schon beschrieben habe:

    double g( double y_ip1, double y_i, double t_ip1, double h )
{   /* g(y_{i+1}) = (y_i-y_{i+1}) * (t_{i+1}+y_{i+1})^2 + h*t_{i+1} */

    double tmp;
    tmp = t_ip1+y_ip1;
    return (y_i-y_ip1) * tmp * tmp + h*t_ip1;
}

    Die Signatur von gs(y_{i+1}) ist genauso (bitte allein implementieren)

    double gs( double y_ip1, double y_i, double t_ip1, double h );

    Fehlt dann nur noch die Implementieren der Funktion newton() , die g() und gs() aufruft.

    Gruß
    Werner

    0
  • ?

    Erst mal ein dickes Dankeschön an Werner, ich sehe du verstehst was ich meine. 🙂 Ich hätte vorher vielleicht sagen sollen, dass ich nicht wirklich die Ahnung vom Programieren habe. Mir hat bei dem bisherigen ein Freund geholfen, der mit dem jetztigem Problem aber auch nichts mehr anfangen konnte. Heißt ich weiß auch nicht wie ich das implizite in das Program gieße. Sorry für meine undeutliche Frage.

    Handschriftlich komme ich auf die folgende Ableitung:

    F´(t)=-h/(t_{i+1}+y_{i+1})^2 + 2*h*t_{i+1}/(t_{i+1}+y_{i+1})^3

    Grundform war:

    F(t)= 0 = y_{i+1}- y_i - h * t_{i+1}/(t_{i+1}+y_{i+1})^2

    Gruß Terob.

    0
  • W

    Ja - freut mich, dass ich Dir helfen konnte (konnte ich das?)

    Weißt Du nun alles, um die Aufgabe zu lösen? Du stellst eindeutig zu wenig Fragen. Weißt Du doch aus der Sesamstraße: "Wer nicht fragt bleibt dumm." 😉

    Zu Deiner Gleichung: Du darfst nicht nach t sondern musst nach y_{i+1} ableiten. Der Wert für t_{i+1} ist doch bekannt (t_{i+1}=t_i + h).
    Lies Dir bitte seldons und meine Postings nochmal genau durch - steht eigentlich alles drin.

    Kleine Bemerkung, die ich mir nicht verkneifen kann: Wenn es so dringend ist, könntest Du auch schneller auf Antworten reagieren. Umso mehr Feedback Du gibst, umso mehr kann Dir geholfen werden.

    Gruß
    Werner

    0
  • ?

    So, habe es jetzt mal impementiert, kriege aber in der Ausgabe nicht mein gewünschtes Ergebniss von 1,097... . Ich habe keine Ahnung, wie ich, beim Aufruf der Funktionen aus der Funktion newton(), jeweils die ersten Parameter der Funktion g() und gs() verwende.

    Was bringt mir der Parameter eps?
    Habe ich die Funktion gs() korrekt implementiert
    (du sagtest ja, dass sie genau so umgesetzt werden muss, wie es bereits in g() geschah)

    Ist es möglicherweise ein mathematische Problem?
    Stehe total auf dem Schlauch!

    Dein Rückgabewert der Funktion g() habe ich ebenfalls abgeändert, da ich diese bei mir nicht so wiederfinden kann.

    Ist meine Implementierung von newton() denn korrekt?

    Besten Dank für die Mühen und nochmals Sorry! 😞

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h> /* fabs */

    double g( double y_ip1, double y_i, double t_ip1, double h )
    { /* g(y_{i+1}) = (y_i-y_{i+1}) * (t_{i+1}+y_{i+1})^2 + h*t_{i+1} */

    double tmp;
    tmp = t_ip1+y_ip1;
    return (y_ip1-y_i) -(h * t_ip1/( tmp * tmp));
    }

    double gs( double y_ip1, double y_i, double t_ip1, double h ){
    /*double tmp;
    tmp = t_ip1 + y_ip1;
    return ( 1- ((2*h) * t_ip1) / (tmp * tmp * tmp));
    /
    double tmp;
    tmp = t_ip1+y_ip1;
    return (y_i-y_ip1) * tmp * tmp + h    t_ip1;
    }

    double newton( double y, double y_i, double t_ip1, double h, double eps ){
    return y = y_i - (g(y_i, y_i, t_ip1, h) / gs(y_i , y_i, t_ip1, h));
    }
    int main()
    {
    /* f=ti/(ti+yi)^2
    Anfangsbedingung:y(t0=1.0)=y0=0.0
    Schrittweite:h=2
    */

    double y, t, K1, h;
    int i, schrittzahl, endwert, schluss;

    printf("Bitte Anzahl Schritte eingeben\n");
    scanf("%d", &schrittzahl);

    printf("Bitte Endwert eingeben\n");
    scanf("%d", &endwert);

    i = 0;
    t = 1;
    y = 0;
    h = (endwert - t) / schrittzahl;
    printf("i = %d, t = %f, y = %f, h = %f\n", i, t, y, h);

    while (i != schrittzahl) {
    /* K1 = h * t/((t+y)*(t+y));
    y = y + K1; /
    y = newton( y, y, t+h, h, 1.E-9 ); /     y alias y_i wird als Startwert für die Näherung gewählt */
    ++i;
    t = t + h;
    }
    printf("i = %d, t = %f, y = %f, h = %f\n", i, t, y, h);

    printf("\nSchluss?\n");
    scanf("%d", &schluss);
    return 0;
    }

    0
  • W

    Terob. schrieb:

    So, habe es jetzt mal impementiert, kriege aber in der Ausgabe nicht mein gewünschtes Ergebniss von 1,097... . Ich habe keine Ahnung, wie ich, beim Aufruf der Funktionen aus der Funktion newton(), jeweils die ersten Parameter der Funktion g() und gs() verwende.

    Was bringt mir der Parameter eps?

    Es scheint mir, dass Du nicht weißt wie das Newtonverfahren funktioniert. Schau Dir bitte auf der Wikipediaseite zum Newton-Verfahren die erste Graphik an, die ist recht anschaulich und das Beispiel ist auch hilfreich.

    Terob. schrieb:

    Habe ich die Funktion gs() korrekt implementiert

    Nein - schon das g() stimmt nicht. Vergleiche mal Dein g() mit dem g() aus meinem Beitrag. Du hast aus einem '+' ein '-' gemacht. Die Ableitung im Kommentar scheint zu stimmen mit Ausnahme des '-'. Der Kommentar stimmt aber nicht mit dem Code überein!

    Ich habe jetzt leider keine Zeit mehr und bin voraussichtlich bis Montag offline.

    Gruß
    Werner

    0
  • W

    Werner Salomon schrieb:

    Nein - schon das g() stimmt nicht.

    Das muss ich korrigieren. Die Funktion für g() ist korrekt. Da Du die Gleichung gegenüber Deinem ersten Ansatz noch mit -1 multipliziert hast, habe ich das nicht gleich erkannt.

    g(yi+1) = yi+1 - yi - h*ti+1/(ti+1 + yi+1)^2

    double g( double y_ip1, double y_i, double t_ip1, double h ) {
    double tmp;
    tmp = t_ip1+y_ip1;
    return (y_ip1-y_i) - (h * t_ip1/( tmp * tmp)); 
}

    Dann ist g'(yi+1) = 1 + 2h*ti+1/(ti+1 + yi+1)^3

    bzw.in Code

    double gs( double y_ip1, double y_i, double t_ip1, double h ) {
    double tmp;
    tmp = t_ip1+y_ip1;
    return 1 + ((2*h) * t_ip1) / (tmp * tmp * tmp);
}

    also bis auf das Vorzeichen, wie in dem Kommentar zu Deinem Entwurf.

    Damit funktioniert es dann auch. Falls Du Schwierigkeiten bei der Codierung des Newton-Verfahrens hast, so melde Dich bitte und versuche zu schildern worin die Schwierigkeit besteht.

    Gruß
    Werner

    0
Anmelden zum Antworten