Applying the Runge-Kutta method to second order ODEs

Question

How can I replace the Euler method by Runge-Kutta 4th order to determine the free fall motion in not constant gravitional magnitude (eg. free fall from 10 000 km above ground)?

So far I wrote simple integration by Euler method:

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

x variable means current position, v means velocity, getMagnitude(x) returns acceleration on x position.

I tried implement RK4:

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

where rk4() function body is:

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

But something is wrong, because I am integrating only once using RK4 (acceleration). Integrating velocity using RK4 dosen't make sense because it's same as v * dt.

Could You tell me how to solve second order differential equations using Runge-Kutta integration? Should I implement RK4 by computing k1, l1, k2, l2 ... l4 coefficients? How can I do that?

Answer 1

There seems to be quite a bit of confusion about how to apply multi-step (e.g. Runge-Kutta) methods to 2nd or higher order ODEs or systems of ODEs. The process is very simple once you understand it, but perhaps not obvious without a good explanation. The following method is the one I find simplest.

In your case, the differential equation you would like to solve is $F = m\ddot{x}$. The first step is to write this second-order ODE as a system of first order ODEs. This is done as

$$ \left[\begin{array}{c} \dot{x} \\ \dot{v} \end{array}\right] = \left[\begin{array}{c} v \\ F/m \end{array}\right] $$

All equations in this system must be solved simultaneously, which is to say that you should not advance $v$ and then advance $x$, they should both be advanced at the same time. In languages that support vector operations without loops this is easily done by making all the necessary terms in your code vectors of length 2. The function that computes the right hand side (the rate of change) of your ODE should return a vector of length 2, k1 to k4 should be vectors of length 2, and your state variable $(x,v)$ should be a vector of length 2. In MATLAB the necessary code for the time stepping can be written as:

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

where X$=(x,v)$ and RHS( t, X ) returns a vector containing $(\dot{x}(t),\dot{v}(t))$. As you can see, by vectorizing things you don't even need to change the syntax of the RK4 code no matter how many equations are in your ODE system.

Unfortunately C++ does not natively support vector operations like this so you need to either use a vector library, use loops, or manually write out the separate parts. In C++ you can use std::valarray to achieve the same effect. Here's a simple working example with constant acceleration.

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

}