There are 5 such points. One of them lies between the two bodies on the connection axis between the bodies. It's called "L1".
In Wikipedia a formula is derived to calculate the position of L1 in the Sun - Earth system (the L1 distance is about 1.5 million km from the earth). Theoretically, a 5th degree polynomial equation must be solved. To overcome this, an approximation is made in the article, reducing the result to a simple 3rd degree root in the formula.
However, L1 is an equilibrium position, but it is an unstable one: a small deviation from it will grow worse. Hence it is interesting to know what size of error is introduced by the simplification in the Wikipedia calculation.
That's why I decided to make a little tool tool to find the zero's (roots) of any n-degree polynomial. It was peanut's, because i already had a "polynomial" functions library.
The solution method is based upon the Newton-Raphson iteration scheme (my favorite). The function does only find real roots, no complex roots.
There are two testcases
1. Calculation of the exact value of L1 (differs about 5000 km from the approximate value by the Wikipedia formula; seems a lot, but is only 0.3 % of the distance between the Sun and the Earth)
2. A 6th degree polynomial, constructed via known roots, from which one pair of complex roots (and hence 4 real roots).
Both cases are presented hereafter, with the output.
I realize, that this is not "hot stuff" for a lot of us. (not-hot-lot , a game idea is coming up)
Code: Select all
' test for poly_roots() function
' L1 Lagrange point for Sun - Earth system is used
n=5 ! n1=n+1
dim eq(n1),roots(n1),diff(n1), xx(n1)
' basic data
M_sun=1.989E+30 ' kg ( mass of Sun )
M_earth=5.972E+24 ' kg ( mass of Earth )
alfa = M_sun / M_earth
d = 1.496E+11 ' m ( distance Sun <-> Earth )
eq(0) = -(d^5)
eq(1) = +2*d^4
eq(2) = -(d^3)
eq(3) = 3*alfa*d^2
eq(4) = -3*alfa*d
eq(5) = alfa
poly_roots(n,eq,roots,d,.001,50)
print "L1 = ";roots(1)/1000;" km"
end
' finding the (real) roots of a n-degree polynomial
' option base 0 is assumed
' n is the degree of the polynomial
' Poly() contains the n+1 coefficients p(i) in the following order:
' Poly(n) = p(0)+p(1).x+.....+p(n).x^n
' roots() contains the roots, returned by the function, starting with
' roots(1). If there are less than n real roots, the series of
' roots is closed with -999
' xo is a starting value for the proces, is not critical generally
' precision is the required absolute +/- precision
' max_iterations is the allowed max # of iterations
'
def poly_roots(n,Poly(),roots(),xo,precision,max_iterations)
dim eq(n+1),out(n+1)
for i=0 to n ! eq(i)=Poly(i) ! next i
for i=n to 1 step -1
x=poly_root(i,eq,xo,precision,max_iterations)
k=n-i+1 ! roots(k)=x ! if x=-999 then return x
poly_divide(i,eq,out,x)
for j=0 to i-1 ! eq(j)=out(j+1) ! next j ! eq(i)=0
next i
return 1
end def
' find one (real) root of p(n)
'
def poly_root(n,cin(),xo,eps,max_iter)
dim diff(n)
poly_diff(n,cin,diff) ! x=xo
for iter=1 to max_iter
dx=poly_value(n,cin,x)/poly_value(n-1,diff,x) ! x-=dx
if abs(dx)<eps then return x
next iter
return -999
end def
def poly_diff(n,cin(),cout())
for i=0 to n-1 ! cout(i)=(i+1)*cin(i+1) ! next i
end def
def poly_value(n,cin(),x)
val=cin(n)
for i=n-1 to 0 step -1 ! val*=x ! val+=cin(i) ! next i
return val
end def
def poly_divide(n,cin(),cout(),a)
for i=0 to n ! cout(i)=cin(i) ! next i
for i=n-1 to 0 step -1 ! cout(i)+=a*cout(i+1) ! next i
return cout(0)
end def
Code: Select all
' test for poly_roots() function
' 6th degree polynomial with 4 real roots
n=6 ! n1=n+1
dim eq(n1),roots(n1),diff(n1), xx(n1)
print " (2x^2-3x+5)(x-1)(x+3)(3x-1)(x-2) = 0"
print "equals : 6x^6-11x^5-24x^4+108x^3-192x^2+143x-30 = 0" ! print
for i=0 to n ! read eq(i) ! next i
data -30, 143, -192, 108, -24, -11, 6
poly_roots(n,eq,roots,d,.001,50)
for i=1 to n ! print "x(";i;") = ";roots(i) ! next i
end
' finding the (real) roots of a n-degree polynomial
' option base 0 is assumed
' n is the degree of the polynomial
' Poly() contains the n+1 coefficients p(i) in the following order:
' Poly(n) = p(0)+p(1).x+.....+p(n).x^n
' roots() contains the roots, returned by the function, starting with
' roots(1). If there are less than n real roots, the series of
' roots is closed with -999
' xo is a starting value for the proces, is not critical generally
' precision is the required absolute +/- precision
' max_iterations is the allowed max # of iterations
'
def poly_roots(n,Poly(),roots(),xo,precision,max_iterations)
dim eq(n+1),out(n+1)
for i=0 to n ! eq(i)=Poly(i) ! next i
for i=n to 1 step -1
x=poly_root(i,eq,xo,precision,max_iterations)
k=n-i+1 ! roots(k)=x ! if x=-999 then return x
poly_divide(i,eq,out,x)
for j=0 to i-1 ! eq(j)=out(j+1) ! next j ! eq(i)=0
next i
return 1
end def
' find one (real) root of p(n)
'
def poly_root(n,cin(),xo,eps,max_iter)
dim diff(n)
poly_diff(n,cin,diff) ! x=xo
for iter=1 to max_iter
dx=poly_value(n,cin,x)/poly_value(n-1,diff,x) ! x-=dx
if abs(dx)<eps then return x
next iter
return -999
end def
def poly_diff(n,cin(),cout())
for i=0 to n-1 ! cout(i)=(i+1)*cin(i+1) ! next i
end def
def poly_value(n,cin(),x)
val=cin(n)
for i=n-1 to 0 step -1 ! val*=x ! val+=cin(i) ! next i
return val
end def
def poly_divide(n,cin(),cout(),a)
for i=0 to n ! cout(i)=cin(i) ! next i
for i=n-1 to 0 step -1 ! cout(i)+=a*cout(i+1) ! next i
return cout(0)
end def
equals : 6x^6-11x^5-24x^4+108x^3-192x^2+143x-30 = 0
x( 1 ) = 0.333333
x( 2 ) = 1
x( 3 ) = 2
x( 4 ) = -3
x( 5 ) = -999
x( 6 ) = 0
