MATLAB中ode45求解四体问题时积分错误的解决方法

1. 修正万有引力常数G的取值 #

你代码中G的取值是6.6742*10^11，但实际万有引力常数的正确值为6.6742e-11 N·m²/kg²。错误的G值会导致物体的引力加速度达到天文数字（量级为1e21 m/s²），ODE求解器为了追踪这种极端变化会不断缩小步长，最终触发最小步长限制报错。这是问题的核心原因。

2. 优化微分方程的索引逻辑 #

原代码的索引计算虽然结果正确，但可读性差且容易出错。可以重构为更直观的索引方式，明确每个物体的位置、速度、加速度对应区间：

function dubdt = eq(~, ub, masses, G)
    dubdt = zeros(24,1);
    % 位置的导数等于速度
    dubdt(1:3) = ub(4:6);
    dubdt(7:9) = ub(10:12);
    dubdt(13:15) = ub(16:18);
    dubdt(19:21) = ub(22:24);
    
    % 计算每个物体的加速度
    for i = 1:4
        pos_i = ub((6*(i-1)+1):(6*(i-1)+3));
        acc = zeros(3,1);
        for j = 1:4
            if i ~= j 
                pos_j = ub((6*(j-1)+1):(6*(j-1)+3));
                r_ij = pos_i - pos_j;
                dist = norm(r_ij);
                acc = acc - G * masses(j) * r_ij / (dist^3);
            end
         end
        dubdt((6*(i-1)+4):(6*(i-1)+6)) = acc;
    end
end

3. 更换刚性ODE求解器 #

你的四体系统中恒星与行星质量差异达24个量级，属于刚性系统。ode45是非刚性求解器，处理这类系统效率极低且容易发散，建议更换为专门的刚性求解器ode15s，同时可调整精度参数提升稳定性：

options = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t, result] = ode15s(@(t,ub) eq(t, ub, masses, G), [0,time], initials, options);

完整修正后的代码 #

clear;clc;close all;
ms1 = 10^30;ms2 = 10^30;mp1 = 10^6;mp2 = 10^6;
masses = [ms1;ms2;mp1;mp2];
G = 6.6742e-11; % 修正万有引力常数
time = 30*86400;
%[r,v]
us1_0 = [10^10; 0; 0; 0; 40000; 0];
us2_0 = [-10^10; 0; 0; 0; -40000; 0];
up1_0 = [1.5*10^10; 0; 0; 0; 160000; 0];
up2_0 = [-1.5*10^10; 0; 0; 0; -160000; 0];
initials = [us1_0; us2_0; up1_0; up2_0];

% 使用刚性求解器并设置精度
options = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t, result] = ode15s(@(t,ub) eq(t, ub, masses, G), [0,time], initials, options);

figure
x1 = result(:,1);y1 = result(:,2);z1 = result(:,3);
x2 = result(:,7);y2 = result(:,8);z2 = result(:,9);
x3 = result(:,13);y3 = result(:,14);z3 = result(:,15);
x4 = result(:,19);y4 = result(:,20);z4 = result(:,21);

grid on;
plot3(x1,y1,z1);hold on
plot3(x2,y2,z2);hold on
plot3(x3,y3,z3);hold on
plot3(x4,y4,z4);
legend('star1','star2','planet1','planet2')
xlabel('X (m)');ylabel('Y (m)');zlabel('Z (m)');
title('四体问题轨迹模拟')


function dubdt = eq(~, ub, masses, G)
    dubdt = zeros(24,1);
    % 位置的导数等于速度
    dubdt(1:3) = ub(4:6);
    dubdt(7:9) = ub(10:12);
    dubdt(13:15) = ub(16:18);
    dubdt(19:21) = ub(22:24);
    
    % 计算每个物体的加速度
    for i = 1:4
        pos_i = ub((6*(i-1)+1):(6*(i-1)+3));
        acc = zeros(3,1);
        for j = 1:4
            if i ~= j 
                pos_j = ub((6*(j-1)+1):(6*(j-1)+3));
                r_ij = pos_i - pos_j;
                dist = norm(r_ij);
                acc = acc - G * masses(j) * r_ij / (dist^3);
            end
         end
        dubdt((6*(i-1)+4):(6*(i-1)+6)) = acc;
    end
end

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