Octave求解非线性ODE系统绘图与论文结果不符的技术求助

我在Octave环境中使用lsode求解一组非线性常微分方程系统，但运行代码后生成的绘图与目标论文第66页的结果不一致。我的代码如下：

function main
    t=linspace(0, 4000, 8000);
    x0 = [1;5*3.14/360;5*3.14/360;0.175;0;0];
    x = lsode (@f, x0, t);
    subplot(2,2,1);plot(t,x(:,1));title("l");legend("fgo")
    subplot(2,2,2);plot(t,x(:,2));title("theta")
    subplot(2,2,3);plot(t,x(:,3));title("phi")
    subplot(2,2,4);plot(t,x(:,4));title("dl")
endfunction;
function xdot = f(x,t);
    T=0.000644;
    O=0.001049;
    xdot=zeros(6,1);
    xdot(1)=x(4);
    xdot(2)=x(5);
    xdot(3)=x(6);
    xdot(4)=T+x(1)*(x(5)+O)^2*cos(x(3))^2+x(1)*x(6)^2-O^2*x(1)*(1-3*cos(x(2))^2*cos(x(3))^2);
    xdot(5)=-2*x(4)/x(1)*(x(5)+O)+2*x(6)*(x(5)+O)*tan(x(3))-3*O^2*cos(x(2))*sin(x(2));
    xdot(6)=-2*x(4)/x(1)*x(6)-(x(5)+O)^2*cos(x(3))*sin(x(3))-3*O^2*cos(x(2))^2*cos(x(3))*sin(x(3));
endfunction;

以下是针对性的排查方向：

• ODE方程实现核对
  逐行对比论文第66页的原始方程，重点检查：

  • xdot(4)中各项的符号、系数（比如-O^2*x(1)*(1-3*cos(x(2))^2*cos(x(3))^2)是否和论文的引力/离心项一致）
  • xdot(5)中的2*x(6)*(x(5)+O)*tan(x(3))项，确认系数与符号是否匹配论文推导
  • xdot(6)中的耦合项-(x(5)+O)^2*cos(x(3))*sin(x(3))，检查三角函数组合是否正确
  • 确认论文中所有角度运算均为弧度制（Octave三角函数默认弧度，与你的初始条件转换一致）

• 参数与初始条件验证

  • 确认参数T、O的数值、物理意义、单位完全匹配论文定义
  • 核对初始条件x0的每个分量（l, theta, phi, dl, dtheta, dphi）是否与论文第66页的设定完全一致，比如dl初始值0.175是否符合要求

• 求解器精度调整
  若系统存在刚性，lsode默认精度可能不足，尝试提高求解精度：

  options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
  x = lsode (@f, x0, t, options);
  

  同时检查时间区间linspace(0, 4000, 8000)是否与论文的仿真时长、采样密度一致

• 可视化匹配修正

  • 论文绘图可能将弧度转换为角度显示，可修改绘图代码验证：

    subplot(2,2,2);plot(t,x(:,2)*360/(2*pi));title("theta (deg)")
    subplot(2,2,3);plot(t,x(:,3)*360/(2*pi));title("phi (deg)")
    

  • 检查坐标轴范围是否与论文一致，避免因绘图缩放导致趋势误判

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