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.

Answer 2

The Matlab ode45 function already includes error control and I think also interpolation (dense output). to compare with boost.odeint you should use the same functionality there. Boost.odeint provides integrate functions that perform step-size control and dense output if the used stepper algorithm provides this functionality. The following code piece shows how you this is used with the default error control parameters from Matlab given by horchler:

#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 damped_osc( const state_type &x , state_type &dx , const double t )
{
    dx[0] =  x[1];
    dx[1] = -x[0] - gam*x[1];
}

void print( const state_type &x, const double t )
{
    cout << x[0] << endl;
}

int main(int argc, char **argv)
{
    cout.precision(16);  // full precision output
    const double dt = 0.1;
    typedef runge_kutta_dopri5<state_type> stepper_type;
    state_type x(2);
    x[0] = 1.0;
    x[1] = 0.0;

    integrate_const(make_dense_output<stepper_type>( 1E-6, 1E-3 ),
                    damped_osc, x, 0.0, 10.0, dt , print);

    return 0;
}

Please note that the results might still not be exactly the same (as in all 16 digits) because the error control in Boost.odeint might not be impemented exactly as in Matlab's ode45.