三体运动matlab代码

三体运动是一个经典的天体力学问题，可以使用 Matlab 来进行数值计算模拟。以下是一个简单的三体运动的 Matlab 代码示例：

% 三体运动模拟
% 假设三个质点的质量相等

% 初始位置和速度 (单位：AU, yr^-1)
x1 = [-0.97000436; 0.24308753; 0];
v1 = [-0.46620368; -0.43236573; 0];

x2 = [0.97000436; -0.24308753; 0];
v2 = [-0.46620368; -0.43236573; 0];

x3 = [0; 0; 0];
v3 = [0.93240737; 0.86473146; 0];

% 时间步长和总时间
dt = 0.01;
T = 20;

% 计算每个质点的加速度
function a = acceleration(x1, x2, x3)
    G = 4*pi^2;
    r1 = norm(x1 - x2);
    r2 = norm(x2 - x3);
    r3 = norm(x3 - x1);
    a1 = G * (x2 - x1) / r1^3 + G * (x3 - x1) / r3^3;
    a2 = G * (x3 - x2) / r2^3 + G * (x1 - x2) / r1^3;
    a3 = G * (x1 - x3) / r3^3 + G * (x2 - x3) / r2^3;
    a = [a1; a2; a3];
end

% 初始化位置和速度数组
x = [x1, x2, x3];
v = [v1, v2, v3];

% 计算加速度数组
a = acceleration(x(:,1), x(:,2), x(:,3));

% 循环计算每个时间步长的位置和速度
for i = 1:T/dt
    % 计算位置和速度
    x = x + v * dt + a * dt^2 / 2;
    v = v + a * dt;
    
    % 计算加速度
    a = acceleration(x(:,1), x(:,2), x(:,3));
    
    % 绘制当前位置
    plot3(x(1,:), x(2,:), x(3,:), 'o');
    xlim([-1.5 1.5]);
    ylim([-1.5 1.5]);
    zlim([-1.5 1.5]);
    drawnow;
end

这段代码会将三个质点的初始位置和速度输入进去，然后使用数值计算的方法来模拟它们的运动轨迹。在循环中，每个时间步长计算位置、速度和加速度，然后绘制当前位置。最后，程序会输出一个三维图形，展示三个质点的运动轨迹。

需要注意的是，这个代码只是一个简单的示例，还有很多可以优化的地方，比如使用更高阶的数值计算方法，增加精度等等。