Runge-Kutta-Verfahren in C++ < Numerik < Hochschule < Mathe < Vorhilfe
|
| Status: |
(Frage) beantwortet  | | Datum: | 20:28 Do 29.01.2015 | | Autor: | binhbob |
| Aufgabe | | Runge-Kutta-Verfahren in C++ implementieren. |
Hallo, ich wollte das Runge-Kutta-Verfahren zur Lösung der Bewegungsdgl. in C++ implementieren.
Bewegungsdgl.: m*y''+d*y'+k*y = p
Umformung in eine DGL-System erster Ordnung.
y = y0
y' = y1
y''= -(d/m)*y1-(k/m)*y0+p/m
Dabei ist der nächste Schritt [mm] y_m+1 [/mm] = [mm] y_m [/mm] +h*y'_m
// c++ Code
//Startwerte
x.at(0)=x0; //x=y0
v.at(0)=v0; //v=y1
a.at(0)=-(d/m)*v.at(0)-(k/m)*x.at(0)+(p/m); //a=y''
// kann man den Startwert von a so berechnen?
for(unsigned int i=1; i<(n); i++)
{
xa.at(i)=(h/2.0)*v.at(i-1); // Koeffizieten für Runge-Kutta
xb.at(i)=(h/2.0)*xa.at(i);
xc.at(i)=(h)*xb.at(i);
x.at(i)=x.at(i-1)+((h/6.0)*(v.at(i-1)+2.0*xa.at(i)+2.0*xb.at(i)+xc.at(i)));
va.at(i)=(h/2.0)*a.at(i-1); // Koeffizieten für Runge-Kutta
vb.at(i)=(h/2.0)*va.at(i);
vc.at(i)=(h)*vb.at(i);
v.at(i)=v.at(i-1)+((h/6.0)*(a.at(i-1)+2.0*va.at(i)+2.0*vb.at(i)+vc.at(i)));
a.at(i)=-(d/m)*v.at(i)-(k/m)*x.at(i)+(p/m); // Hier ist der Knackpunkt. Kann ich a so direkt berechnen? Oder muss ich auf a das Runge-Kutta-Verfahren anwenden. So komme ich theoretisch auf den nächsten Schritt von a.
}
Das Programm läuft. Jedoch liefert die numerische Lösung nicht annähernd die analytische Lösung.
Die Probleme:
1.) Die Frequenz der Schwingung ist zu niedrig um genau zu sein ungefähr um den Faktor 6.
2.) Die Amplitude nimmt zu, wenn die Dämpfung(d) gleich Null ist. Jedoch stimmt die Amplitude der Schwigung.
Ich habe diese Frage in keinem Forum auf anderen Internetseiten gestellt.
|
|
| |
|
| Status: |
(Antwort) fertig  | | Datum: | 12:07 Fr 30.01.2015 | | Autor: | leduart |
Hallo
irgendwas ist falsch an deinen Fprmeln
was di am Anfang darüber schreibst ist falsch das ist das einfache Eulerverfahren.
schreib bitte die Runge Kutta Formel hin, die du verwenden willst.
Sicher falsch ist, dass du zu oft mit h multiplizierst
deine Koeffizienten mit h und dann am Ende nochmal mit h/6 [mm] x_b [/mm] hat schon [mm] h^2 [/mm] am Ende dann [mm] h^3!
[/mm]
Vielleicht solltest du um nicht durcheinander zu kommen die Koeffizienten und die x mit verschiedenen Namen belegen
Gruß leduart
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 18:41 Fr 30.01.2015 | | Autor: | binhbob |
Also das Runge-Kutta-Verfahren, dass ich verwenden will ist das klassische Runge-Kutta-Verfahren vierter Ordnung.
Also es stimmt wenn man es genau nimmt ist es $ [mm] y_m+1=ym+h/6\cdot{}(k1+2k2+2k2+k4) [/mm] $
Ich habe grad rübergeguckt und der Algorithmus stimmt einfach nicht..
Dieser Algorithmus hier klappt aufjedenfall für diese Funktion y'=y
k1[i] := y[i-1];
k2[i] := y[i-1]+(1/2)*h*k1[i];
k3[i] := y[i-1]+(1/2)*h*k2[i];
k4[i] := h*k3[i]+y0[i-1];
y[i] := y[i-1]+(1/6)*h*(k1[i]+2*k2[i]+2*k3[i]+k4[i])
Dieser wurde in Maple geschrieben und liefert auch genau die E-Funktion.
Nur wie verknüpfe ich mein Gleichungssystem erster Ordnung miteinander?
Wenn ich jetzt diesen Algorithmus auf meine beiden ersten DGL anwende
(y = y0; y' = y1), dann kriege ich ja für beide Funktion eine E-Funktion? Doch die sollen ja schwingen.
Ich habe dabei noch eine anderen Annahme wir wollen ja den nächsten Schritt berechnen und nähern uns diesen so, indem wir unseren derzeitigen Funktionswert $ [mm] y_m [/mm] $ + der Schrittweite h * der Ableitung $ [mm] y_m' [/mm] $ an der Stelle. Da wir ja unsere Ableitung kennen.
x'=v könnte man doch schlicht für die Ableitung v einsetzen.
So wäre der nächste Schritt $ [mm] y_m+1=y_m+h\cdot{}v_m [/mm] $ oder?
|
|
|
| |
|
| Status: |
(Antwort) fertig | | Datum: | 22:12 Fr 30.01.2015 | | Autor: | leduart |
Hallo
Du hast anscheinend das Runge Kutta verfahren nicht verstanden
was du vorschlägst ist das einfache Eulerverfahren, das aber schnell versagt, probier es mal mit y'=-y y(0=1 sollte den cos liefern aus (auch dein eigenartiges Verfahren, mit deinen k, warum das funktioniert verstehe ich nicht, es kann hüchstens für diese einfachste Dgl funktionieren, du hast in [mm] k_1 [/mm] ja einfach, wenn du die Endformel ansiehst h!6*y(i-1) addiert. mit 2*k2 gehst du dann nochmal 2h/6* (y(i-1) und [mm] h^2/6
[/mm]
mit [mm] k_3 [/mm] hast du h^63 erreich mit [mm] k_4 [/mm] dann [mm] h^4 [/mm] wenn h klein genug ist fällt das wwohl wef.
schreib mal deine k so um, dass du sie nicht ineinander einsetzt, sondern wirklich einsetzt. dann sieht das gräßlich aus. kriegst du wirklich in maple damit DIE e funktion, oder einfach eine so ähnliche? wenn die [mm] h^2 [/mm] klein genug sind spielen sie kaum eine Rolle gegen h und du hast insgesamt auf raffinierte und umständliche Weise einfach yI)=y(i-1)+h*y(i-1) wobei ja y(i-1) auch y'(i-1)ist, also ein Verfahren, das genau deine und nur dein y'=y halbwegs lost. rechne mal mit berschiedenen h y(1) aus, das sollt e ergeben!
jetzt zum richtigen RK
was heisst Ich habe grad rübergeguckt und der Algorithmus stimmt einfach nicht.
worüber? über deinen code oder das richtige Verfahren? das richtige stimmt und wird unglaublich oft benutzt
musst du Runge Kutta? für dgl 2 ter ordnung reicht auch das "Halbschrittverfahren:"
Dgl: y'=f(y) (auch f(t,y) y,f skalar oder Vektoren.
du hast [mm] y_0 [/mm] daraus [mm] y_m=y_0+h/2*F(y) [/mm] d.as heisst du rechnest den Wert nach einem halben Schritt aus indem du mit der Steigung y' das Stück h/2 weiter gehst. mit der Steigung in der Mitte gehst du jetzt vom Anfangspunkt an einen ganzen Schritt weiter. dann hast du [mm] y_1=y_0+h*F(y_m)
[/mm]
beim Runge Kutte rechnet man jetzt die Steigung nicht nur in der Mitte aus, sondern an 3 Zwischenstellen, diese mittelt man gewichtet.
also [mm] y_{m1}=y_0+h/2*F(y_0)
[/mm]
[mm] y_{m2}=y_0+h/2*F(y_(m1)
[/mm]
[mm] y_(m3)=y_o+h*F(y_{m2}
[/mm]
diese 4 Steigungen mittelt man jetzt mit den Gewichten 1,2,2,1
und hat dann [mm] y_1=y_0+h/6*(F((y_{0}+2*F(y:{m1}+2*F(y_{m2}+1*F(y(m_3))
[/mm]
das sollten deine [mm] k_i [/mm] sein Die [mm] k_i [/mm] sind also nichts anderes als Steigungen an den Zwischenstellen.
wird es damit klarer.
du hast eine Vektorfkt
mit deinen Namen [mm] \vektor{x \\ v}'=(\vektor{v\\ f(x,v} [/mm] dein f ist das was du a genannt hast
und musst natürlich immer beide Komponenten parallel rechnen um f an den entsprechenden Stellen zu kennen.
Gruß leduart
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig  | | Datum: | 18:40 Fr 30.01.2015 | | Autor: | binhbob |
Also das Runge-Kutta-Verfahren, dass ich verwenden will ist das klassische Runge-Kutta-Verfahren vierter Ordnung.
Also es stimmt wenn man es genau nimmt ist es [mm] y_m+1=ym+h/6*(k1+2k2+2k2+k4)
[/mm]
Ich habe grad rübergeguckt und der Algorithmus stimmt einfach nicht..
Dieser Algorithmus hier klappt aufjedenfall für diese Funktion y'=y
k1[i] := y[i-1];
k2[i] := y[i-1]+(1/2)*h*k1[i];
k3[i] := y[i-1]+(1/2)*h*k2[i];
k4[i] := h*k3[i]+y0[i-1];
y[i] := y[i-1]+(1/6)*h*(k1[i]+2*k2[i]+2*k3[i]+k4[i])
Dieser wurde in Maple geschrieben und liefert auch genau die E-Funktion.
Nur wie verknüpfe ich mein Gleichungssystem erster Ordnung miteinander?
Wenn ich jetzt diesen Algorithmus auf meine beiden ersten DGL anwende
(y = y0; y' = y1), dann kriege ich ja für beide Funktion eine E-Funktion? Doch die sollen ja schwingen.
Ich habe dabei noch eine anderen Annahme wir wollen ja den nächsten Schritt berechnen und nähern uns diesen so, indem wir unseren derzeitigen Funktionswert [mm] y_m [/mm] + der Schrittweite h * der Ableitung [mm] y_m' [/mm] an der Stelle. Da wir ja unsere Ableitung kennen.
x'=v könnte man doch schlicht für die Ableitung v einsetzen.
So wäre der nächste Schritt [mm] y_m+1=y_m+h*v_m [/mm] oder?
|
|
|
|
