Comparison of odeint runge_kutta4 with matlab ode45
I would like to use a method runge_kutta4
in the odeint C ++ library . I solved the problem in Matlab. My next code in Matlab for a solution x'' = -x - g*x'
with initial values x1 = 1
is x2 = 0
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 have installed the odeint library. My code to use is runge_kutta4
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 shown in the following figure
My question is, why does the result change? 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, I forgot to include the seed in my C ++ code. Thanks for noticing @horchler. After including the correct values and using runge_kutta_dopri5
as suggested by @horchler, 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 my code to reflect these changes.
source to share
The pedometer runge_kutta4
in odeint is not like Matlab ode45
, which is an adaptive circuitry based on Dormand-Prince . To reproduce the Matlab results you should probably try stepper runge_kutta_dopri5
. Also, make sure that your code is in C ++ using the same absolute and relative tolerances that ode45
(by default - 1e-6
and 1e-3
, respectively). Finally, it looks like you cannot store / print your initial condition in your C ++ results.
If you are confused as to why it ode45
doesn't take fixed steps even if you specified t = 0:0.1:10;
, see my detailed answer here .
If you have to use a stepping step runge_kutta4
, you will need to reduce the integration step in your C ++ code to match the Matlab output.
source to share
Matlab's ode45 feature already includes error management, and I think interpolation (tight inference) does too. for comparison with boost.odeint you should use the same functionality. Boost.odeint provides functions integrate
that perform stepping size control and dense output if the stepping algorithm used provides this functionality. The following piece of code shows you how you use this with the default error management options 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;
}
Note that the results may not be exactly the same (as with all 16 digits), because the error management in Boost.odeint may not be written exactly the same as in Matlab ode45.
source to share