Error tolerance not working in C++, when plotted against ODE45 in MATLAB

Question

I am comparing the ODE solver of C++ (BOOST ODEINT) against MATLAB's ODE45. Each set of code has the same equations and initial conditions. When error tolerance is added to the C++ code, the difference between the two values (C++ and MATLAB) is quite large. The values of the two should be roughly the same. It does not make sense that the error tolerance addition would cause a large difference in values.

    #define _CRT_SECURE_NO_WARNINGS
    #include <iostream>
    #include <fstream>
    #include <sstream>
    #include <boost/numeric/odeint.hpp>
    #include <cmath>
    #include <math.h>
    #include <boost/array.hpp>
    #include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>


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

    const double sigma = 10.0;
    const double R = 28.0;
    const double b = 8.0 / 3.0;
    const float mu = 3.986e5;

    typedef boost::array< double, 6 > state_type;

    ofstream opfile;

    void lorenz(const state_type& x, state_type& dxdt, double t)
   {

        float r_2 = (pow(x[0],2) + pow(x[1],2) + pow(x[2],2));
        float r = sqrt(r_2);

        dxdt[0] = x[3];
        dxdt[1] = x[4];
        dxdt[2] = x[5];
        dxdt[3] = (-mu / pow(r, 3)) * x[0];
        dxdt[4] = (-mu / pow(r, 3)) * x[1];
        dxdt[5] = (-mu / pow(r, 3)) * x[2];

    }

     void write_lorenz(const state_type& x, const double t)
    {

      opfile << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2]<< '\t' << x[3] << '\t' 
      << x[4] << '\t' << x[5] << '\n'<< endl;

     }

   int main(int argc, char** argv)
    {

       double T = 4.2706e4;
       typedef runge_kutta_dopri5<state_type> stepper_type;
       state_type x = {3.1117e3 , 0 , 4.565e3, 3.0, 4.5, 6.8}; // initial conditions
       const double dt = 0.1;
       cout.precision(16);  // full precision output
    
       opfile.open("test_rk_dopri.txt");



      integrate_const(make_dense_output<stepper_type>(1E-9, 1E-9), lorenz, x, 0.0, T, dt, write_lorenz);

//integrate_const(make_dense_output(1.0e-6, 1.0e-6, runge_kutta_dopri5< state_type >()), lorenz, x, 0.0, T, dt, write_lorenz);

     opfile.close();

return 0;
 }

Here is the MATLAB code:

      %% ODE45 plotted against C++ odeint solvers
      %% initial conditions
      r0_vec = [3.1117e3;0;4.565e3];
      v0_vec = [3.0; 4.5; 6.8];
      mew = 3.986e5;


     %% Runge Kutta cash karp
     fid = fopen('test_cashkarp_Rk.txt', 'rt');
     a = textscan(fid, '%f%f%f%f%f%f%f');
     fclose(fid);
     rk_cashkarp = cell2mat(a);
     % t = linspace(0.1,10,length(rk_cashkarp));

     fid_1 = fopen('test_rk4.txt', 'rt');
     b = textscan(fid_1, '%f%f%f%f%f%f%f');
     fclose(fid_1);
     rk_4 = cell2mat(b);
     % t = linspace(0.1,10,length(rk_4));

     fid_2 = fopen('test_rk_dopri.txt', 'rt');
      c = textscan(fid_2, '%f%f%f%f%f%f%f');
      fclose(fid_2);
     rk_dopri = cell2mat(c);
     % t = linspace(0.1,10,length(rk_dopri));

     fid_3 = fopen('test_rk_fehlberg.txt', 'rt');
     d = textscan(fid_3, '%f%f%f%f%f%f%f');
     fclose(fid_3);
     rk_fehlberg = cell2mat(d);
     % t = linspace(0.1,10,length(rk_fehlberg));
     T =  4.2706e4;


    %%ODE45 function
     options = odeset('relTol',1e-6,'absTol',1e-6);
     tspan = rk_dopri(:,1);
       [t,y] = ode45(@(t,y) ODE_eqnsofMot(t,y,mew),tspan, [r0_vec; v0_vec],options);
       close all

      figure ('name','Postion Ode45')
      plot(rk_dopri(:,2)-y(:,1))
      % plot3(y(:,1),y(:,2),y(:,3),'.')


      % 
      % figure ('name','Postion rk4')
      % plot3(rk_4(:,2),rk_4(:,3),rk_4(:,4),'.')
          % 
      % figure ('name','Postion rkcashkarp')
      % plot3(rk_cashkarp(:,2),rk_cashkarp(:,3),rk_cashkarp(:,4),'.')
      % 
      % figure ('name','Postion rkdopri')
      % plot3(rk_dopri(:,2),rk_dopri(:,3),rk_dopri(:,4),'.')
      % 
      % figure ('name','Postion  rkfehlberg')
      % plot3(rk_fehlberg(:,2),rk_fehlberg(:,3),rk_fehlberg(:,4),'.')
      % 


      %%ODE45 function

      function dydt = ODE_eqnsofMot(t, y, mew)

      dydt = zeros(6,1); 
      r = sqrt(y(1)^2+(y(2))^2+(y(3))^2);

      dydt(1) = y(4);
      dydt(2) = y(5);
      dydt(3) = y(6);
      dydt(4) = (-mew/(r)^3)*y(1);
      dydt(5) = (-mew/(r)^3)*y(2);
      dydt(6) = (-mew/(r)^3)*y(3);


  end

Answer 1

This is not a difference between the different implementations (I think they all follow the standard that Hairer and Shampine established). I get the same graph in python with the following short lines

def grav(t,u): p,v = u[:3],u[3:]; a = -mu*p*sum(p**2)**-1.5; return *v, *a

mu, p0, v0 = 3.986e5, [3.1117e3 , 0 , 4.565e3], [3.0, 4.5, 6.8]
u0 = [*p0, *v0]
t = np.arange(0,4.2706e4, 0.1)

u6 = solve_ivp(grav, t[[0,-1]], u0, t_eval=t, atol=1e-6, rtol=1e-6)
u9 = solve_ivp(grav, t[[0,-1]], u0, t_eval=t, atol=1e-9, rtol=1e-9)

plt.plot(t,u9.y[0]-u6.y[0]); plt.grid(); plt.show()

Of course relative to the range [-1000, 4000] of that component this is an error of 0.2 % at the extreme, so below a visual difference if the graphs of the components are plotted together.

Still, it is far above the desired error level. Let's change that by providing a scale for the absolute error

sc = np.concatenate([[1e3]*3,[1]*3])
def solve(tol): return solve_ivp(grav, t[[0,-1]], u0, t_eval=t, atol=sc*tol, rtol=tol)
u6 = solve(1e-6)
u7 = solve(1e-7)
u8 = solve(1e-8)
u9 = solve(1e-9)

That does not work as expected, between u6 and u9 the range of the error reduces by half, not by orders of magnitude. The error from u8 to u9 has a range to 0.2, from u7 to u9 to 1.5.

In a log-plot of the full position difference

plt.semilogy(t,norm(u6.y[:3]-u9.y[:3])/norm(u9.y[:3])); 
plt.grid(); plt.show()

one can see that the desired error level holds to about t=3000

After that it is more like the time axes are shifted, so that the position difference is more of being on the same orbit at different times than errors of the trajectory. Exactly the same pattern is obtained when comparing u9 to u12, the latter at tol=1e-12, only with the error level reduced by 1e-3.

This is likely an effect of long-term integration, either that the method step error has a noticeable component in direction of the trajectory, changing the speed along it, or that accumulated errors in the time step summation become noticeable.