Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Initial condition of the restricted problem of three bodies

  1. Feb 26, 2005 #1
    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 lets 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 [Broken])

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

    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
    Last edited by a moderator: May 1, 2017
  2. jcsd
  3. Feb 26, 2005 #2


    User Avatar
    Staff Emeritus
    Science Advisor

    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


    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 equilbirium 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

    \dot{p} = -\frac{\partial H}{\partial q}
    \dot{q} = \frac{\partial H}{\partial p}

    For this particular problem, we are particularly interested in the locaiton of L1. As an equilbirium 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 equilbirium 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

    Last edited by a moderator: May 1, 2017
  4. Feb 27, 2005 #3
    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
  5. Feb 27, 2005 #4
    didn't read your post well enough, but I think that book also has stable hamiltonians in them.
  6. Feb 28, 2005 #5
    hum, thx for all your answers.

    I found a web site 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 :

    [tex]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})[/tex]

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

    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 [Broken]
    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 [tex]\mu = 0.01[/tex] and without relative velocity, where it was supposed not to move... and it moved away very fast... :/

    thanks if you can help me :-D
    Last edited by a moderator: May 1, 2017
  7. Feb 28, 2005 #6


    User Avatar
    Staff Emeritus
    Science Advisor

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

    [tex]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})[/tex]

    because the gravitational potential energy is negative, not positive.

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

    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 [tex]\frac{\partial H}{\partial x}[/tex] is zero for x=q1,q2,p1,p2

    I don't know if this fixes all your problems, but it should be a start.
  8. Mar 1, 2005 #7
    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 :

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

    yes I have the same results

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

    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 [Broken]

    and have you got an idea about my question on the jacobi's integral ?
    Last edited by a moderator: May 1, 2017
  9. Mar 1, 2005 #8


    User Avatar
    Staff Emeritus
    Science Advisor

    More typos?

    I get
    q1' = p1 + q2

    (for ease of typing I'm writing [tex]\frac{d q_1}{d t} [/tex] 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.

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

    Representing [tex]\frac{\partial H}{\partial p1}[/tex] 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

    p1 + q2

    p2 - q1

    -p2 + (1 - u) (q1 + u)/R1^3

    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.

    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


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


    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
    Last edited: Mar 1, 2005
  10. Mar 1, 2005 #9

    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.

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

    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.

    yep :)

    yep thank you :)

    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.

    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 [Broken]

    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 [Broken]
    Last edited by a moderator: May 1, 2017
  11. Mar 1, 2005 #10
    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 [Broken]
    Last edited by a moderator: May 1, 2017
  12. Mar 2, 2005 #11


    User Avatar
    Staff Emeritus
    Science Advisor

    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.
    Last edited by a moderator: May 1, 2017
  13. Mar 2, 2005 #12
    thank you very much for your help !!
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook