matlab太阳系内行星运动轨迹

在MATLAB中，可以使用数值计算工具箱中的ode45函数来求解太阳系内行星运动轨迹问题。以下是一个示例代码，可以绘制太阳系内行星的运动轨迹：

G = 6.674e-11; % 万有引力常数
M = [1.989e30, 3.301e23, 4.869e24, 5.972e24, 6.419e23, 1.899e27, 5.685e26, 8.683e25]; % 太阳、水星、金星、地球、火星、木星、土星、天王星、海王星的质量
r0 = [147e9, 0, 0, 0, -35.98e6, 0, 778.57e9, 0, 0, -227.94e9, 0, 0, 0, 1427.0e9, 0, 0, 0, 4497.1e9, 0, 0]; % 初始位置和速度（假设所有行星在平面上运动，单位是米和米/秒）

% 定义ODE
f = @(t, r) nbody_ode(t, r, M, G);

% 求解ODE
[t, y] = ode45(f, [0, 165*365*24*3600], r0);

% 绘制轨迹
figure;
hold on;
colors = ['r', 'b', 'g', 'm', 'c', 'y', 'k', 'r', 'b'];
for i = 1:9
    plot(y(:, 4*i-3)/1e9, y(:, 4*i-2)/1e9, colors(i));
end
axis equal;
xlabel('x (10^9 m)');
ylabel('y (10^9 m)');
title('Solar System');
legend('Sun', 'Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune');

其中，nbody_ode是一个自定义的函数，可以计算所有天体的运动状态。如果你需要使用该代码，可以在同一文件夹下创建一个名为nbody_ode.m的文件，并将以下代码复制到该文件中：

function drdt = nbody_ode(t, r, M, G)
% 计算太阳系内所有天体的运动状态
n = length(M); % 天体数目
drdt = zeros(4*n, 1); % 初始化运动状态

% 计算速度和加速度
v = r(n+1:end);
a = zeros(3*n, 1);
for i = 1:n
    for j = 1:n
        if i ~= j
            rij = r(4*(j-1)+1:4*(j-1)+3) - r(4*(i-1)+1:4*(i-1)+3);
            a(3*(i-1)+1:3*i) = a(3*(i-1)+1:3*i) + G*M(j)*rij/norm(rij)^3;
        end
    end
end

% 计算运动状态
drdt(1:n) = v;
drdt(n+1:end) = a;
end

运行该代码，即可得到太阳系内行星的运动轨迹。你可以根据需要调整初始位置和速度、ODE以及绘图参数。