在matlab 上用欧拉法拟合日月地三体的运动轨迹的程序

以下是在MATLAB上使用欧拉法拟合日月地三体运动轨迹的示例程序：

% Parameters
G = 6.6743e-11; % Gravitational constant
m1 = 5.97e24; % Mass of Earth
m2 = 7.34e22; % Mass of Moon
m3 = 1.989e30; % Mass of Sun
r1 = 6378.1e3; % Radius of Earth
r2 = 1737.4e3; % Radius of Moon
r3 = 695700e3; % Radius of Sun
d = 384400e3; % Distance between Earth and Moon
t0 = 0; % Initial time
tf = 10*24*3600; % Final time (10 days in seconds)
h = 60; % Step size (60 seconds)

% Initial conditions
x1 = [r1+d; 0; 0]; % Position of Earth
v1 = [0; sqrt(G*m3/d); 0]; % Velocity of Earth
x2 = [d; 0; 0]; % Position of Moon
v2 = [0; sqrt(G*m3/d)+sqrt(G*m1/(d+r1)); 0]; % Velocity of Moon
y = [x1; v1; x2; v2]; % State vector

% Euler method integration
t = t0:h:tf;
n = length(t);
Y = zeros(12,n);
Y(:,1) = y;
for i = 2:n
    ydot = rhs(t(i-1),y,G,m1,m2,m3);
    y = y + h*ydot;
    Y(:,i) = y;
end

% Plotting
figure;
plot(Y(1,:),Y(2,:),'b','LineWidth',2);
hold on;
plot(Y(7,:),Y(8,:),'r','LineWidth',2);
axis equal;
xlabel('x (m)');
ylabel('y (m)');
title('Three-body problem: Earth-Moon-Sun');

% Right-hand side of the ODE system
function ydot = rhs(t,y,G,m1,m2,m3)
    r12 = norm(y(1:3)-y(7:9));
    r13 = norm(y(1:3)-y(13:15));
    r23 = norm(y(7:9)-y(13:15));
    a1 = -G*m2*(y(1:3)-y(7:9))/r12^3 - G*m3*(y(1:3)-y(13:15))/r13^3;
    a2 = -G*m1*(y(7:9)-y(1:3))/r12^3 - G*m3*(y(7:9)-y(13:15))/r23^3;
    a3 = -G*m1*(y(13:15)-y(1:3))/r13^3 - G*m2*(y(13:15)-y(7:9))/r23^3;
    ydot = [y(4:6); a1; y(10:12); a2; y(16:18); a3];
end

在这个程序中，我们首先定义了一些参数，如引力常数、太阳、地球和月球的质量、太阳、地球和月球的半径以及月球与地球之间的距离。接下来，我们定义了初始条件，包括地球和月球的位置和速度，并将它们合并为一个状态向量。然后，我们使用欧拉法来计算ODE系统的右手边，并在一定时间范围内进行积分。最后，我们将结果绘制在一个图形中。

请注意，欧拉法不是最准确的方法来计算ODE系统，因为它的误差会随着时间的增加而增加。更好的方法是使用更高阶的数值积分方法，如四阶龙格-库塔法（RK4）。