Implizites RKV < Matlab < Mathe-Software < Mathe < Vorhilfe
|
Status: |
(Frage) beantwortet | Datum: | 09:48 Mi 01.07.2009 | Autor: | jojo2502 |
Hallo,
für ein Simulationsprojekt muss ich ein implizites Runge-Kutta-Verfahren 4. Ordung implementieren. Leider finde ich im Netz keine für mich verständlichen Ansätze zur Implementierung. Hat hier vielleicht jemand ein Beispiel oder brauchbare Links zum Thema steife DGL´s und impliziten Runge-Kutta-Verfahren. Bei bedarf kann ich auch gerne weitere Informationen zur Aufgabenstellung geben.
Gruß Jojo
P.S: Leider haben mir die Antworten in anderen Foren noch nicht weitergeholfen.
Ich habe diese Frage auch in folgenden Foren auf anderen Internetseiten gestellt:
http://www.matheplanet.com/matheplanet/nuke/html/viewtopic.php?topic=125244
http://www.matheplanet.com/matheplanet/nuke/html/viewtopic.php?topic=125334
|
|
|
|
Hallo
> Hallo,
>
> für ein Simulationsprojekt muss ich ein implizites
> Runge-Kutta-Verfahren 4. Ordung implementieren. Leider
> finde ich im Netz keine für mich verständlichen Ansätze
> zur Implementierung. Hat hier vielleicht jemand ein
> Beispiel oder brauchbare Links zum Thema steife DGL´s und
> impliziten Runge-Kutta-Verfahren. Bei bedarf kann ich auch
> gerne weitere Informationen zur Aufgabenstellung geben.
>
> Gruß Jojo
>
> P.S: Leider haben mir die Antworten in anderen Foren noch
> nicht weitergeholfen.
> Ich habe diese Frage auch in folgenden Foren auf anderen
> Internetseiten gestellt:
>
> http://www.matheplanet.com/matheplanet/nuke/html/viewtopic.php?topic=125244
>
> http://www.matheplanet.com/matheplanet/nuke/html/viewtopic.php?topic=125334
Ich habe mal ein Runge-Kutta-Verfahren implementiert. Hier der Matlab-Code:
function [t,y] = RungeKutta (f ,a ,b ,ya ,h)
%f Funktion, a und b Grenzen, ya Anfangswert, h Schrittgrösse
%Anfangswerte und Laufvariabel
t = [a];
y = [ya];
i = 1;
while t(i) + h <= b
i = i + 1;
t(i) = t(i-1) + h;
%Definition der ki's
k1 = h * f(t(i-1),y(i-1));
k2 = h * f(t(i-1) + h / 2,y(i-1) + k1 / 2);
k3 = h * f(t(i-1) + h / 2,y(i-1) + k2 / 2);
k4 = h * f(t(i-1) + h, y(i-1) + k3);
%Runge-Kutta Schritt
y(i) = y(i-1) + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6;
end
Hoffe, das hilft dir weiter..
Grüsse, Amaro
|
|
|
|
|
Hallo,
vielen Dank für die schnell Rückmeldung.
Leider ist der von dir gepostete Algorithmus das explizite Runge-Kutta-Verfahren 4. Ordnung, welches sich leider nicht zum lösen steifer DGL´s eignet. Also bin ich immer noch auf der Suche nach Beispielen.
Trotzdem vielen Dank Amaro.
Gruß Jojo
|
|
|
|
|
> Hallo,
> vielen Dank für die schnell Rückmeldung.
> Leider ist der von dir gepostete Algorithmus das explizite
> Runge-Kutta-Verfahren 4. Ordnung, welches sich leider nicht
> zum lösen steifer DGL´s eignet. Also bin ich immer noch
> auf der Suche nach Beispielen.
> Trotzdem vielen Dank Amaro.
>
> Gruß Jojo
Hallo Jojo,
könntest du für uns mal kurz erklären:
1.) Was ist der wesentliche Unterschied
zwischen "explizitem" und "implizitem" RKV ?
2.) Was versteht man unter einer "steifen" DGL ?
... und gib doch bitte ein (einfaches) Beispiel einer
DGL von der Sorte an, die du lösen willst !
Gruß Al-Chw.
|
|
|
|
|
Status: |
(Mitteilung) Reaktion unnötig | Datum: | 12:37 Mi 01.07.2009 | Autor: | jojo2502 |
1.) Bei impliziten RK-Verfahren k¨onnen die Stufen ki nicht einfach
sukzessiv bestimmt werden, sondern sind als Lösung eines (i.a.
nichtlinearen) Gleichungssystems definiert.
2.) Steife Probleme tauchen bei Systemen von Differentialgleichungen auf,
deren Jacobi-Matrix durchweg Eigenwerte mit negativem Realteilen
aufweisen, die darüberhinaus noch von stark unterschiedlicher
Größenordnung sind.
Beide Erklärungen sind aus dem Netz...
Bsp.:
[mm]\ddot{x} = \bruch{1}{m}*(P*A - c(x)*x-h(\dot{x})-F_{AN})[/mm]
Um das Beispiel zu lösen muss dieses ja in ein System von Gleichungen 1. Ordnung umgewandelt werden.
Gruß Jojo
|
|
|
|
|
> 1.) Bei impliziten RK-Verfahren können die Stufen [mm] k_i
[/mm]
> nicht einfach sukzessiv bestimmt werden, sondern
> sind als Lösung eines (i.a. nichtlinearen) Gleichungs-
> systems definiert.
>
> 2.) Steife Probleme tauchen bei Systemen von
> Differentialgleichungen auf, deren Jacobi-Matrix
> durchweg Eigenwerte mit negativen Realteilen
> aufweisen, die darüberhinaus noch von stark
> unterschiedlicher Größenordnung sind.
>
> Bsp.:
>
> [mm]\ddot{x} = \bruch{1}{m}*(P*A - c(x)*x-h(\dot{x})-F_{AN})[/mm]
>
> Um das Beispiel zu lösen muss dieses ja in ein System von
> Gleichungen 1. Ordnung umgewandelt werden.
>
> Gruß Jojo
Hi Jojo,
Danke für die Erläuterungen. Weshalb "steife" DGLn
nach dieser Definition nicht gut mit dem "normalen"
RK lösbar sind, habe ich zwar noch nicht geschnallt,
aber lassen wir dies mal.
Auf einer Webseite über Runge-Kutta
http://en.wikipedia.org/wiki/Runge–Kutta_methods
habe ich folgende Aussage gefunden:
>>> The simplest example of an implicit Runge–Kutta
>>> method is the backward Euler method:
>>> $\ [mm] y_{n+1}\ [/mm] =\ [mm] y_n\,+\,h*f(t_n+h,y_{n+1})$
[/mm]
Ich würde es mal damit versuchen. Um ein Prinzip
zu begreifen, sind meistens die einfachsten Reali-
sierungen die nützlichsten, weil da die Rechnungen
(hoffentlich) noch einigermaßen überschaubar bleiben.
Um aus der gegebenen DGL 2.Ordnung ein
System von DGLn 1.Ordnung zu machen,
können wir definieren:
$\ [mm] y(t):=\vektor{u(t)\\v(t)}$
[/mm]
mit $\ u(t)=x(t)$ und $\ [mm] v(t)=\dot{x}(t)$
[/mm]
Das System sieht dann so aus:
[mm] $\begin{cases} (1) & \dot{u}(t)=v(t) \\ (2) & \dot{v}(t)= \bruch{1}{m}*(P*A - u*c(u)-h(v)-F_{AN}) \end{cases}$
[/mm]
Ich versuche mal (ohne Gewähr !), daraus die
impliziten RK- Gleichungen abzuleiten:
[mm] $\begin{cases} (1) & u_{n+1}=u_n+h*v_{n+1} \\ (2) & v_{n+1}=v_n+ \bruch{h}{m}*(P*A - u_{n+1}*c(u_{n+1})-h(v_{n+1})-F_{AN}) \end{cases}$
[/mm]
Ich wäre froh, wenn da jemand genau drüber
schauen würde ...
Heute habe ich sonst nicht mit (instabilen oder
stabileren) Lösungsverfahren zu kämpfen, aber
leider mit einer hoffnungslos instabilen ADSL-
Verbindung. Es ist ein Übel, wie das Netz mehr
und mehr mit einer unheimlichen Menge von
Datenmüll überschwemmt wird. Dabei nehme
ich einmal an, dass das was wir hier schreiben,
kein Müll ist und auch quantitativ kaum
groß ins Gewicht fällt ...
LG Al-Chw.
|
|
|
|
|
Status: |
(Mitteilung) Reaktion unnötig | Datum: | 15:07 Mi 01.07.2009 | Autor: | Arcesius |
> >>> The simplest example of an implicit Runge–Kutta
> >>> method is the backward Euler method:
>
> >>> [mm]\ y_{n+1}\ =\ y_n\,+\,h*f(t_n+h,y_{n+1})[/mm]
Weswegen ich an einen Eulerschritt gedacht habe, um die neuen ki's zu berechnen. Allerdings bräuchte man dafür einen Euler-Forward-Schritt.. hmm...
|
|
|
|
|
Status: |
(Mitteilung) Reaktion unnötig | Datum: | 15:28 Mi 01.07.2009 | Autor: | jojo2502 |
Hallo Al-Chwarizmi
erstmal ein großes Dankeschön an dich für deine Mühen.
Ich glaube ich habe eine Lösung für mein Problem gefunden und denke auch das ich es soweit verstehe. Mit der Vorschrift für implizite Runge-Kutta-Verfahren:
[mm]
k_{r} := f(t_{j} + \alpha_{r} h, u_{j} + h(\beta_{r1} k_{1} + \cdots + \beta_{rm} k{m}))
für r = 1, \cdots, m
u_{j+1} := u_{j} + h(\gamma_{1} k_{1} + \cdots + y_{m} k_{m})
[/mm]
und dem Butcher-Diagramm
[Dateianhang nicht öffentlich]
[Dateianhang nicht öffentlich]
komme ich auf eine Matlabimplementierung. Dies ist die Lobato IIIB-Form
Die Lösungen für die [mm]k_{i}[/mm] werden zur Zeit mit einem einfachen Fixpunktiterationsverfahren gefunden.
1: |
| 2: | function [y,t]=rungekutta_implizit(f,ta,ya,te,N,M)
| 3: | h=(te-ta)/N;
| 4: | h6=h/6;
| 5: | t = ta:h:te; % Gitter anlegen
| 6: | y=zeros(length(ya),N+1);
| 7: | y(:,1) = ya; % Startwert schreiben
| 8: |
| 9: | for i=1:N % durchlaufe Gitterpunkte
| 10: | %setze Startwerte fuer Iterationsverfahren
| 11: | k1=zeros(length(ya),1);
| 12: | k2=zeros(length(ya),1);
| 13: | %starte Iterationsverfahren
| 14: | for j=1:M
| 15: | k1alt=k1; % alte Werte merken fuer Iteration
| 16: | k2alt=k2;
| 17: | k1=f(t(i),y(:,i)+h6*(k1alt-k2alt));
| 18: | k2=f(t(i)+h/2,y(:,i)+h6*(k1alt+2*k2alt));
| 19: | end
| 20: | k3=f(t(i)+h,y(:,i)+h6*(k1+5*k2));
| 21: | y(:,i+1)=y(:,i)+h6*(k1+4*k2+k3);
| 22: | end
| 23: | end
|
Dateianhänge: Anhang Nr. 1 (Typ: png) [nicht öffentlich] Anhang Nr. 2 (Typ: png) [nicht öffentlich]
|
|
|
|
|
> Hallo Al-Chwarizmi
>
> erstmal ein großes Dankeschön an dich für deine Mühen.
>
> Ich glaube ich habe eine Lösung für mein Problem gefunden
> und denke auch das ich es soweit verstehe.
> .....
> .....
Hallo Jojo, Arcesius
und alle anderen,
jetzt hätte ich eigentlich den Wunsch, die
Funktionsweise eines impliziten Verfahrens
auch zu verstehen - obwohl mein Ansatz
mit dem "reluE", also dem backwards-
Euler nicht vierter Ordnung ist. Ich denke
allerdings nicht, dass dabei auch ein Euler-
Vorwärtsschritt involviert ist.
Mit Butcher-Diagrammen habe ich mich
noch nie beschäftigt und will das auch nicht
unbedingt.
Könnte jemand meine in diesem Artikel
dargestellte Idee überprüfen ?
LG Al
|
|
|
|
|
Hallo
> Hallo Jojo, Arcesius
> und alle anderen,
>
> jetzt hätte ich eigentlich den Wunsch, die
> Funktionsweise eines impliziten Verfahrens
> auch zu verstehen - obwohl mein Ansatz
> mit dem "reluE", also dem backwards-
> Euler nicht vierter Ordnung ist. Ich denke
> allerdings nicht, dass dabei auch ein Euler-
> Vorwärtsschritt involviert ist.
Doch doch.. schaue weiter unten, ich habe ein Code geschrieben, da wird es erläutert:
> Mit Butcher-Diagrammen habe ich mich
> noch nie beschäftigt und will das auch nicht
> unbedingt.
Auch noch nie was davon gehört...
>
> Könnte jemand meine in
> diesem Artikel
> dargestellte Idee überprüfen ?
>
> LG Al
>
>
Das Problem ist, dass die Eulermethoden und das explizite RKV NUR für DGL erster Ordnung gut sind... Ich habe mich bisher nicht mit weiteren numerischen Methoden für DGL zweiter Ordnung auseinandergesetzt...
Aber hier mal den Backward-Euler:
function [t,y] = BackwardEuler(f,a,b,ya,h)
%Anfangswerte und Laufvariabel.
t = [a];
y = [ya];
i = 1;
while t(i) + h <= b
i = i+1;
t(i) = t(i-1) + h;
%Forward-Euler-Schritt
q = y(i-1) + h*f(t(i-1),y(i-1));
%Backward-Euler-Schritt
y(i) = y(i-1) + h*f(t(i),q);
end
(Die Formatierfarben gehen grade nicht.. naja ^^)
Die explizite Forward-Euler-Methode nimmt die Steigung am Punkt y und geht so zum Punkt y+1.
Die implizite Backward-Euler-Methode nimmt die Steigung am Punkt y+1, geht dann zurück zum Punkt y und von dort mit der Steigung von y+1 zum nächsten Punkt. D.h, die Funktion muss zuerst einen Schritt vorwärts machen (Euler Forward), und dann wieder einen zurück (Euler Backward).
Also sind implizite Methoden immer vom nächsten Schritt abhängig.
Mein Beitrag ist jetzt vielleicht nicht das, was du hören wolltest, aber leider ist mein Wissen in Numerik im Moment noch ein bisschen begränzt.. :)
Vielleicht kann jemand anders auf deine Idee eingehen!
Grüsse, Amaro
|
|
|
|
|
> Hallo
>
>
> > Hallo Jojo, Arcesius
> > und alle anderen,
> >
> > jetzt hätte ich eigentlich den Wunsch, die
> > Funktionsweise eines impliziten Verfahrens
> > auch zu verstehen - obwohl mein Ansatz
> > mit dem "reluE", also dem backwards-
> > Euler nicht vierter Ordnung ist. Ich denke
> > allerdings nicht, dass dabei auch ein Euler-
> > Vorwärtsschritt involviert ist.
>
> Doch doch.. schaue weiter unten, ich habe ein Code
> geschrieben, da wird es erläutert:
> >
> > Könnte jemand meine in
> > diesem Artikel
> > dargestellte Idee überprüfen ?
> >
>
> Das Problem ist, dass die Eulermethoden und das explizite
> RKV NUR für DGL erster Ordnung gut sind... Ich habe mich
> bisher nicht mit weiteren numerischen Methoden für DGL
> zweiter Ordnung auseinandergesetzt...
Nun, man kann ja eben die DGL 2.Ord. auf
ein System erster Ordnung für einen
Zweiervektor zurückführen, wie ich das
dargestellt habe.
> Aber hier mal den Backward-Euler:
>
> function [t,y] = BackwardEuler(f,a,b,ya,h)
>
> %Anfangswerte und Laufvariabel.
> t = [a];
> y = [ya];
>
> i = 1;
>
> while t(i) + h <= b
> i = i+1;
>
> t(i) = t(i-1) + h;
>
> %Forward-Euler-Schritt
> q = y(i-1) + h*f(t(i-1),y(i-1));
>
> %Backward-Euler-Schritt
> y(i) = y(i-1) + h*f(t(i),q);
> end
>
> (Die Formatierfarben gehen grade nicht.. naja ^^)
>
> Die explizite Forward-Euler-Methode nimmt die Steigung am
> Punkt y und geht so zum Punkt y+1.
>
> Die implizite Backward-Euler-Methode nimmt die Steigung am
> Punkt y+1, geht dann zurück zum Punkt y und von dort mit
> der Steigung von y+1 zum nächsten Punkt. D.h, die Funktion
> muss zuerst einen Schritt vorwärts machen (Euler Forward),
> und dann wieder einen zurück (Euler Backward).
>
> Also sind implizite Methoden immer vom nächsten Schritt
> abhängig.
Ich stell mir diese implizite Methode anders vor:
Wir haben den Punkt [mm] y_{i-1} [/mm] und die Funktion f(t,y).
Nun wird für das neue y, also [mm] y_{i} [/mm] die Gleichung
aufgestellt:
$\ [mm] y_{i}\ [/mm] =\ [mm] y_{i-1}+h*f(t_{i},y_{i})$
[/mm]
Nun hat man diese implizite Gleichung für [mm] y_{i}
[/mm]
nach [mm] y_{i} [/mm] aufzulösen - womöglich mit einer
geeigneten numerischen Methode.
Man berechnet also nicht zuerst einen "provisorischen"
neuen Wert wie etwa bei Euler-Heun.
Ich meine, dass das was du in dem Code dar-
gestellt hast, eine etwas abgewandelte Heun-
Methode ist, und keine "implizite" Methode.
Gruß Al-Chw.
|
|
|
|
|
Status: |
(Mitteilung) Reaktion unnötig | Datum: | 17:47 Mi 01.07.2009 | Autor: | Arcesius |
Hallo Al-Chw> > Hallo
> Ich stell mir diese implizite Methode anders vor:
>
> Wir haben den Punkt [mm]y_{i-1}[/mm] und die Funktion f(t,y).
> Nun wird für das neue y, also [mm]y_{i}[/mm] die Gleichung
> aufgestellt:
>
> [mm]\ y_{i}\ =\ y_{i-1}+h*f(t_{i},y_{i})[/mm]
>
> Nun hat man diese implizite Gleichung für [mm]y_{i}[/mm]
> nach [mm]y_{i}[/mm] aufzulösen - womöglich mit einer
> geeigneten numerischen Methode.
> Man berechnet also nicht zuerst einen "provisorischen"
> neuen Wert wie etwa bei Euler-Heun.
> Ich meine, dass das was du in dem Code dar-
> gestellt hast, eine etwas abgewandelte Heun-
> Methode ist, und keine "implizite" Methode.
>
> Gruß Al-Chw.
>
Ich habe ein schönes Dokument gefunden:
Klicke hier
Es kann sein, dass die implizite Methode nicht dadurch charakterisiert wird, dass sie den nächsten Schritt braucht und dies nur bei der Euler-Backward-Methode halt notwendig ist...
Grüsse, Amaro
|
|
|
|
|
Status: |
(Mitteilung) Reaktion unnötig | Datum: | 09:40 Do 02.07.2009 | Autor: | jojo2502 |
He interessantes Dokument,
das erklärt einiges. Ich merke immer mehr das mein werter Herr Professor
mir ne schöne Aufgabe aufgedrückt hat. Eigentlich nicht im Fach Numerik sonder
in Simulation dynamischer System. Naja werde mich dort durchkämpfen.
Gruß Jojo
|
|
|
|
|
Status: |
(Mitteilung) Reaktion unnötig | Datum: | 17:02 Mi 01.07.2009 | Autor: | jojo2502 |
Hallo,
leider kann ich mich heute nicht mehr wirklich mit der Thematik beschäftigen
und werde mich erst morgen wieder dazu äussern können.
Erst mal euch beiden vielen Dank.
Vielleicht steig ich ja morgen dort durch.
Gruß Jojo
|
|
|
|
|
Status: |
(Mitteilung) Reaktion unnötig | Datum: | 11:42 Mi 01.07.2009 | Autor: | Arcesius |
Hallo
Vielleicht eine dumme Idee, aber wie wäre es, wenn man den nächsten y-Wert mit einem Eulerschritt berechnen würde, um den Wert in den ki's dann anstelle des y(i-1) zu benutzen?
(Es würde die geniale genauigkeit des RKV einfach ein bisschen zurückschrauben...)
Oder sehen die ki's der impliziten Version ganz anders aus?
Grüsse, Amaro
|
|
|
|
|
Status: |
(Mitteilung) Reaktion unnötig | Datum: | 12:43 Mi 01.07.2009 | Autor: | jojo2502 |
Siehe
hier
Gruß Jojo
|
|
|
|
|
Status: |
(Mitteilung) Reaktion unnötig | Datum: | 11:20 Do 16.07.2009 | Autor: | matux |
$MATUXTEXT(ueberfaellige_frage)
|
|
|
|