Stability of collinear Lagrange points

Question

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?

Answer 1

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.012277471`20;

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.04`20}, {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]

Answer 2

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.