klassisches Runge-Kutta-Verfahren

Wally

Erstmal hallöchen und nen schönen Abend,

ich hab hier folgendes explizites Runge-Kutta-Verfahren als C-Code:

const double dt_by_2 = 0.5 * dt;
  const double dt_by_6 = (1./6.) * dt;

  // determine acceleration at time t
  double v_0     = v;                 // slope of x at time t
  double accel_0 = a( x, v, t , dt);  
//der LCP braucht hier dt (und später dt/4,dt/2,dt) weil er nicht zeitdiskret
//arbeitet. Das ist aber ok und funzt so.                   

  // at time t + 0.5 dt using information just obtained
  double x_1     = x + dt_by_2*v_0;
  double v_1     = v + dt_by_2*accel_0;
  double accel_1 = a( x_1, v_1, t + dt_by_2 , dt/4);

  // again at time t + 0.5 dt, but now with the information at the midpoint
  double x_2     = x + dt_by_2*v_1;
  double v_2     = v + dt_by_2*accel_1;
  double accel_2 = a( x_2, v_2, t + dt_by_2 , dt/2);

  // and at the end of the time interval
  double x_3     = x + dt*v_2;
  double v_3     = v + dt*accel_2;
  double accel_3 = a( x_3, v_3, t + dt, dt );

  // the final update is performed using a weighted mean of all the above
  // weighting the center points more than those at the extreme points of 
  // the interval.
  x += dt_by_6*(v_0+v_3+2.*(v_1+v_2));
  v += dt_by_6*(accel_0+accel_3+2.*(accel_1+accel_2));

Das funktioniert soweit auch ganz gut, nur ein Problem hab ich natürlich.
Die DGL a() berechnet bei mir durch lösen linearer Gleichungssysteme eine "Zwangsbeschleunigung". Diese wird auch richtig berechnet! Nur die Position stimmt hinterher nicht, die müsste dann nämlich x_i+1=x_i+dt*v_i+1 sein, stattdessen nähert sie sich dem geforderten Wert wie durch ein PT-1 Glied verzögert.

Wäre es sinnvoll x_i+1 implizit aus v_i+1 zu gewinnen, oder eher Quark (da ich durch a() ja auch "nicht Zwangsbeschleunigungen", also gewöhnliche Bewegungsgleichungen berechne) und diese dann nicht mehr vernünftig integreiert würden?

Wenns sinnvoll ist, wie implementiere ich das dann am besten?

Danke schon mal im voraus , und denunziert mich nicht zu sehr, wegen meiner vielleicht blöden Frage habe nichts Studiert oder so (ist nur Hobby) und spät ist ja auch schon.

PaulM

Hi,

was sind diese "Zwangsbeschleunigungen"? Versuchst du damit x auf einer bestimmten Bahn zu halten?

Es sieht für mich aus, als würdest du eine Differentialgleichung 2. Ordnung mit dem Runge-Kutta-Verfahren versuchen zu lösen. Das geht nur, wenn man sie in eine DGL 1. Ordnung (aber dafür zweidimensional) umschreibt.

im moment hast du:

x' = v
v' = a(x,v,t)

für runge-kutta brauchst du:

s(t) = (x,v) // Orts-Geschwindigkeits-Vektor
s'(t) = F(t,s(t)) = (v,a(x,v,t))

wenn du jetzt die Formel einsetzt, kriegst du

k1 = F(t[i],s[i]) = (v[i],a(x[i],v[i],t[i]))
k2 = F(t + dt/2,s[i] + dt/2*k1) = (x[i] + dt/2*v[i], v[i] + dt/2*a(...))
...
s[i+1] = durchschnitt_aller_ki-s(k1,k2,k3,k4)

Die Rechnung für v(t) ist dieselbe wie oben, aber x[i+1] berechnest du aus x[i] auch mit dem Runge-Kutta-Verfahren, und nicht mit der brutalen euler-methode.

Kann sein, dass ich dich jetzt falsch verstanden habe, ist wirklich schon spät.

Wally

Ok, ich habe die Implementierung mal aufgeräumt und die Zwangskräfte mal beiseite gelassen. Viel wichtiger ist nämlich das Verhalten bei "normalen" mechanischen Anfangswertproblemen.

const double dt_by_2 = 0.5 * dt;
  const double dt_by_6 = (1./6.) * dt;

  double v_0     = v;                 
  double accel_0 = a( x, v);  //dt wird nicht benötigt, da a() zeitdiskret ist

  double x_1     = x + dt_by_2*v_0;
  double v_1     = v + dt_by_2*accel_0;
  double accel_1 = a( x_1, v_1);

  double x_2     = x + dt_by_2*v_1;
  double v_2     = v + dt_by_2*accel_1;
  double accel_2 = a( x_2, v_2);

  // and at the end of the time interval
  double x_3     = x + dt*v_2;
  double v_3     = v + dt*accel_2;
  double accel_3 = a( x_3, v_3);

  v_neu = v + dt_by_6*(accel_0+accel_3+2.*(accel_1+accel_2));
  x_neu = x + dt_by_6*(v_0+v_3+2.*(v_1+v_2));

:das sieht ja für mich auch genauso aus, wie das Verfahren von PaulM, oder bin ich zu blöd dafür?

Wenn a() jetzt eine gefederte Masse unter Erdschwerkraft beschreibt:

double a(double x, double v)
{
    double federkraft = 0.0;
    if (x < 1.0) federkraft = (1.0-federkraft)*FEDERHAERTE;
    return federkraft-9.81*MASSE;
}

dann wächst die Energie im System fast wie beim expliziten Euler-Verfahren.

Wie kann ich das nun vermeiden?

eine Idee wäre ja die letzte Zeile im Verfahren durch:

x_neu = x + dt*v;

zu ersetzen, nur glaub ich da nicht dran, dass das die Lösung sein soll.

PaulM

Wally schrieb:

:das sieht ja für mich auch genauso aus, wie das Verfahren von PaulM, oder bin ich zu blöd dafür?

sorry, stimmt, ich hab' mich im quelltext nicht zurechtgefunden

Wenn a() jetzt eine gefederte Masse unter Erdschwerkraft beschreibt:
double a(double x, double v)
{
    double federkraft = 0.0;
    if (x < 1.0) federkraft = (1.0-federkraft)*FEDERHAERTE;
    return federkraft-9.81*MASSE;
}
dann wächst die Energie im System fast wie beim expliziten Euler-Verfahren.

sollte das nicht (1.0 - x) sein? soweit ich das sehe, braucht RK auch die Differenzierbarkeit von a , und das ist bei diesem Potenzial nicht gegeben. Hast du es mal mit a = -x probiert?

eine Idee wäre ja die letzte Zeile im Verfahren durch:
x_neu = x + dt*v;
zu ersetzen, nur glaub ich da nicht dran, dass das die Lösung sein soll.

ich denke auch, dass das nur durch zufall etwas bringen könnte.