| main | class hierarchy | alphabetical | list | functions | files | tutorial | programs |

4th order Runge-Kutta ODE solver

Reference	Tutorial	Example Code
smtk::rk4	4th order Runge-Kutta ODE solver	sine_rk4.cpp damped_SHO.cpp damped_pendulum.cpp lotka_volterra.cpp multi_plot.cpp

Example: damped simple harmonic oscillator

This example solves a damped simple harmonic oscillator. We have chosen this very simple system to solve to make this example simple.

First we include some header files:

#include <smtk.h>
#include <smtk/rk4.h>

The differential equations we wish to solve using the smtk::rk4 class. This smtk::rk4 template class is not declared in smtk.h so it does not gum-up the works when it's not being used. So the smtk::rk4 template class is declared in smtk/rk4.h. We wish to solve the second order ondinary differitial equation (ODE):

d²x/dt² = - x - b dx/dt

which we can we write as the following two first order ODEs

dx/dt = v

dv/dt = - x - b v .

We set up the two differential equations to be solved in the call-back function difeq().

static const float b = 0.02; // damping constant

static void difeq(float *xdot, const float *x, double t)
{
  xdot[0] = x[1];
  xdot[1] = - x[0] - b * x[1];
}

We declare main() and a smtk::rk4 object that uses an array of floats for state variables, a double for time, and a float for the time step:

int main(void)
{
  smtk::rk4<float,double,float>
    rk4(difeq, /* differential equations call-back function */
        2,     /* 2 equations */
        0.05f  /* time step */);

We declare the floats for the state variables, and a double for the time:

  float  x[2] = { 1.0f, 0.0f }; // 2 state variables  x, v
  double t    = 0.0;            // time

Now we go to time 100.0 printing between steps:

  while(t < 100.0)
    {
      printf("%.15g %.7g %.7g\n", t, x[0], x[1]);
      rk4.go(x, &t); // go one time step
    }

The printing will not slow down the program at every loop unless it's being flushed every loop, like if it writes to a terminal. If you pipe the output to a file it will not flush every loop. That's buffered I/O. It's that same if you use the C++ cout output class.

Print the last values calculated and exit with status success.

  printf("%.15g %.7g %.7g\n", t, x[0], x[1]);

  return 0;
}

You can quickly see the results of the simulation by piping the output of the program to quickplot by running

./sine_rk4 | quickplot

You can get quickplot at http://quickplot.sourceforge.net/.

There are many other way to use smtk::rk4(), see the smtk::rk4 reference manual for more details.

---