Anwendung der Runge-Kutta-Methode auf ODEs zweiter Ordnung

11

Wie kann ich die Euler-Methode durch Runge-Kutta 4. Ordnung ersetzen, um die Bewegung des freien Falls in nicht konstanter Gravitationsgröße (z. B. freier Fall aus 10 000 km Höhe) zu bestimmen?

Bisher habe ich eine einfache Integration nach der Euler-Methode geschrieben:

while()
{
    v += getMagnitude(x) * dt;
    x += v * dt;
    time += dt;
}

x Variable bedeutet aktuelle Position, v bedeutet Geschwindigkeit, getMagnitude (x) gibt die Beschleunigung auf x Position zurück.

Ich habe versucht, RK4 zu implementieren:

while()
{
    v += rk4(x, dt) * dt; // rk4() instead of getMagintude()
    x += v * dt;
    time += dt;
}

Dabei ist der Funktionskörper von rk4 ():

inline double rk4(double tx, double tdt)
{
   double k1 = getMagnitude(tx);
   double k2 = getMagnitude(tx + 0.5 * tdt * k1);
   double k3 = getMagnitude(tx + 0.5 * tdt * k2);
   double k4 = getMagnitude(tx + tdt * k3);

   return (k1 + 2*k2 + 2*k3 + k4)/6.0;
}

Aber etwas stimmt nicht, weil ich nur einmal mit RK4 (Beschleunigung) integriere. Die Integration der Geschwindigkeit mit RK4 ist nicht sinnvoll, da sie mit v * dt identisch ist.

Können Sie mir sagen, wie man Differentialgleichungen zweiter Ordnung mithilfe der Runge-Kutta-Integration löst? Sollte ich RK4 implementieren, indem ich die Koeffizienten k1, l1, k2, l2 ... l4 berechne? Wie kann ich das machen?

Marcin W.
quelle
Hallo @Marcin, ich habe Ihren Titel bearbeitet, um besser zu reflektieren, was Ihrer Meinung nach Ihr Problem tatsächlich ist. Ich denke, wir erhalten möglicherweise nützlichere Antworten und es wird für andere, die diese Frage in Zukunft mit dem neuen Titel sehen, besser durchsuchbar sein. Fühlen Sie sich frei, es zurück zu ändern, wenn Sie nicht einverstanden sind.
Doug Lipinski

Antworten:

17

Es scheint ziemlich verwirrend zu sein, wie mehrstufige (z. B. Runge-Kutta) Methoden auf ODEs 2. oder höherer Ordnung oder ODE-Systeme angewendet werden sollen. Der Prozess ist sehr einfach, wenn Sie ihn erst einmal verstanden haben, aber ohne eine gute Erklärung vielleicht nicht offensichtlich. Die folgende Methode finde ich am einfachsten.

F=mx¨

[x˙v˙]=[vF/m]

vxk1k4(x,v)

while (t<TMAX)
    k1 = RHS( t, X );
    k2 = RHS( t + dt / 2, X + dt / 2 * k1 );
    k3 = RHS( t + dt / 2, X + dt / 2 * k2 );
    k4 = RHS( t + dt, X + dt * k3 );
    X = X + dt / 6 * ( k1 + 2 * k2 + 2 * k3 + k4 );
    t = t + dt;
end

X=(x,v)RHS( t, X )(x˙(t),v˙(t))

Leider unterstützt C ++ solche Vektoroperationen nicht nativ, sodass Sie entweder eine Vektorbibliothek verwenden, Schleifen verwenden oder die einzelnen Teile manuell ausschreiben müssen. In C ++ können Sie std::valarrayden gleichen Effekt erzielen. Hier ist ein einfaches Arbeitsbeispiel mit konstanter Beschleunigung.

#include <valarray>
#include <iostream>

const size_t NDIM = 2;

typedef std::valarray<double> Vector;

Vector RHS( const double t, const Vector X )
{
  // Right hand side of the ODE to solve, in this case:
  // d/dt(x) = v;
  // d/dt(v) = 1;
  Vector output(NDIM);
  output[0] = X[1];
  output[1] = 1;
  return output;
}

int main()
{

  //initialize values

  // State variable X is [position, velocity]
  double init[] = { 0., 0. };
  Vector X( init, NDIM );

  double t = 0.;
  double tMax=5.;
  double dt = 0.1;

  //time loop
  int nSteps = round( ( tMax - t ) / dt );
  for (int stepNumber = 1; stepNumber<=nSteps; ++stepNumber)
  {

    Vector k1 = RHS( t, X );
    Vector k2 = RHS( t + dt / 2.0,  X + dt / 2.0 * k1 );
    Vector k3 = RHS( t + dt / 2.0, X + dt / 2.0 * k2 );
    Vector k4 = RHS( t + dt, X + dt * k3 );

    X += dt / 6.0 * ( k1 + 2.0 * k2 + 2.0 * k3 + k4 );
    t += dt;
  }
  std::cout<<"Final time: "<<t<<std::endl;
  std::cout<<"Final position: "<<X[0]<<std::endl;
  std::cout<<"Final velocity: "<<X[1]<<std::endl;

}
Doug Lipinski
quelle
6
" Leider unterstützt C ++ solche Vektoroperationen nicht nativ " Ich denke, dass dies in der Standardbibliothek sogar der Fall ist, aber nicht unbedingt einfach mit anderen linearen Algebra-Bibliotheken zu verwenden: en.cppreference.com/w/cpp/numeric/valarray, denke ich Gängige lineare Algebra-Bibliotheken wie Eigen sollten ebenfalls als "Unterstützung" gelten.
Kirill
1
@Kirill, Danke für den Tipp. Ich bin noch relativ neu in C ++ und habe Valarray noch nie benutzt. Ich habe gerade auch etwas Nützliches gelernt! Bearbeiten zum Hinzufügen.
Doug Lipinski
1
Vielleicht ist dieser Rat auch dann hilfreich: 1) Verwenden Sie das Clang-Format zum automatischen Formatieren Ihres Codes. Es ist wirklich Standard und einheitlich. 2) Verwendung typedef std::valarray<double> Vectorfür häufig verwendete Typen. 3) Verwenden Sie const int NDIM = 2anstelle von #definefür Typensicherheit und Korrektheit. 4) Seit C ++ 11 können Sie den Körper von RHS einfach durch ersetzen return {X[1], 1}. 5) In C ++ ist es wirklich ungewöhnlich (im Gegensatz zu C), zuerst Variablen zu deklarieren und später zu initialisieren. Variablen lieber an derselben Stelle zu deklarieren, an der Sie sie initialisieren ( double t = 0.usw.)
Kirill
1
@MarcinW. RHS()berechnet die rechte Seite der Differentialgleichung. Der Zustandsvektor X ist (x, v), also ist dX / dt = (dx / dt, dv / dt) = (v, a). Für Ihr Problem (wenn a = G * M / x ^ 2) sollte RHS zurückkehren { X[1], G*M/(X[0]*X[0]) }.
Doug Lipinski
1
@Kirill Ich weiß, aber das funktioniert nur seit C ++ 11, was bedeutet, dass es nicht mit Standard-Compileroptionen auf den beliebtesten Compilern funktioniert. Ich habe mich dafür entschieden, dies zugunsten von etwas wegzulassen, das auch mit den alten Standards funktioniert, und hoffentlich die Verwirrung zu verringern, die durch die Unfähigkeit verursacht wird, den Code zu kompilieren.
Doug Lipinski