MATLAB-Code Hilfe. Rückwärts-Euler-Verfahren

StackOverflow https://stackoverflow.com/questions/2937183

Frage

Hier ist die MATLAB / FreeMat Code ich habe eine ODE numerisch die rückwärts-Euler Methode. Allerdings sind die Ergebnisse nicht mit meinem Lehrbuch Ergebnissen, und manchmal sogar lächerlich inkonsistent. Was ist falsch mit dem Code? 

function [x,y] = backEuler(f,xinit,yinit,xfinal,h)

    %f - this is your y prime
    %xinit - initial X
    %yinit - initial Y
    %xfinal - final X
    %h - step size

    n = (xfinal-xinit)/h; %Calculate steps

    %Inititialize arrays...
    %The first elements take xinit and yinit corespondigly, the rest fill with 0s.
    x = [xinit zeros(1,n)];
    y = [yinit zeros(1,n)];

    %Numeric routine
    for i = 1:n
        x(i+1) = x(i)+h;
        ynew = y(i)+h*(f(x(i),y(i)));
        y(i+1) = y(i)+h*f(x(i+1),ynew);
    end
end
War es hilfreich?

Lösung

Ihre Methode ist eine Methode, einer neuen Art . Es ist weder rückwärts noch vorwärts Euler. : -)

Vorwärts-Euler: y1 = y0 + h*f(x0,y0)

Rückwärts-Euler solve in y1: y1 - h*f(x1,y1) = y0

Ihre Methode: y1 = y0 +h*f(x0,x0+h*f(x0,y0))

Ihre Methode ist nicht Rückwärts-Euler.

Sie lösen nicht in y1, die Sie gerade schätzen y1 mit dem Vorwärts Euler-Verfahren. Ich möchte nicht, die eine Analyse der Methode verfolgen, aber ich glaube, es ist in der Tat schlecht verhalten wird, auch im Vergleich mit nach vorn Euler, da Sie die Funktion f am falschen Punkt zu bewerten.

Hier ist die nächste Methode, um Ihre Methode, dass ich mir vorstellen kann, explizit als auch, was viel bessere Ergebnisse geben sollte. Es ist Heun-Verfahren :

y1 = y0 + h/2*(f(x0,y0) + f(x1,x0+h*f(x0,y0)))

Andere Tipps

Die Frage kann ich nur erkennen ist, dass die Zeile: 

n=(xfinal-xinit)/h

Sollte sein: 

n = abs((xfinal-xinit)/h)

Um eine negative Schrittzahl zu vermeiden. Wenn Sie in der negativen x-Richtung bewegen, stellen Sie sicher, dass der Funktion eine negative Schrittgröße zu geben.

Ihre Antworten wahrscheinlich wegen abweichen, wie grob Sie Ihre Antwort sind annähert. Um ein halb genaues Ergebnis zu erhalten, muss deltaX sehr als sehr klein und Ihre Schrittgröße muss sehr sehr sehr klein ist.

PS. Dies ist nicht die „Rückwärts-Euler-Methode,“ es ist nur regelmäßige alte Euler-Methode ist.

Wenn dies Hausaufgaben bitte Tag es so.

Hier finden Sie aktuelle numerischen Rezepte , insbesondere Kapitel 16, Integration gewöhnlichen Differential Gleichungen. Euler-Verfahren ist bekannt, Probleme haben:

Es gibt mehrere Gründe dafür, dass die Eulersche Methode nicht für den praktischen Gebrauch wird empfohlen, unter ihnen, (i) das Verfahren ist nicht sehr genau, wenn im Vergleich zu anderen, schicker, laufen Verfahren bei der äquivalenten Schrittgrßenregister, und (ii) weder ist es sehr stabil

Also, wenn Sie wissen, dass Ihr Lehrbuch Euler-Methode verwendet, würde ich nicht die Ergebnisse Spiel erwarten. Auch wenn es, haben Sie wahrscheinlich eine identische Schrittgröße zu erhalten ein identisches Ergebnis verwenden.

Wenn Sie wirklich wollen eine ODE über Euler-Methode zu lösen, dass Sie selbst geschrieben haben, sollten Sie einen Blick auf Einbau-ODE-Solver .

Auf einer Bemerkung am Rande: Sie brauchen nicht x(i) innerhalb der Schleife wie diese zu erstellen: x(i+1) = x(i)+h;. Stattdessen können Sie einfach x = xinit:h:xfinal; schreiben. Außerdem können Sie n = round(xfinal-xinit)/h); zu vermeiden Warnungen schreiben.

Hier sind die Solver von MATLAB implementiert.

ode45 basiert auf einer expliziten Runge-Kutta (4,5) Formel, die Dormand-Prince Paar. Es ist ein einstufiges Solver - bei der Berechnung von y (tn), es braucht nur die Lösung bei der unmittelbar vorhergehenden Zeitpunkt, y (tn-1). Im Allgemein ode45 ist die beste Funktion gilt als ein erster Versuch für die meisten Probleme.

ode23 ist eine Implementierung eines explizites Runge-Kutta (2,3) Paar Bogacki und Shampine. Es kann sein, mehr effizienter als ode45 bei Rohöl Toleranzen und in Gegenwart von moderate Steifigkeit. Wie ode45, ode23 ist ein einstufiges Löser.

ode113 ist eine Variable, um Adams-Bashforth-Moulton PECE Löser. Es kann effizienter sein als ode45 bei strenge Toleranzen und wenn die ODE Datei-Funktion ist besonders teuer zu bewerten. ode113 ist ein Multistep-Löser - es normalerweise Bedürfnisse die Lösungen bei mehreren vorhergehenden Zeitpunkt um den Strom zu berechnen, Lösung.

Die obigen Algorithmen sollen lösen steife Systeme. Wenn sie erscheinen mäßig langsam zu sein, versuchen Sie eine der Verwendung die steifen Solver unten.

ode15s ist eine Variable, um Solver auf der Grundlage der numerischen Differentiation Formeln (NDF). Optional verwendet es das BDF-Verfahren (BDF, das auch als Zahnrad Methode bekannt) Das ist in der Regel weniger effizient. Mögen ode113, ode15s ist ein mehrstufiger Löser. Versuchen Sie ode15s wenn ode45 ausfällt, oder ist sehr ineffizient, und Sie vermuten, dass das Problem ist steif, oder wenn die Lösung ein differentiell-algebraische Problem.

ode23s basiert auf einer modifizierten Rosenbrock- Formel der Ordnung 2. Weil es ist ein Ein-Schritt-Solver, kann es sein, effizienter als ode15s bei Rohöl Toleranzen. Es kann einige Arten von lösen steife Probleme, für die nicht ode15s Wirksam.

ode23t ist eine Implementierung der Trapezregel eines „freien“ mit Interpolant. Verwenden Sie diesen Solver, wenn die Problem ist nur mäßig steif und Sie brauchen eine Lösung, ohne numerische Dämpfung. ode23t kann DAEs lösen.

ode23tb ist eine Implementierung von TR-BDF2, ein impliziter Runge-Kutta Formel mit einer ersten Stufe, daß a Trapezregel und ein zweiter Schritt Bühne, die eine rückständig Differenzierung Formel der Ordnung zwei. Durch die Konstruktion ist die gleiche Iteration Matrix wird bei der Bewertung der beiden verwendet Stufen. Wie ode23s kann dieser Solver effizienter sein als ode15s bei Rohöl Toleranzen.

Ich denke, dieser Code funktionieren könnte. Versuchen Sie dies. 

for i =1:n
    t(i +1)=t(i )+dt;
    y(i+1)=solve('y(i+1)=y(i)+dt*f(t(i+1),y(i+1)');
    end

Der Code ist in Ordnung. Nur Sie eine weitere Schleife innerhalb der for-Schleife hinzufügen. Um das Niveau der Konsistenz zu überprüfen. 

if abs((y(i+1) - ynew)/ynew) > 0.0000000001 
    ynew = y(i+1);
    y(i+1) = y(i)+h*f(x(i+1),ynew);
end

Ich habe für eine Dummy-Funktion und die Ergebnisse waren vielversprechend.

Lizenziert unter: CC-BY-SA mit Zuschreibung
Nicht verbunden mit StackOverflow
scroll top