Comparison of odeint's runge_kutta4 with Matlab's ode45

Question

I would like to use runge_kutta4 method in the odeint C++ library. I've solved the problem in Matlab. My following code in Matlab to solve x'' = -x - g*x', with initial values x1 = 1, x2 = 0, is as follows

main.m

clear all
clc
t = 0:0.1:10;
x0 = [1; 0];
[t, x] = ode45('ODESolver', t, x0);
plot(t, x(:,1));
title('Position');
xlabel('time (sec)');
ylabel('x(t)');

ODESolver.m

function dx = ODESolver(t, x)
dx = zeros(2,1);
g  = 0.15;
dx(1) =  x(2);
dx(2) = -x(1) - g*x(2);
end

I've installed the odeint Library. My code for using runge_kutta4 is as follows

#include <iostream>

#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;

/* The type of container used to hold the state vector */
typedef std::vector< double > state_type;

const double gam = 0.15;

/* The rhs of x' = f(x) */
void lorenz( const state_type &x , state_type &dx ,  double t )
{
    dx[0] =  x[1];
    dx[1] = -x[0] - gam*x[1];
}


int main(int argc, char **argv)
{
    const double dt = 0.1;
    runge_kutta_dopri5<state_type> stepper;
    state_type x(2);
    x[0] = 1.0;
    x[1] = 0.0;


   double t = 0.0;
   cout << x[0] << endl;
   for ( size_t i(0); i <= 100; ++i){
       stepper.do_step(lorenz, x , t, dt );
       t += dt;
       cout << x[0] << endl;
   }


    return 0;
}

The result is in the following picture

My question is why the result varies? Is there something wrong with my C++ code?

These are the first values of both methods

Matlab    C++
-----------------
1.0000    0.9950
0.9950    0.9803
0.9803    0.9560
0.9560    0.9226
0.9226    0.8806
0.8806    0.8304
0.8304    0.7728
0.7728    0.7084
0.7083    0.6379

---

Update:

The problem is that I forgot to include the initial value in my C++ code. Thanks for @horchler for noticing it. After including the proper values and using runge_kutta_dopri5 as @horchler suggested, the result is

Matlab    C++
-----------------
1.0000    1.0000
0.9950    0.9950
0.9803    0.9803
0.9560    0.9560
0.9226    0.9226
0.8806    0.8806
0.8304    0.8304
0.7728    0.7728
0.7083    0.7084

I've updated the code to reflect these modifications.

Answer 1

The runge_kutta4 stepper in odeint is nothing like Matlab's ode45, which is an adaptive scheme based on the Dormand-Prince method. To replicate Matlab's results, you should probably try the runge_kutta_dopri5 stepper. Also, make sure that your C++ code uses the same absolute and relative tolerances as ode45 (defaults are 1e-6 and 1e-3, respectively). Lastly, it looks like you may not be saving/printing your initial condition in your C++ results.

If you're confused at why ode45 is not taking fixed steps even though you specified t = 0:0.1:10;, see my detailed answer here.

If you must use the fixed steprunge_kutta4 stepper, then you'll need to reduce the integration step-size in your C++ code to match Matlab's output.