11
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Hi @Marcin, I edited your title to better reflect what I think your problem actually is. I think we may get more useful answers and it will be more searchable for others that see this question in the future with the new title. Feel free to change it back if you disagree. $\endgroup$ – Doug Lipinski Mar 1 '15 at 3:35
16
$\begingroup$

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;

}
$\endgroup$
  • 6
    $\begingroup$ "Unfortunately C++ does not natively support vector operations like this" I think it does, in the standard library even, but not necessarily easy to use with other linear algebra libraries: en.cppreference.com/w/cpp/numeric/valarray I think common linear algebra libraries like Eigen, also should count as "support". $\endgroup$ – Kirill Mar 1 '15 at 8:18
  • 1
    $\begingroup$ @Kirill, Thanks for the tip. I'm still relatively new to C++ and I've not used valarray before, I just learned something useful too! Editing to add. $\endgroup$ – Doug Lipinski Mar 1 '15 at 13:10
  • 1
    $\begingroup$ Maybe this advice will also be helpful then: 1) Use clang-format for automatically formatting your code, it's really standard and uniform. 2) Use typedef std::valarray<double> Vector for commonly used types. 3) Use const int NDIM = 2 instead of #define for type safety and correctness. 4) Since C++11 you can replace the body of RHS simply with return {X[1], 1} . 5) It's really uncommon in C++ (unlike C) to first declare variables, then later initialize them, prefer declaring variables at the same place where you initialize them (double t = 0., etc.) $\endgroup$ – Kirill Mar 1 '15 at 18:09
  • 1
    $\begingroup$ @MarcinW. RHS() computes the Right Hand Side of the differential equation. The state vector X is (x, v) so dX/dt = ( dx/dt, dv/dt ) = ( v, a ). For your problem (if a=G*M/x^2) RHS should return { X[1], G*M/(X[0]*X[0]) }. $\endgroup$ – Doug Lipinski Mar 2 '15 at 0:04
  • 1
    $\begingroup$ @Kirill I know, but that only works since C++11 which means it does not work with default compiler options on the most popular compilers. I chose to leave that out in favor of something that works with the old standards as well and hopefully reduce confusion caused by inability to compile the code. $\endgroup$ – Doug Lipinski Mar 2 '15 at 0:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.