用matlab 的ode45模拟太阳系的轨迹

好的，模拟太阳系轨迹是一个非常有趣的问题。由于太阳系中天体的数量很多，模拟太阳系轨迹需要考虑这些天体之间的相互作用。下面是一个简单的示例，展示如何使用 Matlab 中的 ODE45 求解太阳系中行星的轨迹问题。

首先，我们需要设置太阳、行星和初始条件。在这个示例中，我们将考虑太阳、地球和火星之间的相互作用。我们假设太阳是一个固定的点，不考虑它对行星的引力影响。我们还需要定义每个天体的质量和初始位置。

% 设置太阳、地球和火星的初始位置、质量等参数
G = 6.67428e-11; % 重力常数
Ms = 1.9891e30; % 太阳质量
Me = 5.9742e24; % 地球质量
Mm = 6.4191e23; % 火星质量

% 设置太阳、地球和火星的初始位置（x, y, z）和速度（vx, vy, vz）
Rs = [0, 0, 0];
Vs = [0, 0, 0];
Re = [1.496e11, 0, 0];
Ve = [0, 2.9783e4, 0];
Rm = [2.2794e11, 0, 0];
Vm = [0, 2.4077e4, 0];

% 连接所有天体的初始位置和速度
y0 = [Rs, Vs, Re, Ve, Rm, Vm];

% 计算天体之间的距离
dse = norm(Re-Rs);
dsm = norm(Rm-Rs);
dem = norm(Rm-Re);

% 计算地球和火星的初始速度
Ve = Ve + [0, 0, G*Ms/dse^2];
Vm = Vm + [0, 0, G*Ms/dsm^2];

接下来，我们需要定义一个函数，用于计算天体之间的相互作用。在这个函数中，我们需要计算每个天体受到其他天体引力的加速度。

% 定义函数，计算天体之间的相互作用
function dydt = planet_orbit(t, y)
    % 设置太阳、地球和火星的质量
    G = 6.67428e-11;
    Ms = 1.9891e30;
    Me = 5.9742e24;
    Mm = 6.4191e23;
    
    % 从 y 中提取每个天体的位置和速度
    Rs = y(1:3);
    Vs = y(4:6);
    Re = y(7:9);
    Ve = y(10:12);
    Rm = y(13:15);
    Vm = y(16:18);
    
    % 计算每个天体受到的引力和加速度
    Fse = G*Me*Ms/norm(Re-Rs)^2;
    Fes = -Fse;
    Fsm = G*Mm*Ms/norm(Rm-Rs)^2;
    Fms = -Fsm;
    Fem = G*Mm*Me/norm(Rm-Re)^2;
    Fme = -Fem;
    
    ase = Fse/Me;
    aes = Fes/Ms;
    asm = Fsm/Mm;
    ams = Fms/Ms;
    aem = Fem/Me;
    ame = Fme/Mm;
    
    % 设置返回值
    dydt = [Vs; ase; aes; Ve; aem; ame; Vm; asm; ams];
end

最后，我们使用 ODE45 函数来求解行星的轨迹。在这个示例中，我们将计算 10 年的时间，每隔 1 小时计算一次行星的位置和速度。

% 求解行星的轨迹
tspan = [0, 10*365*24*3600];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[T, Y] = ode45(@planet_orbit, tspan, y0, options);

% 绘制行星的轨迹
figure;
plot3(Y(:, 1), Y(:, 2), Y(:, 3), 'b');
hold on;
plot3(Y(:, 7), Y(:, 8), Y(:, 9), 'g');
plot3(Y(:, 13), Y(:, 14), Y(:, 15), 'r');
axis equal;
xlabel('X');
ylabel('Y');
zlabel('Z');
title('太阳系行星轨迹');
legend('地球', '火星', '太阳');

执行上述代码后，会得到一个三维图形，显示地球、火星和太阳之间的相对位置关系。