Problems with simulation of Earth-Moon trajectories around the Sun

Question

Julien Barrat 2017-5-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/340883-problems-with-simulation-of-earth-moon-trajectories-around-the-sun

评论： Sara 2022-10-30

Hi everybody! My name is Julien, I'm new on the forum :) I am trying (for a homework in computational physics) to plot the trajectories of the Earth and of the Moon around the Sun in a xy-plane. The Sun is a fixed point in (0,0), therefore I end up having 4 coupled differential equations (2 for the Earth, 2 for the Moon) with initial conditions, and I am trying to solve them numerically using ode45 (which I have never used before). I read the documentation, but something is still not working: my Earth and Moon are strangely flying away from the Sun in an (almost) straight line. Would you mind taking a look at my script and telling me if this is a conceptual problem (maybe the equations are wrong?) or if it is a programming issue?

The equations are pure Newton's 2nd law and are indicated in the code. I have a feeling that the problem might lie in an incorrect handling of the 2nd order equations in my definition of f. Please let me know what you guys think.

close all, clear all
G=6.67408e-11; %gravitational constant in m^3/kg/s^2 [CODATA]
m=[1.989e30... %mass of Sun (kg)
   5.972e24... %mass of Earth (kg)
   7.349e22]; %mass of Moon (kg)
r0=[1.496e11 0;... %initial position of Earth (m)
    1.496e11 -3.844e8]; %initial position of Moon (m)
v0=[0 29.78e3;... %initial velocity of Earth (m/s)
    1.023e3 0]; %initial velocities of Moon (m/s)
ic=[r0(1,1);r0(1,2);v0(1,1);v0(1,2);... %initial conditions for Earth
    r0(2,1);r0(2,2);v0(2,1);v0(2,2)]; %initial conditions for Moon
f=@(t,x) [x(3);-G*(m(1)*x(1)/(x(1)^2+x(2)^2)^(3/2)+m(3)*(x(1)-x(5))/((x(1)-x(5))^2+(x(2)-x(6))^2)^(3/2));... %equation for Earth (x)
        x(4);-G*(m(1)*x(2)/(x(1)^2+x(2)^2)^(3/2)+m(3)*(x(2)-x(6))/((x(1)-x(5))^2+(x(2)-x(6))^2)^(3/2));... %equation for Earth (y)
        x(7);-G*(m(1)*x(5)/(x(5)^2+x(6)^2)^(3/2)+m(2)*(x(5)-x(1))/((x(5)-x(1))^2+(x(6)-x(2))^2)^(3/2));... %equation for Moon (x)
        x(8);-G*(m(1)*x(6)/(x(5)^2+x(6)^2)^(3/2)+m(2)*(x(6)-x(2))/((x(5)-x(1))^2+(x(6)-x(2))^2)^(3/2))]; %equation for Moon (y)
tspan=0:1000:10000; %time span (s)
[t,x]=ode45(f,tspan,ic); %solve coupled differential equations with initial conditions
figure
plot(0,0,'o','Color',[0.8 0.8 0]) %plots the Sun at (0,0)
hold on
plot(x(:,1),x(:,2),'o','Color','blue') %plots the Earth
hold on
plot(x(:,5),x(:,6),'o','Color',[0.3 0.3 0.3]); %plots the Moon
xlim([-2.5e11 2.5e11])
ylim([-2.5e11 2.5e11])

The problem is that the y-coordinates of both the Earth (x(:,2)) and the Moon (x(:6)) are barely changing, while the velocities in x-direction (x(:,3) for the Earth, x(:,7) for the Moon) evolve rapidly!

Thank you very much in advance for your answers.

Julien.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Julien Barrat 2017-5-18

Okay I think I have circled the problem. It seems that because I did (unwittingly) a substitution for the 2nd order of the differential equation, I had to rewrite the equation differently (and possibly add an equation for each coordinate).

I am working on a fix at the moment, I still would appreciate advices if someone comes around here.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

James Tursa 2017-5-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/340883-problems-with-simulation-of-earth-moon-trajectories-around-the-sun#answer_267533

在 MATLAB Online 中打开

I haven't looked at your derivative equations yet, but the order is definitely wrong. E.g., take the first two elements:

ic=[r0(1,1);r0(1,2);v0(1,1);v0(1,2);... %initial conditions for Earth
    r0(2,1);r0(2,2);v0(2,1);v0(2,2)]; %initial conditions for Moon

Clearly, the first two elements of your state vector are the position elements for the Earth. But then in your derivative function:

f=@(t,x) [x(3);-G*(m(1)*x(1)/(x(1)^2+x(2)^2)^(3/2)+m(3)*(x(1)-x(5))/((x(1)-x(5))^2+(x(2)-x(6))^2)^(3/2));... %equation for Earth (x)
        x(4);-G*(m(1)*x(2)/(x(1)^2+x(2)^2)^(3/2)+m(3)*(x(2)-x(6))/((x(1)-x(5))^2+(x(2)-x(6))^2)^(3/2));... %equation for Earth (y)
        x(7);-G*(m(1)*x(5)/(x(5)^2+x(6)^2)^(3/2)+m(2)*(x(5)-x(1))/((x(5)-x(1))^2+(x(6)-x(2))^2)^(3/2));... %equation for Moon (x)
        x(8);-G*(m(1)*x(6)/(x(5)^2+x(6)^2)^(3/2)+m(2)*(x(6)-x(2))/((x(5)-x(1))^2+(x(6)-x(2))^2)^(3/2))]; %equation for Moon (y)

in that 2nd spot you have the velocity derivative, not the position derivative. This does not match your state vector ic. So at the very least, change the order of your derivative. E.g.,

f=@(t,x) [x(3);...
           x(4);...
          -G*(m(1)*x(1)/(x(1)^2+x(2)^2)^(3/2)+m(3)*(x(1)-x(5))/((x(1)-x(5))^2+(x(2)-x(6))^2)^(3/2));... %equation for Earth (x)
          -G*(m(1)*x(2)/(x(1)^2+x(2)^2)^(3/2)+m(3)*(x(2)-x(6))/((x(1)-x(5))^2+(x(2)-x(6))^2)^(3/2));... %equation for Earth (y)
         x(7);...
         x(8);...
        -G*(m(1)*x(5)/(x(5)^2+x(6)^2)^(3/2)+m(2)*(x(5)-x(1))/((x(5)-x(1))^2+(x(6)-x(2))^2)^(3/2));... %equation for Moon (x)
        -G*(m(1)*x(6)/(x(5)^2+x(6)^2)^(3/2)+m(2)*(x(6)-x(2))/((x(5)-x(1))^2+(x(6)-x(2))^2)^(3/2))]; %equation for Moon (y)

And then as a side comment I would suggest you double check your initial conditions for the moon velocity. Seems to me that it should be close to the earth velocity but that is not the case in your ic vector.

13 个评论
显示 11更早的评论隐藏 11更早的评论

Julien Barrat 2017-5-18

编辑：Julien Barrat 2017-5-18

在 MATLAB Online 中打开

You are right, the velocity of the Moon was meant relatively to the Earth. I didn't realise that until now.

Thank you both for your comments, I have corrected the order of the initial conditions and separated the equations in a function to make it easier to edit. Now I have as a main script:

close all, clear all
G=6.67408e-11; %gravitational constant in m^3/kg/s^2 [CODATA]
m=[1.989e30... %mass of Sun (kg)
   5.972e24... %mass of Earth (kg)
   7.349e22]; %mass of Moon (kg)
r0=[1.496e11 0;... %initial position of Earth (m)
    1.496e11 -3.844e8]; %initial position of Moon (m)
v0=[0 29.78e3;... %initial velocity of Earth (m/s)
    1.023e3 29.78e3]; %initial velocities of Moon (m/s)
ic=[r0(1,1);...
    v0(1,1);...
    r0(1,2);...
    v0(1,2);...
    r0(2,1);...
    v0(2,1);...
    r0(2,2);...
    v0(2,2)];
tspan=0:10:100; %time span (s)
[t,x]=ode45(@(t,y) odesolver(t,y),tspan,ic); %solve coupled differential equations with initial conditions
figure
plot(0,0,'o','Color',[0.8 0.8 0]) %plots the Sun at (0,0)
hold on
plot(x(:,1),x(:,3),'-','Color','blue') %plots the Earth
hold on
plot(x(:,5),x(:,7),'-','Color',[0.3 0.3 0.3]); %plots the Moon
xlim([-2.5e11 2.5e11])
ylim([-2.5e11 2.5e11])

and my function is:

function [ dydt ] = odesolver( t,y )
dydt=zeros(8,1);
G=6.67408e-11; %gravitational constant in m^3/kg/s^2 [CODATA]
m=[1.989e30... %mass of Sun (kg)
    5.972e24... %mass of Earth (kg)
    7.349e22]; %mass of Moon (kg)
dydt(1)=y(2); %Eq. of substitution (Earth,x)
dydt(2)=-G*(m(1)*y(1)/(y(1)^2+y(3)^2)^(3/2)+m(3)*(y(1)-y(5))/((y(1)-y(5))^2+(y(3)-y(7))^2)^(3/2)); %Eq, for Earth (x)
dydt(3)=y(4); %Eq. of substitution (Earth,y)
dydt(4)=-G*(m(1)*y(3)/(y(1)^2+y(3)^2)^(3/2)+m(3)*(y(3)-y(7))/((y(1)-y(5))^2+(y(3)-y(7))^2)^(3/2)); %Eq, for Earth (y)
dydt(5)=y(6); %Eq. of substitution (Moon,x)
dydt(6)=-G*(m(1)*y(5)/(y(5)^2+y(7)^2)^(3/2)+m(3)*(y(5)-y(1))/((y(5)-y(1))^2+(y(7)-y(3))^2)^(3/2)); %Eq, for Moon (x)
dydt(7)=y(8); %Eq. of substitution (Moon,y)
dydt(8)=-G*(m(1)*y(7)/(y(5)^2+y(7)^2)^(3/2)+m(3)*(y(7)-y(3))/((y(5)-y(1))^2+(y(7)-y(3))^2)^(3/2)); %Eq, for Moon (y)
end

However, this is still not working. Now both Earth and Moon are flying away in the -y-direction... Maybe the equations are wrong in the 1st place, I am gonna check that out.

Candy Reigo 2019-10-21

Hello Julien :) would you please upload the file once again. I am having problems with solving sun-moon-earth systems. It would be of a great help if you could upload your code again :) thank you.

Sara 2022-10-30

Hi

I know it's a long time ago but can you upload the file again?

thank you!

请先登录，再进行评论。

Problems with simulation of Earth-Moon trajectories around the Sun

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

13 个评论
显示 11更早的评论隐藏 11更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Problems with simulation of Earth-Moon trajectories around the Sun

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

13 个评论 显示 11更早的评论隐藏 11更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

13 个评论
显示 11更早的评论隐藏 11更早的评论