Why Are My 3-Body Problem Orbits Merging Instead of Separating?

[John O' Meara]

I am writing a program in a version of Basic for showing the orbits of 3 bodies. To check that it works I do the following. Each body is of the same mass say 50. I cannot get them to rotate in three separate orbits of equal size, which is the test that the program works.
I place the bodies at the respective corners of an equilateral triangle so that the three bodies are an equal distance apart. The coordinates of the three bodies are:
x1,y1 = (-500,-288.675...),
x2,y2 = (500,-288.675...),
x3,y3 = (0,577.35...),
and their velocity components are:
Vx1,Vy1 = (20,-34.6...),
Vx2,Vy2 = (20,34.6...),
Vx3,Vy3 = (-40,0). This way the change of momentum of the entire system equals zero. The result that I get is that body-2's orbit coalesces with body-3's and three's coalesces with body-1's e.t.c., to give three bodies having the same large anticlockwise orbit, instead of three separate orbits. I know that the starting positions and directions of the initial velocity components are critical.While I have the initial positions correct I do not know if the directions ' are correct. If I put any two bodies in a clockwise orbit and the remaining body in a anticlockwise orbit, I seem not able to fulfill the requirement that the change in momentum of the entire system equals zero, yet I feel that this is required to get the system of masses orbiting correctly.
Help would be appreciated.
Note : In order for these magnitude of figures to give orbits at all I have set G=1.

 

[John O' Meara]

To see what I mean by coalesced and correct orbits please see word document attachment

 

[John O' Meara]

Sorry it did not upload the first time I attempted it.

 

[Ulysees]

Is it proper 4th order Runge Kutta integration you're doing? Or just Euler, which is unstable?

 

[John O' Meara]

After finding the component accelerations from Newton' Law, I first do the following:
Vx=Vx+ax*DT/2
Vy=Vy+ay*DT/2
On all subequent runs the following calculations is done:
Vx=Vx+ax*DT
Vy=Vy+ay*DT
where ax is the acceleration in the x direction, and Vx is the velocity in the x direction and x is the displacement in the x direction and DT is the time interval.
x=x+Vx*DT
y=y+Vy*DT
My problem is to get the correct orbits as indicated in the attachment. Do I need to put one body rotating counterclockwise and the other two rotating clockwise, if so how do I do that keeping the change in the system momentum zero? Because putting the three bodies in say a clockwise orbit gives my "coalesced" orbit see bottom diagram in attachment.

 

[Ulysees]

So it is Euler. You'll notice that if you just use two bodies, one heavy one light, so started that the light goes around the heavier, the ellipsis never closes. That's because of cumulative error.

Cumulative error that Euler does not prevent would be even worse in an unstable configuration, like a 3-body system in certain configurations. To be honest, I don't understand what you're trying to achieve, but almost certainly Euler is not good enough for it.

Why don't you start with reading a bit about Lagrange points, in the 3-body system, to get a taste of what instability you get? Is it a school exercise?

 

[John O' Meara]

It was a school exercise that I did years ago, and was originally a n-body problem, I don't have the school notes that I made then anymore. So I am working from memory. I know that it had a piece of program to close the loop and reduce the cumulative error.
What I am trying to achieve to do the n-body problem after I have got the correct motion of the three bodies.
I was thinking if I put one body rotating counterckockwise and the other two clockwise I should get three separate stable orbits as in the first diagram of the attachment. But I can't and that is my problem. A problem that I still have even if I was to use a non Euler method like a 4th order Runge Kutta for the calculation routine.

 

[Ulysees]

Sorry but you can't make them rotate in different directions on the same circle. This is because at any point the total force experienced by each object should match the curvature centre, direction, and speed. The one going in the wrong direction, does not satisfy this a moment after starting where it is a little closer to one of the others and a little further from the other. So it will go off the circle and mess up the others too.

But if they all go clockwise in the same circle and with equal distances between them and equal masses, then it's ok - except it's unstable, even proper integration will lose its balance after a while. Which is exactly what I have seen with a simulation program. Why don't you download it and play with it, to understand some stable and some unstable configurations.

One stable configuration, is the sun-earth-moon system. Put these masses into your parameters, put the right initial velocities and positions like the real thing, and away you go, you have a stable system. Stable in the software that is, where integration is right - with your Euler integration loops won't close exactly but it will still be reasonable. Here's the software:

http://www.orbitsimulator.com/gravity/articles/lagrange.html

 

[Ulysees]

> But if they all go clockwise in the same circle and with equal distances between them and equal masses, then it's ok - except it's unstable

Actually it's not ok, it can only work for certain combinations of masses, radius and angular velocity. I presume you have solved the equation and have a radius, speed and mass for each object that satisfies the equation:

centripetal acceleration = total gravity

for each of the 3 objects.

 

[John O' Meara]

What I wanted to know is just what you have said "But if they all go clockwise in the same circle and with equal distances between them and equal masses, then it's ok - except it's unstable". Thanks for the site and the tip about Lagrange points. When I have finished this mini-project I will study Lagrange.

 

[Ulysees]

Start with writing down the formula for the total gravity each object experiences. Have you done that? What is it?

Then imagine they have moved a little, but not all in the same direction along the circle. Is the total gravity the same for all?

It's a different vector for each. Therefore it's a different orbit for each object, nothing to do with a circle.

 

[John O' Meara]

Yes, for example the acceleration in the x direction is given in this procedure:
DEFPROCmotion(bx,by,x,y,Vx,Vy)
d=SQR((x-bx)^2+(y-by)^2)
IF (x-bx) <>0 THEN ax=-G*Mt*(x-bx)/d^1.5 ELSE ax=0
IF (y-by) <>0 THEN ay=-G*Mt*(y-by)?d^1.5 ELSE ay=0
IF n > 0 THEN GOTO ...
etc.(see post no 5)
Where x is the displacement of the body-1 under test in the x direction, and bx is the displacement ( in the x direction ) of the combined other two bodies who are acting on body-1 under their gravity. Mt ( which is a global variable) is the mass of the combined other two bodies.

 

[Ulysees]

This is fine. I meant work out numerically the acceleration in the second step and draw vectors on paper. It should be obvious the thing can't continue as it started in the 1st step because now the vectors are not all pointing to the centre, none of them is.

But if they all start clockwise, the vector still points to the centre. Unless you have shot the too fast or too slow. For circular motion the speed must satisfy this:

centrifugal accel = gravitational accel

=>

v^2 / r = GM / r^2

=>

v = sqrt( GM / r )

 

[John O' Meara]

I still cannot get the orbits out properly even though I now give a radial velocity component as accordingly to the last post. Please see the attachment to understand what my other posts on this subject were trying to say.

 

[Ulysees]

You write:
> The first three small circles are the correct orbits that I should get when three bodies are of equal mass are placed an equal distance appart

Three small circles? How on Earth did you come to expect this? Circular motion in each little circle requires a constant force towards the centre of each circle. You don't have such a force, you have a total force that points somewhere else.

Your code is correct but I recommend that you experiment with different speeds, watch what happens, and post a video of what you see on the screen. I hope you have written the code to show objects on the screen.

The speed formula I gave was incorrect, sorry (what's capital M doing there, you should have asked). The correct one comes from vector-addition of the two forces from the other two bodies. You seem to not have learned vector addition in this context so I say just experiment with different values of speed and show the results.

There exists one value of speed, and pointing to the right direction for each body, that will make them go round the centre of the triangle they form, but only for a while because the balance is unstable.

 

[John O' Meara]

I don't know how I expected the impossible with regard to the three small orbits. I have experimented with different speeds and have got one large circular orbit that is correct. I have only written code to show the paths on the screen, not the actual bodies as objects on the screen. At the moment I don't have a video recorder to tape the motion of the bodies.
Thanks for the reply.

 

[John O' Meara]

I had to give this a miss for nearly 8 months as I had not much time to myself. I have changed it to a n-body problem, but it does not work correctly. When I use it as a two-body problem e.g., M(1)=moon's mass=7.35E22, M(2)=Earth's mass=5.98E24, x(1)=-3.844E8, y(1)=0, x(2)=0, y(2)=0, Vx(1)=0, Vy(1)=1.023E3, Vx(2)=0, Vy(2)=-12.57, I get a horizontal line on the graphics screen with the moon oscillating back and forth with an increasing ampitude, not a circular orbit. In fact to get a circular orbit I must make Vy(1)=2.8E9, and Vy(2)=-3.44E7. All velocities are in m/s, coordinates are in metres, masses in kg. I will sumit the relevant parts of the program later today. Thanks for the help.

 

[John O' Meara]

I have had very little time to myself for the last 8 months, so I hadn't time to work on this thread. I have updated it to the n-body problem, lately. On testing the program with the two body problem it fails to work correctly. E.g., when I let M(1)= moon's mass=7.35E22, M(2)=Earth's mass=5.98E24, x(1)=-3.844E8 (which is the distance of the Moon from Earth), y(1)=0, x(2)=0, y(2)=0, Vx(1)=0, Vy(1)=1.023E3 ( the speed of the moon around the earth), Vx(2)=0, Vy(2)=-12.57; I get a horizontal line with the Moon oscillating back and forth with increasing amplitude per cycle. In fact to get a circular orbit I must make Vy(1)=2.8E9 and Vy(2)=-3.44E7, Which are values greater than the speed of light!
Where coordinates are in metres, [speed]=m/s, [mass]=kg. Please see attachments. Thanks for the help.

 

[John O' Meara]

Just a few additions to the attachments:line nos:
191 IF J=no THEN PROCtransfer
192 PROCinial

199 UNTIL orbitiscomplete