# Stability of collinear Lagrange points

I am currently simulating the restricted 3-body problem in Mathematica. I have identified the Lagrange points and their coordinates. Now I want to see how particles move in their vicinity, to check that they are indeed stationary (for some period of time). Using the code below I get instability even if I place my particle exactly at the Lagrange points.

I have the equations of motion

$\ddot{x}-2\dot{y}=-\frac{\partial\Omega}{\partial x }$

$\ddot{y}+2\dot{x}=-\frac{\partial\Omega}{\partial y }$

$\Omega=-\frac{1}{2}\mu r_1^2-\frac{1}{2}(1-\mu)r_2^2-\frac{\mu}{r_1}-\frac{1-\mu}{r_2}$

Where $r_{1,2}$ are the distances between the small body and the two massive bodies , and $\mu$ is such that the ratio of the masses is $\mu:1-\mu$.

I have shown the L points to be at:

[{ {-0.609035110014235, 0}, {1.0416089085755191, 0}, {-1.2596998329049578, 0}, {-0.4000000000569969, 0.866025403864164}, {-0.4000000000569969, -0.866025403864164} }]


Then my code for investigating the particle motion is

r1[x_, y_, u_] := Sqrt[ (x + 1 - u)^2 + y^2]

u = 0.012277471

r2[x_, y_, u_] := Sqrt[(x - u)^2 + y^2]

Om[x_, y_, u_] := -1/2*u*r1[x, y, u]^2 - 1/2*(1 - u)*r2[x, y, u]^2 -
u/r1[x, y, u] - (1 - u)/r2[x, y, u]

sol = NDSolve[{x''[t] - 2 y'[t] == -D[Om[x[t], y[t], u], x[t]],
y''[t] + 2 x'[t] == -D[Om[x[t], y[t], u], y[t]],
x[0] == -0.4000000000569969, y[0] == 0.866025403864164, x'[0] == 0,
y'[0] == 0.09}, {x, y, x', y'}, {t, 30}, WorkingPrecision -> 75,
MaxSteps -> 100000, Method -> "StiffnessSwitching"]


I know that the differential equation is being solved numerically in a sufficiently accurate manner. But when I set my initial conditions such that the particle starts at a Lagrange point, it does not stay there. Why is this so?

I display the motion using this code:

traj = ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 30},
PlotRange -> All, ImageSize -> Large]
cpo = ContourPlot[{Om[x, y, u]}, {x, -2, 2}, {y, -2, 2},
PlotTheme -> "Scientific", PlotLegends -> Automatic,
Contours -> 100, ColorFunction -> "TemperatureMap"] ;
Show[cpo, traj]


EDIT:After reading @BillWatts' answer, instead of using the FindMaximum function I looked for L points by explicitly setting derivatives of $\Omega$ to zero. This gave different points, including (-0.8362925908999327,0), so I changed the initial conditions to

x[0] == -0.8362925908999327, y[0] == 0, x'[0] == 0, y'[0] == 0.0


and I still get particle motion. Why can this be?

• The problem I see is that the Lagrange points ("L points") you give near the beginning of your question do not correspond to the value of mu=0.012277471 you use in your code. The points you give correspond to mu=0.1. Hence, there is motion at these non-Lagrange points as @BillWatts shows. The answer by Alex Trounev calculates the correct Lagrange points "p" corresponding to the input mu=0.012277471. – KennyColnago Sep 11 '18 at 18:21

The Lagrange points and all other calculations are made with insufficient accuracy. The speed is set too high. After correction, the trajectory lies around the point of libration.

r1[x_, y_, u_] := Sqrt[(x + 1 - u)^2 + y^2]

u = 0.01227747120;

r2[x_, y_, u_] := Sqrt[(x - u)^2 + y^2]

Om[x_, y_, u_] := -1/2*u*r1[x, y, u]^2 - 1/2*(1 - u)*r2[x, y, u]^2 -
u/r1[x, y, u] - (1 - u)/r2[x, y, u]
Ux = D[-1/2*u*r1[x, y, u]^2 - 1/2*(1 - u)*r2[x, y, u]^2 -
u/r1[x, y, u] - (1 - u)/r2[x, y, u], x];
Uy = D[-1/2*u*r1[x, y, u]^2 - 1/2*(1 - u)*r2[x, y, u]^2 -
u/r1[x, y, u] - (1 - u)/r2[x, y, u], y];
p = {x, y} /. NSolve[Ux == 0 && Uy == 0, {x, y}];

sol = NDSolve[{x''[t] - 2 y'[t] == -D[Om[x[t], y[t], u], x[t]],
y''[t] + 2 x'[t] == -D[Om[x[t], y[t], u], y[t]],
x[0] == p[[1, 1]], y[0] == p[[1, 2]], x'[0] == 0,
y'[0] == 0.0420}, {x, y, x', y'}, {t, 30},
WorkingPrecision -> 19, MaxSteps -> 10^6];

traj = ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 30},
PlotRange -> All]
cpo = ContourPlot[{Om[x, y, u]}, {x, -2, 2}, {y, -2, 2},
PlotLegends -> Automatic, Contours -> 50,
ColorFunction -> "TemperatureMap"]; Lp =
Graphics[{White, PointSize[.01], Point[p]}];

Show[cpo, traj, Lp]


• To be clear, you increased accuracy by multiplying by 10 the number of steps, and that's all. Correct? What made you think it was Ok to reduce WorkingPrecision? Thanks. – tbfr416 Sep 11 '18 at 8:41
• Also, why did you not define Ux and Uy as D[Om[x,y,u],x] instead of re-typing the whole expression? It seems to indeed make a difference, which I don't understand. – tbfr416 Sep 11 '18 at 8:51
• I increased the accuracy of the representation of all input data to 20, for example, $u = 0.01227747120$, so the Lagrange points are calculated with an accuracy of 19 ..... The initial position and velocity are also specified with this accuracy, so all calculations are done with this accuracy WorkingPrecision -> 19. In your example, you used u = 0.012277471, but required WorkingPrecision -> 75. This is a typical error. You must put $u = 0.01227747175$, and calculate the Lagrange points with this accuracy. – Alex Trounev Sep 11 '18 at 9:11

Check your accelerations at the Lagrangian point and you will see that it is not zero.

Clear["Global*"]

u = 0.012277471

x1 = u - 1
x2 = u

r1[x_, y_, u_] := Sqrt[(x - x1)^2 + y^2]
r2[x_, y_, u_] := Sqrt[(x - x2)^2 + y^2]

Om[x_, y_, u_] := -1/2*u*r1[x, y, u]^2 - 1/2*(1 - u)*r2[x, y, u]^2 -
u/r1[x, y, u] - (1 - u)/r2[x, y, u]


Assuming you start stationary at your L pt. your accelerations are the derivatives of Om.

xaccel = -D[Om[x, 0.866025403864164, u], x] /. x -> -0.4
(*0.0551961*)

yaccel = -D[Om[-0.4, y, u], y] /. y -> 0.866025403864164
(*-0.112649*)
`

Clearly non-zero accelerations are causing the particle to move out. You might recheck your inputs.