用matlab写一个模拟三体运行的代码

以下是使用MATLAB实现三体问题模拟的示例代码，其中使用了RK4数值积分方法：

% 定义常数
G = 6.67408e-11; % 万有引力常数
m1 = 5.97e24; % 地球质量
m2 = 7.34e22; % 月球质量
m3 = 1.99e30; % 太阳质量

% 定义初始状态
r1 = [-1.4710e11, 0]; % 地球初始位置，单位：m
v1 = [0, -3.0287e4]; % 地球初始速度，单位：m/s
r2 = [3.8440e8, 0]; % 月球初始位置，单位：m
v2 = [0, 1023]; % 月球初始速度，单位：m/s
r3 = [0, 0]; % 太阳初始位置，单位：m
v3 = [0, 0]; % 太阳初始速度，单位：m/s

% 定义积分步长和时长
dt = 60; % 秒
tmax = 365*24*3600; % 秒

% 定义初始状态向量
y0 = [r1, v1, r2, v2, r3, v3];

% 定义ODE函数
f = @(t, y) [
    y(4:5), 
    -G*m3*(y(1:2)-y(11:12))/norm(y(1:2)-y(11:12))^3-G*m2*(y(1:2)-y(7:8))/norm(y(1:2)-y(7:8))^3,
    y(8:9),
    -G*m3*(y(3:4)-y(11:12))/norm(y(3:4)-y(11:12))^3-G*m1*(y(3:4)-y(1:2))/norm(y(3:4)-y(1:2))^3,
    y(12:13),
    -G*m2*(y(5:6)-y(7:8))/norm(y(5:6)-y(7:8))^3-G*m1*(y(5:6)-y(1:2))/norm(y(5:6)-y(1:2))^3
];

% 使用RK4数值积分方法
[t, y] = rk4(f, y0, 0, tmax, dt);

% 绘制轨迹
figure;
hold on;
plot(y(:, 1), y(:, 2), 'b');
plot(y(:, 7), y(:, 8), 'r');
plot(y(:, 13), y(:, 14), 'y');
axis equal;
xlabel('X (m)');
ylabel('Y (m)');
legend('Earth', 'Moon', 'Sun');

% RK4数值积分方法的实现
function [t, y] = rk4(f, y0, t0, tmax, dt)
    t = t0:dt:tmax;
    y = zeros(length(t), length(y0));
    y(1, :) = y0;
    for i = 2:length(t)
        k1 = dt * f(t(i-1), y(i-1, :));
        k2 = dt * f(t(i-1) + dt/2, y(i-1, :) + k1/2);
        k3 = dt * f(t(i-1) + dt/2, y(i-1, :) + k2/2);
        k4 = dt * f(t(i-1) + dt, y(i-1, :) + k3);
        y(i, :) = y(i-1, :) + (k1 + 2*k2 + 2*k3 + k4)/6;
    end
end

上述代码将三个天体的初始位置和速度作为输入，并使用RK4数值积分方法进行模拟。最后绘制了三个天体的轨迹。需要注意的是，由于三体问题的复杂性，模拟结果可能存在误差。