Application de la méthode Runge-Kutta aux ODE de second ordre

11

Comment puis-je remplacer la méthode Euler par Runge-Kutta 4e ordre pour déterminer le mouvement de chute libre dans une amplitude gravitationnelle non constante (par exemple, chute libre à partir de 10 000 km au-dessus du sol)?

Jusqu'à présent, j'ai écrit une intégration simple par la méthode Euler:

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

x variable signifie la position actuelle, v signifie la vitesse, getMagnitude (x) renvoie l'accélération sur la position x.

J'ai essayé d'implémenter RK4:

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

où le corps de la fonction rk4 () est:

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;
}

Mais quelque chose ne va pas, car je n'intègre qu'une seule fois en utilisant RK4 (accélération). L'intégration de la vitesse à l'aide de RK4 n'a pas de sens car c'est la même chose que v * dt.

Pourriez-vous me dire comment résoudre des équations différentielles de second ordre en utilisant l'intégration Runge-Kutta? Dois-je implémenter RK4 en calculant les coefficients k1, l1, k2, l2 ... l4? Comment puis je faire ça?

Marcin W.
la source
Bonjour @Marcin, j'ai modifié votre titre pour mieux refléter ce que je pense que votre problème est réellement. Je pense que nous pourrons obtenir des réponses plus utiles et il sera plus consultable pour d'autres qui verront cette question à l'avenir avec le nouveau titre. N'hésitez pas à le modifier en cas de désaccord.
Doug Lipinski

Réponses:

17

Il semble y avoir un peu de confusion sur la façon d'appliquer des méthodes à plusieurs étapes (par exemple Runge-Kutta) à des ODE d'ordre 2 ou supérieur ou à des systèmes d'ODE. Le processus est très simple une fois que vous le comprenez, mais peut-être pas évident sans une bonne explication. La méthode suivante est celle que je trouve la plus simple.

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))

Malheureusement, C ++ ne prend pas en charge nativement les opérations vectorielles comme celle-ci, vous devez donc utiliser une bibliothèque vectorielle, utiliser des boucles ou écrire manuellement les parties distinctes. En C ++, vous pouvez utiliser std::valarraypour obtenir le même effet. Voici un exemple de travail simple avec une accélération constante.

#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
la source
6
" Malheureusement, C ++ ne prend pas en charge nativement les opérations vectorielles comme celle-ci " Je pense que oui, même dans la bibliothèque standard, mais pas nécessairement facile à utiliser avec d'autres bibliothèques d'algèbre linéaire: en.cppreference.com/w/cpp/numeric/valarray je pense les bibliothèques d'algèbre linéaire communes comme Eigen, devraient également compter comme "support".
Kirill
1
@Kirill, merci pour le conseil. Je suis encore relativement nouveau en C ++ et je n'ai jamais utilisé valarray auparavant, je viens d'apprendre quelque chose d'utile aussi! Modification à ajouter.
Doug Lipinski
1
Peut-être que ces conseils seront également utiles: 1) Utilisez le format clang pour formater automatiquement votre code, c'est vraiment standard et uniforme. 2) Utilisez typedef std::valarray<double> Vectorpour les types couramment utilisés. 3) Utiliser const int NDIM = 2au lieu de #definepour la sécurité et l'exactitude du type. 4) Depuis C ++ 11, vous pouvez simplement remplacer le corps de RHS par return {X[1], 1}. 5) Il est vraiment rare en C ++ (contrairement à C) de déclarer d'abord des variables, puis de les initialiser plus tard, préférez déclarer des variables au même endroit où vous les initialisez ( double t = 0., etc.)
Kirill
1
@MarcinW. RHS()calcule le côté droit de l'équation différentielle. Le vecteur d'état X est (x, v) donc dX / dt = (dx / dt, dv / dt) = (v, a). Pour votre problème (si a = G * M / x ^ 2) RHS devrait revenir { X[1], G*M/(X[0]*X[0]) }.
Doug Lipinski
1
@Kirill Je sais, mais cela ne fonctionne que depuis C ++ 11, ce qui signifie qu'il ne fonctionne pas avec les options de compilation par défaut sur les compilateurs les plus populaires. J'ai choisi de laisser cela de côté en faveur de quelque chose qui fonctionne également avec les anciennes normes et, espérons-le, de réduire la confusion causée par l'incapacité de compiler le code.
Doug Lipinski