Initial condition of the restricted problem of three bodies

[Heimdall]

Hi !


I'm trying to figure out what kind of initial conditions I should put in my program in order to have stable trajectories.

Just to remind, The restricted problem of three bodies suppose that you've got two heavy masses in circular orbit around their common center of mass, and a negligeable mass particule orbiting in that system.

well, now we know that let's go further. In the equations describing this system, you probably know that there is a constant called the jacobi's integral.

This constant is J = -2*K, where K his the http://nicolas.aunai.free.fr/cours/magistere/p3cr/cvn/toto.jpg )

you can see a little plot here : http://nicolas.aunai.free.fr/cours/magistere/p3cr/cvn/cvn.gif

representing those areas.

Well, since I'm looking for some stable trajectories, I thought that I could find a family of {Q1,Q2,P1,P2} variables which satisfies the relation J=-2K.

the problem is that I don't know how to do this... since my hamiltonian isn't trivial

is someone here got a idea ?

thx for all
bye

 

[pervect]

Hi !


I'm trying to figure out what kind of initial conditions I should put in my program in order to have stable trajectories.

Just to remind, The restricted problem of three bodies suppose that you've got two heavy masses in circular orbit around their common center of mass, and a negligeable mass particule orbiting in that system.

well, now we know that let's go further. In the equations describing this system, you probably know that there is a constant called the jacobi's integral.

This constant is J = -2*K, where K his the http://nicolas.aunai.free.fr/cours/magistere/p3cr/cvn/toto.jpg )

you can see a little plot here : http://nicolas.aunai.free.fr/cours/magistere/p3cr/cvn/cvn.gif

representing those areas.

Well, since I'm looking for some stable trajectories, I thought that I could find a family of {Q1,Q2,P1,P2} variables which satisfies the relation J=-2K.

the problem is that I don't know how to do this... since my hamiltonian isn't trivial

is someone here got a idea ?

thx for all
bye

I'm not sure what you mean by stable trajectories.

If you are looking for trajectories that don't leave a specific bounded region of space, take a look at

http://www.geocities.com/syzygy303/

For a mass ratio of u=.01, for J>3.2 and intial values "near" x=1, y=0 won't leave the vicinity of x=1, y=0. There are plots for various values of J near the interesting values of J which restrict the orbit of the thrid body to this region of space.

Check the formula for J used to make sure it matches yours, the variables may be different. Note that the "hamiltonian" in the variables used here is the numerically equivalent energy function - which is the Hamiltonian expressed in terms of (q, q') rather than (q,p).

The physical interpretation of the plots on the webpage is that it is a region of space near a "moon" that a body cannot leave - because it doesn't have enough of the conserved energy-like quantiy J to do so.

Now, you may also be interested in the Lagrange points.

The Lagrange points are equilibrium solutions with p=q=constant (so p'=q'=0) for a body. You can find them in principle from your Hamiltonian by solving Hamilton's equaitons


<br /> \dot{p} = -\frac{\partial H}{\partial q}<br />
<br /> \dot{q} = \frac{\partial H}{\partial p}<br />

For this particular problem, we are particularly interested in the locaiton of L1. As an equilibrium point, it's automaticaly an extremeum of the Jacobi intergal, just because it solves the Hamiltonian and x' and y' are zero because it's an equilibrium soluiton.

It turns out that L1 is the important point which marks the boundry between orbits that can never leave the vicinity of the moon, and orbits which can. Any orbit which can reach L1 can escape the vicinity of the moon.

Calculating the position of the Lagrange point L1 isn't possible in a closed form solution, but for small mass ratios u it's about 1-(u/3)^(1/3) where u is the mass ratio. You can find the positoin of L1 as a taylor serios in u^(1/3) if you really want to.

While all the Lagrange points are equilibrium solutions of Hamilton's equations, but only L4 and L5 are stable equilbria, and even then L4 and L5 are only stable for sufficiently high mass ratios. See

http://scienceworld.wolfram.com/physics/LagrangePoints.html

 

[remcook]

Do you have access to a good library?
In the book of Gomez et al. (2001) there's a list of initial conditions for halo orbits around L2.

Here it is...(It's a bit of a mess, because I copied it from a pdf table. only arranged the first few rows, but you get the drill)
x , z , dy/dt
Nr. 1 A 1.00841650141 0.0000000000000 0.0097914931461
B 1.01127600558 0.0000000000000 -0.0089501491171
Nr. 2 A 1.00841645317 -0.0000131712297 0.0097917163651
B 1.01127601837 0.0000165157448 -0.0089503226141
Nr. 3 A 1.00841506015 -0.0001130181527 0.0097964155892
B 1.01127582432 0.0001417315823 -0.0089533363648
Nr. 4 A 1.00522538832 -0.0046431089695 0.0193059074154
B 1.01033360395 0.0074555355146 -0.0135712737234
Nr. 5 A 1.00340906546 -0.0050249372803 0.0244921306580
B 1.00951124354 0.0095866938101 -0.0148198893936
Nr. 6 A 1.00187070402 -0.0047184790519 0.0296955121606
B 1.00860215012 0.0110624814649 -0.0152055165041
Nr. 7 A 1.00151252067 -0.0045380625232 0.0312210555816
B 1.00833004882 0.0113643518122 -0.0151705666816
Nr. 8 A 1.00140488602 -0.0044728053106 0.0317242727716
B 1.00824034396 0.0114517576939 -0.0151461497488
Nr. 9 A 1.00018703160 -0.0030285251761 0.0420515821101
B 1.00653488091 0.0123038027924 -0.0137546617228
Nr. 10 A 1.00016292645 -0.0029686792139 0.0425428553645
B 1.00646443609 0.0123164554124 -0.0136677101959
Nr. 11 A 1.00001497688 -0.0024895543937 0.0469651397961
B 1.00588177665 0.0123818197139 -0.0128808614831
Nr. 12 A 0.999994411581 -0.0023958169776 0.0479605410589
B 1.00576294309 0.0123879012149 -0.0127065079406
Nr. 13 A 0.999950688084 -0.0021420140511 0.0509494338767
B 1.00543105500 0.0123945720219 -0.0121964602104
Nr. 14 A 0.999906016810 -0.0014363159374 0.0629279255105
B 1.00440232257 0.0123446184006 -0.0104150756731
Nr. 15 A 0.999913431613 -0.0011421976120 0.0709224693220
B 1.00390589312 0.0122954359845 -0.0094544451251
Nr. 16 A 0.999996742223 -0.0000144169691 0.6490999947120
B 1.00039321477 0.0114924646536 -0.0011450792680
Nr. 17 A 0.999996913628 -0.0000051561989 1.0857415544000
B 1.00022058771 0.0111589135151 -0.0006877016857


edit...it looks good in this window, but oh well...

By the way...these are non-dimensional coordinates in the rotating Sun-Earth frame

 

[remcook]

didn't read your post well enough, but I think that book also has stable hamiltonians in them.

 

[Heimdall]

hum, thanks for all your answers.

I found a website where someone does what I'm trying to do.. unfortunatly for you it is in french, I'll try to translate and explain the problem.

It talks about horses orbits, which is a trajectory I can obtain in my model of restricted 3 bodies problem. well, with the hamiltonian :

H=\frac{1}{2} (P_1^2+P_2^2) +P_1 Q_2 - P_2Q_1 + (\frac{1-\mu}{R_1} + \frac{\mu}{R_2})

with :
R_1 = \sqrt{(Q_1+\mu)^2+Q_2^2}
and R_2 = \sqrt{(Q_1+\mu-1)^2+Q_2^2}

The guy wants to study such closed trajectories.

remebering there is the jacobi's integral J=-2K, he tells that he makes it equal to 3.0015051185.

and in all his study, P2 = 0 and Q2 = 0.002.

then he'll change P1 et Q1.

in his first example (you can see it

http://www.oamp.fr/lasco/pages_perso/danjard/3corps/P1=0.001.html
P1 = 0.001 and Q1 = 0.02.

well if I calcul J whit theses values I find arount 94, and not around 3 !

so I don't understand... how can he tell that J equals 3 and change the value of Q1 and P1 ?

therfore, I tried to put my particule at the L4 point, with a \mu = 0.01 and without relative velocity, where it was supposed not to move... and it moved away very fast... :/

thanks if you can help me :-D

 

[pervect]

hum, thanks for all your answers.

I found a website where someone does what I'm trying to do.. unfortunatly for you it is in french, I'll try to translate and explain the problem.

It talks about horses orbits, which is a trajectory I can obtain in my model of restricted 3 bodies problem. well, with the hamiltonian :

H=\frac{1}{2} (P_1^2+P_2^2) +P_1 Q_2 - P_2Q_1 + (\frac{1-\mu}{R_1} + \frac{\mu}{R_2})

I ran this into maple and stared at it quite a while until I finally saw your sign error. H should be

H=\frac{1}{2} (P_1^2+P_2^2) +P_1 Q_2 - P_2Q_1 - (\frac{1-\mu}{R_1} + \frac{\mu}{R_2})

because the gravitational potential energy is negative, not positive.

So with H=
<br /> 1/2\,{{\it p1}}^{2}+1/2\,{{\it p2}}^{2}+{\it p1}\,{\it q2}-{\it p2}\,{<br /> \it q1}-{\frac {1-u}{\sqrt {{{\it q1}}^{2}+2\,{\it q1}\,u+{u}^{2}+{{<br /> \it q2}}^{2}}}}-{\frac {u}{\sqrt {{{\it q1}}^{2}+2\,{\it q1}\,u-2\,{<br /> \it q1}+{u}^{2}-2\,u+1+{{\it q2}}^{2}}}}<br />

cut directly from the workspace with the latex output function to eliminate any possible errors, though the computer does some silly things) we now have a stable Lagrange point at

q1 := .5-u
q2 := sqrt(3)/2
p1 = -q2
p2 = q1

and \frac{\partial H}{\partial x} is zero for x=q1,q2,p1,p2

I don't know if this fixes all your problems, but it should be a start.

 

[Heimdall]

I ran this into maple and stared at it quite a while until I finally saw your sign error. H should be

H=\frac{1}{2} (P_1^2+P_2^2) +P_1 Q_2 - P_2Q_1 - (\frac{1-\mu}{R_1} + \frac{\mu}{R_2})

because the gravitational potential energy is negative, not positive.

oh thank you to check... but it seems unfortunatly that it was just a mistake that I made here... since my calculs seems to be correct.

I have the following system for the equations of movement :

\dot{Q_1} = P_1 + Q_1
\dot{Q_2} = P_2 - Q_1
\dot{P_1} =P_2 - (1-\mu) \frac{Q_1 + \mu}{R_1^3} - \mu \frac{Q_1+\mu-1}{R_2^3}
\dot{P_2} =-P_1 - (1-\mu) \frac{Q_2}{R_1^3} - \mu \frac{Q_2}{R_2^3}

q1 := .5-u
q2 := sqrt(3)/2
p1 = -q2
p2 = q1

yes I have the same results

and \frac{\partial H}{\partial x} is zero for x=q1,q2,p1,p2

yes... since all the \frac{\partial q_i,p_i}{\partial t}=0

I don't know if this fixes all your problems, but it should be a start.

well unfortunatly no... :-)
I checked my equations but I saw 0 mistakes, the same for my C functions. Problem must but elsewhere.

well in fact your post made me think of my initial conditions.
I was writting for m=0.001

Q1 = 0.499
Q2 = 0.8660254 /*sqrt(3)/2*/
P1 = 0
P2 = 0

0 for P1,2 because I thought there was no velocity at the lagrange points, thanks to you I remember know that the condition of the lagrange points is P1 = -Q2 and P2 = Q1

so know I write :
Q1 = 0.499
Q2 = 0.8660254 /*sqrt(3)/2*/
P1 = -0.8660254
P2 = 0.499

look what I have :http://nicolas.aunai.free.fr/cours/magistere/p3cr/test.jpg


and have you got an idea about my question on the jacobi's integral ?

 

[pervect]

oh thank you to check... but it seems unfortunatly that it was just a mistake that I made here... since my calculs seems to be correct.

I have the following system for the equations of movement :

\dot{Q_1} = P_1 + Q_1
\dot{Q_2} = P_2 - Q_1
\dot{P_1} =P_2 - (1-\mu) \frac{Q_1 + \mu}{R_1^3} - \mu \frac{Q_1+\mu-1}{R_2^3}
\dot{P_2} =-P_1 - (1-\mu) \frac{Q_2}{R_1^3} - \mu \frac{Q_2}{R_2^3}

More typos?

I get
q1' = p1 + q2


(for ease of typing I'm writing \frac{d q_1}{d t} as q1', and similarly for other variables - a prime indicates differeintiation wrt time).

The conditions at the Lagrange point are that q and p are constant, this is why q1', q2', p1' and p2' must all be zero at the lagrange point. You should check that they actually are zero in your code when you evaluate your functions at the Lagrange point. If they are not zero, you have a typo in your code.

0 for P1,2 because I thought there was no velocity at the lagrange points, thanks to you I remember know that the condition of the lagrange points is P1 = -Q2 and P2 = Q1

Note that these initial conditions follows simply from the fact that q1' and q2' are zero.


Representing \frac{\partial H}{\partial p1} as diff(H,p1), I get the following expressions that may be useful as a double-check, though I didn't see any difference from what you wrote except for the one error I already pointed out

diff(H1,p1);
p1 + q2

diff(H1,p2);
p2 - q1

diff(H,q1);
-p2 + (1 - u) (q1 + u)/R1^3

diff(H,q2);
p1 + (1-u) q2 / R1^3 + u q2 / R2^3


As I said earlier, iif your code is predecting motion at the Lagrange point, you must have the equations wrong, because we've just shown analytically that at the Lagrange point, p1' = p2' = q1' = q2' = 0.

and have you got an idea about my question on the jacobi's integral ?

I'm not sure I understood what your exact question was with respect to the Jacobi - could you repeat it? Since J = -2H, and we have the expression for the Hamiltonian written down, it should just be a matter of plugging the values of p and q into the equation to get J. I think you probably need to set m=1, see for instance

http://scienceworld.wolfram.com/physics/JacobiIntegral.html

their expression for J matches our expression for H only when m=1.

[add]

Actually, the equations we've written down so far are only valid if m=1. If you are assuming the test mass has a mass other than 1, that might be the soruce of more than one problem

 

[Heimdall]

More typos?

I get
q1' = p1 + q2

AARRRGH


no, not a typo this time. I mean not really. I explain you... I first made the calculus, then when it was all good theoritically I typed it and printed it. What I didn't know at this moment, is that I had this typo q1' = p1+q1 instead of p1' = p1+q2 !

I made x times the calculus, and found x time q2, but it seems my eyes just saw what they wanted to see...

thank you sometimes it is good to have an outside point of view of a stupid problem like this.

(for ease of typing I'm writing \frac{d q_1}{d t} as q1', and similarly for other variables - a prime indicates differeintiation wrt time).

The conditions at the Lagrange point are that q and p are constant, this is why q1', q2', p1' and p2' must all be zero at the lagrange point.

yes I know, this is how I found the lagrange points too.

You should check that they actually are zero in your code when you evaluate your functions at the Lagrange point. If they are not zero, you have a typo in your code.

well that ain't important really here, because the x's values are just initialized from the equation of movement, no matter of what they have in them first. Then they're just used in some algorithms like fourth order runge kutta for exemple.

Note that these initial conditions follows simply from the fact that q1' and q2' are zero.

yep :)

Representing \frac{\partial H}{\partial p1} as diff(H,p1), I get the following expressions that may be useful as a double-check, though I didn't see any difference from what you wrote except for the one error I already pointed out

yep thank you :)

As I said earlier, iif your code is predecting motion at the Lagrange point, you must have the equations wrong, because we've just shown analytically that at the Lagrange point, p1' = p2' = q1' = q2' = 0.

hum, I'm not sure to understand very well what you're saying... (sorry, I'm french, it doesn't help :D)
do you say that to study the movement of a body at a stable lagrange point I have to use another set of equations ?

or are you saying that I have to vary a little bit the initial condition to put my particule close to the lagrange poitn but not exactly there... ? if it is what you mean, I thought that my approximate values fo Q_i P_i for my initial condition would make big enough error to observe a mouvement.

I'm not sure I understood what your exact question was with respect to the Jacobi - could you repeat it? Since J = -2H, and we have the expression for the Hamiltonian written down, it should just be a matter of plugging the values of p and q into the equation to get J. I think you probably need to set m=1, see for instance

http://scienceworld.wolfram.com/physics/JacobiIntegral.html

their expression for J matches our expression for H only when m=1.

hum, yeah the fact that J=-2H, I agree. for exemple with :
Q1 = 0.499
Q2 = 0.8660254
P1 = -0.8660254
P2 = 0.499

so with mu = 0.001, I have J = 2.999001.

now look at this page http://www.oamp.fr/lasco/pages_perso/danjard/3corps/numeric3.html

the guy tells :
"Nous garderons dans tous les cas q2=0,002 et p2=0. La masse réduite sera de 0,001 et l'intégrale de Jacobi J=3,0015051185."

translation : we'll kepp in all cases q2=0.002 and p2=0. ratio mass will be 0.001 et jacobi's integral J = 3,0015051185.


well, he gives us q2 and p2, and J ! (J so to have closed trajectories), and then he gives different values for q1 and p1 ! since J value is determined by the combinaition of the q_i,p_i variables, how can he change p1 and q1 freely keeping the J he wants ??

and therefore, this is what he does : http://www.oamp.fr/lasco/pages_perso/danjard/3corps/P1=0.001.html

 

[Heimdall]

ok I misunderstood... 0.001 etc... were the distance form the lagrange point !

it is better since I know this : http://nicolas.aunai.free.fr/cours/magistere/p3cr/test.jpg

 

[pervect]

ok I misunderstood... 0.001 etc... were the distance form the lagrange point !

it is better since I know this : http://nicolas.aunai.free.fr/cours/magistere/p3cr/test.jpg

OK, I think I understand my confusion now - I thought the .001 was a mass, not a distance.

It sounds like you have everything figured out now - good luck with your project.

 

[Heimdall]

thank you very much for your help !