C++ BOOST ODEINT与MATLAB ODE45对比时误差容限异常问题

我正在对比C的Boost ODEINT常微分方程（ODE）求解器与MATLAB的ODE45求解器，二者使用完全相同的微分方程与初始条件。但为C代码添加误差容限后，两者计算结果差异极大——理论上二者结果应大致一致，误差容限的添加不应导致如此显著的差异。

C++代码 #

#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

问题排查关键点 #

• 浮点数精度不匹配：C代码里mu、r_2、r用了单精度float，而MATLAB默认是双精度double。单精度有效位数仅7位左右，远低于双精度的15-16位，积分过程中会积累大量截断误差，长期运行后结果必然偏差严重。把C中所有float替换为double，包括mu的定义和r_2、r的变量类型。
• 输出格式干扰：C++的write_lorenz函数同时用了'\n'和endl，endl会强制刷新缓冲区并额外输出换行，导致输出文件出现空行，MATLAB读取时可能解析错误。删掉<< endl，只保留<< '\n'。
• 误差容限不一致：当前C用的是1E-9的严格容限，MATLAB用的是1e-6。虽然更严格的容限理论上更准确，但如果因为单精度问题导致求解器步长调整异常，反而会引入问题。先把C的容限改成和MATLAB一致的1e-6测试。
• 积分函数调用逻辑：integrate_const是固定步长输出，结合make_dense_output的自适应步长求解器时，内部插值输出逻辑和MATLAB ODE45的插值可能有细微差异，但这不是核心问题，先解决精度问题再验证。

内容的提问来源于stack exchange，提问作者Jon Snow