近地轨道航天器发动机运动ODE仿真仅显示点的问题求助

我用Matlab实现了近地轨道航天器的运动方程求解，但绘制ODE解的结果时只显示一个点，无法定位问题。相关代码如下：

function sol = orbital_motion()
% Define the function that returns the derivative of the state vector
    function dydt = ode(t, y)
        % Extract the state variables
        X = y(1);
        Y = y(2);
        Z = y(3);
        Xdot = y(4);
        Ydot = y(5);
        Zdot = y(6);
     
        Xm = 278,481; Ym = 378,481; Zm = 268,481;
        Xs = 147098,074 ; Ys = 696,166 ; Zs = 1515,776 ;
        mu_e = 3.986135e14; %sets the gravitational parameter for the Earth.
        mu_m = 4.89820e12;%sets the gravitational parameter for the Moon.
        mu_s = 1.3253e20;%sets the gravitational parameter for the Sun.
        a = 6.37826e6; %the equatorial radius of the Earth.
        J = 1.6246e-3;%sets the Earth's second dynamic form factor.
        % Calculate the distances and other parameters
        r = sqrt(X^2 + Y^2 + Z^2);
        rm = sqrt(Xm^2 + Ym^2 + Zm^2);
        delta_m = sqrt((X - Xm)^2 + (Y - Ym)^2 + (Z - Zm)^2);
        delta_s = sqrt((X - Xs)^2 + (Y - Ys)^2 + (Z - Zs)^2);
        Y_s = Ys;
        Z_s = Zs;
        % Calculate the derivatives
        X2dot = -((mu_e*Y)/r^3) * (1 + (J*(a/r)^2)*(1 - 5*Z^2/r^2)) ...
            - (mu_m*(X - Xm)/delta_m^3) ...
            - (mu_m*Xm/rm^3) ...
            - (mu_s*(X - Xs)/delta_s^3) ...
            - (mu_s*Xs/r^3);
        Y2dot = -((mu_e*Y)/r^3) * (1 + (J*(a/r)^2)*(1 - 5*Z^2/r^2)) ...
            - (mu_m*(Y - Ym)/delta_m^3) ...
            - (mu_m*Ym/rm^3) ...
            - (mu_s*(Y - Y_s)/delta_s^3) ...
            - (mu_s*Ys/r^3);
        Z2dot = -((mu_e*Z)/r^3) * (1 + (J*(a/r)^2)*(3 - 5*Z^2/r^2)) ...
            - (mu_m*(Z - Zm)/delta_m^3) ...
            - (mu_m*Zm/rm^3) ...
            - (mu_s*(Z - Z_s)/delta_s^3) ...
            - (mu_s*Zs/r^3);
     
        % Return the derivative vector
        dydt = [Xdot; Ydot; Zdot; X2dot; Y2dot; Z2dot];
    end

% Set the initial conditions
X0 = 0;
Y0 = 0;
Z0 = 0;
Xdot0 = 11436.6;
Ydot0 = 0;
Zdot0 = 0;
y0 = [X0; Y0; Z0; Xdot0; Ydot0; Zdot0];

% Define the time span
tf = 20;
tspan = [0, tf];

% Call ODE45
%options = odeset('RelTol', 1e-8, 'AbsTol
# [t,y] = ode45(@ode,[0 20],y0);
plot(t,y(:,1),'-o',t,y(:,2),'-o',t,y(:,3),'-o')
title('Solution of the Motion Equation with gODE45');
end

关键错误分析 #

• ODE求解代码未执行：代码中# [t,y] = ode45(@ode,[0 20],y0);的#不是Matlab注释符号，会触发语法错误，导致求解逻辑未运行；同时前面的%options = odeset('RelTol', 1e-8, 'AbsTol语句未闭合，也会引发语法中断。
• 初始位置致命错误：初始条件设为X0=Y0=Z0=0，即航天器位于地球质心，此时r=sqrt(X^2+Y^2+Z^2)=0，后续计算r^3会出现除以0的错误，导致ode45仅能计算初始点后就终止。
• 数值格式错误：变量赋值中的逗号（如Xm = 278,481）在Matlab中是语句分隔符，会被解析为Xm=278，多余的481会引发未定义变量错误；小数点需用英文点号，原147098,074应改为147098.074。
• 运动方程逻辑错误：X2dot的第一项错误使用Y而非X，导致X方向引力加速度计算完全错误，正确应为-((mu_e*X)/r^3) * ...。

修复建议 #

1. 修正ODE调用语句：删除#以执行ode45，同时补全未闭合的odeset配置（如options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);）。
2. 调整初始位置：设置合理的近地轨道初始坐标，比如X0 = 6.37826e6 + 400e3（地球赤道半径+400km轨道高度），Y0=Z0=0。
3. 修正数值格式：去掉所有千位分隔的逗号，将小数点统一改为英文点号。
4. 修复运动方程：将X2dot中的mu_e*Y替换为mu_e*X，确保各坐标轴的引力加速度分量对应正确。

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