我正在比较C的ODE解算器(Boost ODEINT)和MatLab的ODE45。每组代码都有相同的方程和初始条件。当在C代码中添加容错时，两个值(C++和MatLab)之间的差异相当大。两者的价值应该大致相同。增加误差容限会导致值的巨大差异是没有意义的。

#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;
 }

以下是MatLab代码：

%% 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