# Three-body problem question

1. Jul 28, 2012

### Dmalyuta

Hi all,

I'm using Matlab to try and smulate the orbits of Mercury and Venus aroudn the Sun, allowing all three to disturbe each others' orbits.

However, I am constantly running into this problem- all of my orbits look like in the pictures below- in effect, the planets spin out of orbit and will not stay in one ellipse. Even when I Venus out of the equation, mercury still displays this wild behavior.

This is obviosly the result of a problem in the code I've written for this simulation; my question si whether someone could explain to me for what reason such bizarre motion would result (what force is missing, etc.?) so that I can start looking in the right places in my code. For reference, the code I wrote is as follows. For only Sun and Mercury:

Code (Text):
% Mapping planet coordinates for 3 bodies- Sun, Mercury, Venus

clear
clf

%% Asking user to input some information and defining some constants

tim=10000;
interva=0.3;

G=6.67384e-11;

%% Sun initial coordinates

M_S=1.989e+30; % mass of sun
x(1,1)=0.0; % x-coordinate. (a,b) where a is the planet number (Sun is #1) and b is the time interval number (1 means time=0s)
y(1,1)=0.0;
V_x(1,1)=0.0; % velocity in x-direction
V_y(1,1)=0.0; % velocity in y-direction

%% Mercury initial coordinates

M_M=3.3e+23; % mass of mercury
x(2,1)=6.98e+10; % note that here a=2, which is mercury's number from the Sun. x-coordinate derived from planet's aphelion
y(2,1)=0.0; % all planets in my simulation start from x-axis
V_x(2,1)=0.0;
V_y(2,1)=38860; % note that all planets start from minimum velocity as such is the case at the aphelion (furthest) position from the sun

%% Venus initial coordinates

M_V=4.87e+24; % mass of venus
x(3,1)=1.09e+11;
y(3,1)=0.0;
V_x(3,1)=0.0;
V_y(3,1)=34790;

%% We now run a loop which will drop out coordiantes of planet positions at different times

timeend=tim*24*3600;
interval=interva*24*3600;

cell=1;

for iter=0:interval:timeend

for i=1:2

switch i

case 1

r_SM=sqrt((x(2,cell)-x(1,cell))^2+(y(2,cell)-y(1,cell))^2);
r_SV=sqrt((x(3,cell)-x(1,cell))^2+(y(3,cell)-y(1,cell))^2);

a_x(1,cell)=-(G*M_M*(x(1,cell)-x(2,cell)))/(r_SM^3);
a_y(1,cell)=-(G*M_M*(y(1,cell)-y(2,cell)))/(r_SM^3);

V_x(1,cell+1)=V_x(1,cell)+interval*a_x(1,cell);
V_y(1,cell+1)=V_y(1,cell)+interval*a_y(1,cell);

x(1,cell+1)=x(1,cell)+interval*V_x(1,cell);
y(1,cell+1)=y(1,cell)+interval*V_y(1,cell);

case 2

r_SM=sqrt((x(2,cell)-x(1,cell))^2+(y(2,cell)-y(1,cell))^2);
r_MV=sqrt((x(3,cell)-x(2,cell))^2+(y(3,cell)-y(2,cell))^2);

a_x(2,cell)=-(G*M_S*(x(2,cell)-x(1,cell)))/(r_SM^3);
a_y(2,cell)=-(G*M_S*(y(2,cell)-y(1,cell)))/(r_SM^3);

V_x(2,cell+1)=V_x(2,cell)+interval*a_x(2,cell);
V_y(2,cell+1)=V_y(2,cell)+interval*a_y(2,cell);

x(2,cell+1)=x(2,cell)+interval*V_x(2,cell);
y(2,cell+1)=y(2,cell)+interval*V_y(2,cell);

end

end

cell=cell+1;

end

hold on

for i=1:2

plot(x(i,:),y(i,:))

end

hold off
For Sun, Mercury, and Venus:

Code (Text):
% Mapping planet coordinates for 3 bodies- Sun, Mercury, Venus

clear
clf

%% Asking user to input some information and defining some constants

tim=10000;
interva=0.3;

G=6.67384e-11;

%% Sun initial coordinates

M_S=1.989e+30; % mass of sun
x(1,1)=0.0; % x-coordinate. (a,b) where a is the planet number (Sun is #1) and b is the time interval number (1 means time=0s)
y(1,1)=0.0;
V_x(1,1)=0.0; % velocity in x-direction
V_y(1,1)=0.0; % velocity in y-direction

%% Mercury initial coordinates

M_M=3.3e+23; % mass of mercury
x(2,1)=6.98e+10; % note that here a=2, which is mercury's number from the Sun. x-coordinate derived from planet's aphelion
y(2,1)=0.0; % all planets in my simulation start from x-axis
V_x(2,1)=0.0;
V_y(2,1)=38860; % note that all planets start from minimum velocity as such is the case at the aphelion (furthest) position from the sun

%% Venus initial coordinates

M_V=4.87e+24; % mass of venus
x(3,1)=1.09e+11;
y(3,1)=0.0;
V_x(3,1)=0.0;
V_y(3,1)=34790;

%% We now run a loop which will drop out coordiantes of planet positions at different times

timeend=tim*24*3600;
interval=interva*24*3600;

cell=1;

for iter=0:interval:timeend

for i=1:3

switch i

case 1

r_SM=sqrt((x(2,cell)-x(1,cell))^2+(y(2,cell)-y(1,cell))^2);
r_SV=sqrt((x(3,cell)-x(1,cell))^2+(y(3,cell)-y(1,cell))^2);

a_x(1,cell)=-(G*M_M*(x(1,cell)-x(2,cell)))/(r_SM^3)+-(G*M_V*(x(1,cell)-x(3,cell)))/(r_SV^3);
a_y(1,cell)=-(G*M_M*(y(1,cell)-y(2,cell)))/(r_SM^3)+-(G*M_V*(y(1,cell)-y(3,cell)))/(r_SV^3);

V_x(1,cell+1)=V_x(1,cell)+interval*a_x(1,cell);
V_y(1,cell+1)=V_y(1,cell)+interval*a_y(1,cell);

x(1,cell+1)=x(1,cell)+interval*V_x(1,cell);
y(1,cell+1)=y(1,cell)+interval*V_y(1,cell);

case 2

r_SM=sqrt((x(2,cell)-x(1,cell))^2+(y(2,cell)-y(1,cell))^2);
r_MV=sqrt((x(3,cell)-x(2,cell))^2+(y(3,cell)-y(2,cell))^2);

a_x(2,cell)=-(G*M_S*(x(2,cell)-x(1,cell)))/(r_SM^3)+-(G*M_V*(x(2,cell)-x(3,cell)))/(r_SV^3);
a_y(2,cell)=-(G*M_S*(y(2,cell)-y(1,cell)))/(r_SM^3)+-(G*M_V*(y(2,cell)-y(3,cell)))/(r_SV^3);

V_x(2,cell+1)=V_x(2,cell)+interval*a_x(2,cell);
V_y(2,cell+1)=V_y(2,cell)+interval*a_y(2,cell);

x(2,cell+1)=x(2,cell)+interval*V_x(2,cell);
y(2,cell+1)=y(2,cell)+interval*V_y(2,cell);

case 3

r_SV=sqrt((x(3,cell)-x(1,cell))^2+(y(3,cell)-y(1,cell))^2);
r_MV=sqrt((x(3,cell)-x(2,cell))^2+(y(3,cell)-y(2,cell))^2);

a_x(3,cell)=-(G*M_S*(x(3,cell)-x(1,cell)))/(r_SM^3)+-(G*M_M*(x(3,cell)-x(2,cell)))/(r_SV^3);
a_y(3,cell)=-(G*M_S*(y(3,cell)-y(1,cell)))/(r_SM^3)+-(G*M_M*(y(3,cell)-y(2,cell)))/(r_SV^3);

V_x(3,cell+1)=V_x(3,cell)+interval*a_x(3,cell);
V_y(3,cell+1)=V_y(3,cell)+interval*a_y(3,cell);

x(3,cell+1)=x(3,cell)+interval*V_x(3,cell);
y(3,cell+1)=y(3,cell)+interval*V_y(3,cell);

end

end

cell=cell+1;

end

hold on

plot(x(1,:),y(1,:),'+')
plot(x(2,:),y(2,:),'r')
plot(x(3,:),y(3,:),'b')

hold off

2. Jul 28, 2012

### Filip Larsen

Welcome to PF!

The acceleration you calculate in each step should not be accumulated, that is, you should have "=" and not "=-" for the lines that calculate acceleration.

In addition, you are using the very simple Euler integration [1] which will introduce false energy into the system. A better approach would be to use Leapfrog or a something method [2].

I have not made a detailed examination of your code, so I may easily have missed any other errors.

[1] http://en.wikipedia.org/wiki/Euler_method
[2] http://en.wikipedia.org/wiki/Leapfrog_method

3. Jul 28, 2012

### D H

Staff Emeritus
That's not the problem. This is not accumulating. That would be -=. His "=-" is the same as "acc = -(computation)".

That's not the problem, either. He's using symplectic Euler, not Euler's method. Admittedly there are far better techniques, but symplectic Euler is an okay start. (Aside: Leapfrog is barely a step above symplectic Euler. There are far, far better.)

The problem is that the accelerations are being computed sequentially. The accelerations need to be computed all at once, followed by state propagation.

4. Jul 28, 2012

### Dmalyuta

Could you please use simpler language in explaining that last problem about calculating accelerations all at once? I'm just trying to make sense of why the heck the planets would constantly be gaining energy and be moving into a higher orbit- there is nothing in my code which would imply that. Otherwise, it might be a problem of me loosing some force when the planets go below the x-axis; but then again, where in my code am I loosing that force??

PS. The code is extrapolated from Feynman's "Mechanics 1"- he described a numerical calculation for the sun and a planet (2 bodies) and touched upon 3 bodies (which is reflected in my code).

5. Jul 28, 2012

### AlephZero

You need to calculate all the accelerations before you update the position of the objects.

Code (Text):

r_SM=sqrt((x(2,cell)-x(1,cell))^2+(y(2,cell)-y(1,cell))^2);
r_SV=sqrt((x(3,cell)-x(1,cell))^2+(y(3,cell)-y(1,cell))^2);

a_x(1,cell)=-(G*M_M*(x(1,cell)-x(2,cell)))/(r_SM^3)+-(G*M_V*(x(1,cell)-x(3,cell)))/(r_SV^3);
a_y(1,cell)=-(G*M_M*(y(1,cell)-y(2,cell)))/(r_SM^3)+-(G*M_V*(y(1,cell)-y(3,cell)))/(r_SV^3);

r_SM=sqrt((x(2,cell)-x(1,cell))^2+(y(2,cell)-y(1,cell))^2);
r_MV=sqrt((x(3,cell)-x(2,cell))^2+(y(3,cell)-y(2,cell))^2);

a_x(2,cell)=-(G*M_S*(x(2,cell)-x(1,cell)))/(r_SM^3)+-(G*M_V*(x(2,cell)-x(3,cell)))/(r_SV^3);
a_y(2,cell)=-(G*M_S*(y(2,cell)-y(1,cell)))/(r_SM^3)+-(G*M_V*(y(2,cell)-y(3,cell)))/(r_SV^3);

r_SV=sqrt((x(3,cell)-x(1,cell))^2+(y(3,cell)-y(1,cell))^2);
r_MV=sqrt((x(3,cell)-x(2,cell))^2+(y(3,cell)-y(2,cell))^2);

a_x(3,cell)=-(G*M_S*(x(3,cell)-x(1,cell)))/(r_SM^3)+-(G*M_M*(x(3,cell)-x(2,cell)))/(r_SV^3);
a_y(3,cell)=-(G*M_S*(y(3,cell)-y(1,cell)))/(r_SM^3)+-(G*M_M*(y(3,cell)-y(2,cell)))/(r_SV^3);

V_x(1,cell+1)=V_x(1,cell)+interval*a_x(1,cell);
V_y(1,cell+1)=V_y(1,cell)+interval*a_y(1,cell);

V_x(2,cell+1)=V_x(2,cell)+interval*a_x(2,cell);
V_y(2,cell+1)=V_y(2,cell)+interval*a_y(2,cell);

V_x(3,cell+1)=V_x(3,cell)+interval*a_x(3,cell);
V_y(3,cell+1)=V_y(3,cell)+interval*a_y(3,cell);

x(1,cell+1)=x(1,cell)+interval*V_x(1,cell);
y(1,cell+1)=y(1,cell)+interval*V_y(1,cell);

x(2,cell+1)=x(2,cell)+interval*V_x(2,cell);
y(2,cell+1)=y(2,cell)+interval*V_y(2,cell);

x(3,cell+1)=x(3,cell)+interval*V_x(3,cell);
y(3,cell+1)=y(3,cell)+interval*V_y(3,cell);

(I don't think your "for i=1:3", "switch i", and "case" statements do anything useful here).

6. Jul 28, 2012

### Filip Larsen

Good thing you are more awake than I apparently was Hope the OP didn't get too confused about my wild shots there ...

7. Jul 28, 2012

### haruspex

I don't think that's it either. They are effectively calculated in parallel. The new acceleration is stored in cell+1. This isn't used until the next time around the loop.

It could be just accumulation of error through using finite intervals. Is it better if you reduce the interval size? That's always a useful test - reduce the interval size until it stops having a discernible effect on the result. A possible improvement is to use the average of the old and new accelerations when updating the velocities, and the average of old and new velocities when updating positions. This is roughly equivalent to reducing the interval size a bit.

8. Jul 29, 2012

### Dmalyuta

haruspex, I don't think this is accumulation of error as I used the same method to simulate the motion of just the Sun and mercury without Venus and Mercury's orbit was perfectly fine.

I think we should look into the signs of my calculated forces/accelerations in each quadrant when the x and y-coordinate signs change. I tried but I am blinded by the fact that I wrote the code, so it's harder for me to fidn the error than if somebody took a shot at it...

9. Jul 29, 2012

### D H

Staff Emeritus
You're right, and so was Filip Larsen in post #2. The OP is using plain old Euler, which spirals out.

Even if the OP used the best integrator in the world, the ellipses would still move. The reason is that the Sun is starting with zero velocity. This makes the center of mass of the system as a whole have a non-zero velocity, about 6.45 milimeters per second, or 203 km/year, in the OP's first example of the Sun and Mercury.

10. Jul 29, 2012

### Dmalyuta

I re-wrote the code, and am still getting this bizarre effect. Here is the new code:

Code (Text):
%% The Three Body Problem

clear
clf

%% Asking user to input some information and defining some constants

tim=10000;
interva=0.3;

G=6.67384e-11;

%% Sun initial coordinates

M(1)=1.989e+30; % mass of sun
x(1,1)=0.0; % x-coordinate. (a,b) where a is the planet number (Sun is #1) and b is the time interval number (1 means time=0s)
y(1,1)=0.0;
V_x(1,1)=0.0; % velocity in x-direction
V_y(1,1)=0.0; % velocity in y-direction

%% Mercury initial coordinates

M(2)=3.3e+23; % mass of mercury
x(2,1)=6.98e+10; % note that here a=2, which is mercury's number from the Sun. x-coordinate derived from planet's aphelion
y(2,1)=0.0; % all planets in my simulation start from x-axis
V_x(2,1)=0.0;
V_y(2,1)=38860; % note that all planets start from minimum velocity as such is the case at the aphelion (furthest) position from the sun

%% Venus initial coordinates

M(3)=4.87e+24; % mass of venus
x(3,1)=1.09e+11;
y(3,1)=0.0;
V_x(3,1)=0.0;
V_y(3,1)=34790;

%% We now run a loop which will drop out coordinates of planet positions at different times

timeend=tim*24*3600;
interval=interva*24*3600;

cell=1;

for time=0:interval:timeend

for i=1:3

switch i

case 1

% Calculate the next x-coordinate

r(i,2,cell)=sqrt((x(i,cell)-x(2,cell))^2+(y(i,cell)-y(2,cell))^2);
r(i,3,cell)=sqrt((x(i,cell)-x(3,cell))^2+(y(i,cell)-y(3,cell))^2);

a_x(i,cell)=-(G*M(2)*(x(i,cell)-x(2,cell)))/(r(i,2,cell)^3)+-(G*M(3)*(x(i,cell)-x(3,cell)))/(r(i,3,cell)^3);

V_x(i,cell+1)=V_x(i,cell)+interval*a_x(i,cell);

x(i,cell+1)=x(i,cell)+interval*V_x(i,cell);

% Calculate the next y-coordinate

a_y(i,cell)=-(G*M(2)*(y(i,cell)-y(2,cell)))/(r(i,2,cell)^3)+-(G*M(3)*(y(i,cell)-y(3,cell)))/(r(i,3,cell)^3);

V_y(i,cell+1)=V_y(i,cell)+interval*a_y(i,cell);

y(i,cell+1)=y(i,cell)+interval*V_y(i,cell);

case 2

% Calculate the next x-coordinate

r(i,1,cell)=sqrt((x(i,cell)-x(1,cell))^2+(y(i,cell)-y(1,cell))^2);
r(i,3,cell)=sqrt((x(i,cell)-x(3,cell))^2+(y(i,cell)-y(3,cell))^2);

a_x(i,cell)=-(G*M(1)*(x(i,cell)-x(1,cell)))/(r(i,1,cell)^3)+-(G*M(3)*(x(i,cell)-x(3,cell)))/(r(i,3,cell)^3);

V_x(i,cell+1)=V_x(i,cell)+interval*a_x(i,cell);

x(i,cell+1)=x(i,cell)+interval*V_x(i,cell);

% Calculate the next y-coordinate

a_y(i,cell)=-(G*M(1)*(y(i,cell)-y(1,cell)))/(r(i,1,cell)^3)+-(G*M(3)*(y(i,cell)-y(3,cell)))/(r(i,3,cell)^3);

V_y(i,cell+1)=V_y(i,cell)+interval*a_y(i,cell);

y(i,cell+1)=y(i,cell)+interval*V_y(i,cell);

case 3

% Calculate the next x-coordinate

r(i,2,cell)=sqrt((x(i,cell)-x(2,cell))^2+(y(i,cell)-y(2,cell))^2);
r(i,1,cell)=sqrt((x(i,cell)-x(1,cell))^2+(y(i,cell)-y(1,cell))^2);

a_x(i,cell)=-(G*M(2)*(x(i,cell)-x(2,cell)))/(r(i,2,cell)^3)+-(G*M(1)*(x(i,cell)-x(1,cell)))/(r(i,1,cell)^3);

V_x(i,cell+1)=V_x(i,cell)+interval*a_x(i,cell);

x(i,cell+1)=x(i,cell)+interval*V_x(i,cell);

% Calculate the next y-coordinate

a_y(i,cell)=-(G*M(2)*(y(i,cell)-y(2,cell)))/(r(i,2,cell)^3)+-(G*M(1)*(y(i,cell)-y(1,cell)))/(r(i,1,cell)^3);

V_y(i,cell+1)=V_y(i,cell)+interval*a_y(i,cell);

y(i,cell+1)=y(i,cell)+interval*V_y(i,cell);

end

end

cell=cell+1;

end

hold on

plot(x(1,:),y(1,:),'+');
plot(x(2,:),y(2,:),'r');
plot(x(3,:),y(3,:),'b');

hold off
And here is a picture of the resulting orbits:

I am thinking that perhaps the planets spin out of orbit because the sun is allowed to move around? If you zoom into the Sun, you get this pattern for its motion:

Note how much it moves (7 orders of magnitude in the y-direction!!!); but then again, the real sun is allowed to move and the planets don't spin out of orbit- so what's the problem!?!?! This is really starting to annoy me....

11. Jul 29, 2012

### nasu

The problem may be that your initial conditions are such that the total momentum of the system is not zero. The center of mass is in motion. If you want nice closed trajectories you need to calculate in the center of mass system or arrange the initial conditions so that the speed of the CM is zero in the system of reference used.
Maybe you do calculate in the CM? I did not go through your code in detail.

12. Jul 29, 2012

### TurtleMeister

Yes, you are correct. However, that will cause the simulation to creep, not spiral. This effect is much more noticable when doing a simulation of the earth moon system. When I do a simulation of the earth moon system I must give the moon a positive initial velocity and the earth a negative initial velocity, otherwise it will creep. All bodies should have velocities relatve to the center of mass of the system (as nasu pointed out).

I cannot find a problem with this code. however, I am not very familiar with Matlab. So I don't know what's causing this. It does look like an accumulative problem, but I don't see it in the code. I do agree with AlephZero, that your use of the "I" loop is ineffective. However, I don't think that has anything to do with this problem. But it may be something you will want to fix if you decide to use more n-bodies.

13. Jul 29, 2012

### D H

Staff Emeritus
You are still using Euler's method, and you are still intermixing calculations of accelerations and using them for propagation, and you probably still have the Sun starting with zero velocity.

Don't use Euler's method. It will give bad results. It is important to learn Euler's method because this method is the basis for many integration techniques. However, because the results are inevitably poor, it is also important to learn that it is best to not use this method.

Don't intermingle the calculations of accelerations and using them for propagation. It's not hurting you here, but that's only because you are using Euler's method. You can't get away with this for more advanced techniques. Calculate all of your accelerations before you do anything with them. One good way to do this is to put your acceleration computations (all of them) in a separate subroutine.

Don't make the Sun's initial velocity zero. You need to make the velocity of the system center of mass zero. If you don't the system as a whole will drift with this velocity.

Your wide images are starting to annoy me. Scale them down to a reasonable width. 800 pixels or 640 pixels are a good values. Don't go wider than 800, please.

Last edited: Jul 29, 2012
14. Jul 29, 2012

### Dmalyuta

Thank you for the help. I took a shot at something, and now I feel stupid for having spent so much of your and my own time on this spiral dilemma- the problem was in, like many of you said, the method i was using and, in particular, the interval I was using. With 0.3 days, I get the spiral effect like in the pictures in above posts; I now changed to 0.001 days, and now I get this graph:

Now, my next question is- which method do you suggest I use which will give me a much more accurate result with an interval larger than 0.001 days (which takes a long time for my computer to compute!!)?

Thank you!

15. Jul 29, 2012

### D H

Staff Emeritus
First off, did you try symplectic Euler? In a nutshell, you use the updated velocity to propagate position rather than the pre-update velocity. It seems wrong, but it has a solid mathematical basis.

Filip gave you pointers to leapfrog, which has a number of variants. Another simple technique is velocity verlet. Leapfrog and verlet integration are often used in simulations of a large number (thousands or more) of bodies. These techniques aren't used so much because they're good; they're only good when compared to Euler / symplectic Euler. They're used because of how they come close to conserving energy.

There's an almost inescapable tradeoff here between long-term and short-term accuracy. Energy and angular momentum must be conserved to attain long-term accuracy (many million years for solar system simulations, or much shorter span for galactic-scale simulations). However, the techniques that do conserve energy / angular momentum are low order and only conserve these quantities approximately. (Better said: They preserve something that is a close analog to energy.)

If you want short-term accuracy you need to use higher order techniques. (Note that short term can mean millions of years with a good technique.) A good place to start is fourth order Runge Kutta. This technique takes four internal steps to achieve the desired external step size. Another approach is to take advantage of accelerations and velocities from previous time steps as well as the current step. This leads to the Adams-Bashforth-Moulton family of techniques. Gauss Jackson is very similar to the ABM family, but it is very hairy. Yet another approach is to specify a very long time step and let the integration technique split this up however it sees fit. You use Matlab; this is exactly what ode45() does.

A very different approach is to integrate something other than cartesian position and velocity. The orbits of the planets are very close to being ellipses. So propagate a reference ellipse, and compute your actual state as a perturbation from this. This approach was used with the Apollo missions. An alternative is to represent state as orbital elements rather than position and velocity and compute the time derivatives of these orbital elements. These derivatives are called the planetary equations of motion. There are a number of these, depending on which orbital elements one is using. A similar approach is to use quaternions to represent planetary state via the Karhunen–Loève transform.

16. Jul 29, 2012

### mikeph

I came across this problem a while ago and ended up using the methods of Quinlan and Tremaine (1990). You basically find the new position as a function of the previous n positions/velocities/accelerations. A good approach is to generate the first few positions using a leapfrog method with very small timestep, then switch to a symmetric multistep method once you've gathered enough consecutive positions, increasing the length of the timestep to make it run faster.

They generally seem to conserve energy fairly well, and are fairly easy to code.

http://adsabs.harvard.edu/full/1990AJ....100.1694Q

17. Jul 31, 2012

### DrStupid

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook