MATLAB三维二阶ODE引力方程积分轨道绘制异常排查求助

问题背景 #

我需要在MATLAB中积分经典引力方程，求解物体的三维位置（本质是求解三维二阶ODE）。预期得到弹簧状的周期性重复轨道，但实际运行代码后得到的是类似抛射体的开轨迹，自查代码未发现明显错误，希望有人帮忙排查问题。

相关图示说明 #

• 引力方程：经典平方反比形式的引力加速度公式
• 预期正确轨道：呈弹簧状的闭合/周期性三维轨道
• 当前错误轨道：类似抛射体的非闭合轨迹

完整代码 #

clc; clear;

%% 参数设置
G = 6.67430e-11;  % J*m/kg^2 引力常数
M = 5.9722e24; % kg, 地球质量

%% 初始条件
% x0 y0 z0 xdot0 ydot0 zdot0
initial_state = [6.0827e6; 3.4907e6; 4.6517e6; -651.3041; 1082.1; 7426.4];  % 示例初始状态

%% 时间区间
n = 0.001078; % rad/s
p = 2*pi/n; % 理论周期
tspan = [0 p];  % 仿真时间范围

%% 求解ODE
[t, result] = ode45(@(t, y) unconstranined_ECI_ode(t, y, G, M), tspan, initial_state);

% 提取ECI坐标系位置分量
X = result(:, 1);
Y = result(:, 2);
Z = result(:, 3);
ECImatrix = [X Y Z];

%% 绘图
figure;

subplot(1, 2, 1);
plot(X, Z);
title('无约束运动xz平面投影');
xlabel('X');
ylabel('Z');
grid on;

subplot(1, 2, 2)
plot(Y, Z);
title('无约束运动yz平面投影');
xlabel('Y');
ylabel('Z');
grid on;

%% ODE定义函数
function state = unconstranined_ECI_ode(t, y, G, M)

    vel = y(4:6); % 速度分量
    r = y(1:3); % 位置分量
    rnorm = norm(r);
    rhat = r/rnorm;
    accel = -(G*M/rnorm^2)*rhat;
    state = [vel;accel];
    
end

补充说明 #

引力方程核心为平方反比形式：$\boldsymbol{a} = -\frac{GM}{r^2}\hat{\boldsymbol{r}}$，初始状态采用给定的位置与速度值，时间区间设置为一个理论周期，但结果仍未呈现周期性轨道。

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