MATLAB: Problems with simulation of Earth-Moon trajectories around the Sun

differential equationsode45solar system

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.

Best Answer

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.
Related Question