Aplicación del método Runge-Kutta a las EDO de segundo orden

¿Cómo puedo reemplazar el método de Euler por el cuarto orden de Runge-Kutta para determinar el movimiento de caída libre en una magnitud gravitacional no constante (por ejemplo, caída libre desde 10 000 km sobre el suelo)?

Hasta ahora escribí integración simple por el método de Euler:

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

x variable significa la posición actual, v significa velocidad, getMagnitude (x) devuelve la aceleración en la posición x.

Intenté implementar RK4:

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

donde el cuerpo de la función rk4 () es:

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

Pero algo está mal, porque me estoy integrando solo una vez usando RK4 (aceleración). La velocidad de integración con RK4 no tiene sentido porque es lo mismo que v * dt.

¿Podría decirme cómo resolver ecuaciones diferenciales de segundo orden utilizando la integración Runge-Kutta? ¿Debo implementar RK4 calculando los coeficientes k1, l1, k2, l2 ... l4? ¿Cómo puedo hacer eso?

Parece haber bastante confusión acerca de cómo aplicar los métodos de varios pasos (por ejemplo, Runge-Kutta) a las EDO de segundo orden o superiores o los sistemas de EDO. El proceso es muy simple una vez que lo entiendes, pero quizás no sea obvio sin una buena explicación. El siguiente método es el que encuentro más simple.

$F = m\ddot{x}$

$v$$x$k1k4$(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 )$(\dot{x}(t),\dot{v}(t))$

Desafortunadamente, C ++ no admite de manera nativa operaciones de vectores como esta, por lo que debe usar una biblioteca de vectores, usar bucles o escribir manualmente las partes separadas. En C ++ puedes usar std::valarraypara lograr el mismo efecto. Aquí hay un ejemplo de trabajo simple con aceleración 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;

}