Most treatments of Lagrange points assume two orbiting masses and a test point at one of the five Lagrange points. One of the masses is considered infinitely smaller than the other and then the test point infinitely smaller than either. There are thus two limiting processes in play here and this leads to pit-falls that confuse me at least, and perhaps accounts for some of the discrepancies about the Lagrange points that are found on the web. The plan here is to deal only with the limit as the mass of the test particle goes to zero. The two larger masses are comparable, or even equal.

Let r1 and r2 be the positive radii of the two masses, M1 and M2, about their common center of gravity. Equating the gravitational force to the centrifugal forces generally gives

M1M2G/(r1+r2)2 = ω2r1M1 = ω2r2M2

To shorten the equations I will assume that the angular velocity ω = 1 and that the sum of the two masses is 1—indeed that the masses are m and (1− m). I also assume units that make the gravitational constant G = 1. In our case this becomes

m(1−m)/(r1+r2)2 = r1m = r2(1−m)

r1 and r2 are measured from the center of gravity and their ratio is the inverse of the mass ratio:

r1/r2 = (1−m)/m,
or
r1 = r2(1−m)/m

Substituting this r1 into the first equation gives

m(1−m)/(r2(1−m)/m+r2)2 = r2(1−m)

and multiplying by r22/(1−m) gives

m/((1−m)/m+1)2 = r23
or
m3 = r23
or
r2 = m

and by symmetry

r1 = 1−m

I find this suspiciously simple! We will take the position of the first particle to be negative, for they are on opposite sides of the origin. r1 is the distance from the origin so the r coordinate of the first particle is m − 1. Now the acceleration of a test particle at radius r on the line thru the two masses:

a = r + m/(r − (m − 1))2 + (1 − m)/(r − m)2

Alas the plus signs in the above equation must be changed to minus depending on which Lagrange point we seek as the masses always attract (push towards the mass, whichever direction) whereas the terms are always positive and push towards positive r. The above form is for L2, nearer the secondary than the primary. As the expression stands this leads to a polynomial in r which must be 0 for equilibrium: This scheme code finds the following:
p = ((1)r5 + (−4m + 2)r4 + (6m2 + −6m + 1)r3 + (−4m3 + 6m2 + −2m + 1)r2 + (1m4 + −2m3 + 1m2 + −4m + 2)r + (3m2 + −3m + 1)).
First I try m = 0.1 to calculate a real number for L1:

r5 + 1.6r4 + .46r3 + .856r2 + 1.6081r + .73 = 0

With the kind help of Polynomial root finder, Newton’s method and Mathematica we find a root at − 1.2596998329023314374. That solution makes the expression for the acceleration vanish and seems physically reasonable.

Now we move to the interval between the two masses to seek L1 and flip the sign of the attraction of the first mass getting:

a = r − m/(r − (m − 1))2 + (1 − m)/(r − m)2

from which our program derives:

((1)r5 + (−4m + 2)r4 + (6m2 + −6m + 1)r3 + (−4m3 + 6m2 + −4m + 1)r2 + (1m4 + −2m3 + 5m2 + −4m + 2)r + (−2m3 + 3m2 + −3m + 1))

And in our case for m = 0.1: (1r5 + 8/5r4 + 23/50r3 + 82/125r2 + 16481/10000r + 91/125)
which has a root at −0.60903511002320256380.

And finally:

a = r − m/(r − (m − 1))2 − (1 − m)/(r − m)2

for L3 yielding polynomial

((1)r5 + (−4m + 2)r4 + (6m2 + −6m + 1)r3 + (−4m3 + 6m2 + −2m + −1)r2 + (1m4 + −2m3 + 1m2 + 4m + −2)r + (−3m2 + 3m + −1))

which for m = 0.1 is

(1r5 + 8/5r4 + 23/50r3 − 143/125r2 − 15919/10000r − 73/100)

which has a root at 1.07164421940572131 .

For m = 0.000002, as for the Sun-Earth combination, we get solutions:

L2: −1.0087591722713581551
L1: −0.99128770335330933876
L3: 1.0000008333333333507

The units are AU’s. Keep in mind that these locations are relative not to the Sun but to the center of gravity. The distance from the Sun to L3 is 0.9999988333333335.

### The Whole Picture

This Java program produces this image: We call the three points at the top tl, tc and tr, and the three at the bottom, bl, bc and br. Each horizontal cut thru this picture locates the two masses and the three Lagrange points. The center of gravity is the always the center of the cut. The vertical coordinate is m with m=0 at the top and m=1 at the bottom. The horizontal coordinate is r and r=0 at tc and bc, the center of gravity. The first mass, M1 is on the straight line between tl and bc while M2 is on the line between tc and br. The Lagrange points are at the curved boundaries between black and white. L1 is between the two masses and L2 is closer to the smaller mass, and L3 is closer to the larger. The Sun is closer to L3 than to the Earth, but the center of gravity is closer to the Earth than to L3. The black areas are where a test particle would accelerate to the left.

### L4 and L5

The two points in the orbit’s plane that are each at distance 1 from each of the two masses are also at equilibrium. This includes the classic case found by Lagrange. If we decompose the acceleration at such a point into two orthogonal components, one parallel to the line between the masses, then the acceleration perpendicular to that line is constant due symmetry and the fact that sum of those two masses is constant. The other component is proportional to m−1/2 which is proportional to the difference between the masses. But these two components are just the same as the components of the vector displacement from the center of gravity to the equilibrium point.

The surprise is that for m=1/2 those points are unstable.
 distance from center gravitational force centrifugal force .86 .8736 .86 .87 .8610 .87
It would be unstable, that is, if the particle were constrained to slide along a line perpendicular to the line between the masses. Here is a fairly good resolution of the question.

Here is a promising applet with source. See too this compound pendulum together with source code.