Planets and solution of n-degree polynomials

Post Reply
Henko
Posts: 814
Joined: Tue Apr 09, 2013 12:23 pm
My devices: iPhone,iPad
Windows
Location: Groningen, Netherlands
Flag: Netherlands

Planets and solution of n-degree polynomials

Post by Henko »

"Lagrange points" around two celestial bodies are positions where a satelite can remain in the same position in relation to both bodies, given a certain speed.
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
L1 = 1.49141E+06 km

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
(2x^2-3x+5)(x-1)(x+3)(3x-1)(x-2) = 0
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

Henko
Posts: 814
Joined: Tue Apr 09, 2013 12:23 pm
My devices: iPhone,iPad
Windows
Location: Groningen, Netherlands
Flag: Netherlands

Re: Planets and solution of n-degree polynomials

Post by Henko »

A textual incorrectness:
The 5000 km difference is not 0.3 % of the distance between the sun and the earth, but 0.3 % of the calculated L1 value. Hence it is the relative error in the L1 calculation when using the approximate formula. L1 itself is about 1 % of the distance between the sun and the earth.

User avatar
GeorgeMcGinn
Posts: 435
Joined: Sat Sep 10, 2016 6:37 am
My devices: IPad Pro 10.5in
IMac
Linux i386
Windows 7 & 10
Location: Venice, FL
Flag: United States of America
Contact:

Re: Planets and solution of n-degree polynomials

Post by GeorgeMcGinn »

For me, this is right in my wheel house. Henko explains it well, well enough to get those interested in what and where these points are.

Also for those more advanced, you could create Lagrange points between any two bodies with mass that orbit. That is beyond the scope of these posts.

As Henko said, most satellites put in an L1 to L5 point stay in place. However, with the James Webb telescope, it will orbit around the Sun at L2, which is about 1.5 million miles further from the Earth, which is a more stable orbit and will follow the Earth as it orbits the Sun.

If you are going to develop a game where objects (satellites and even asteroids) orbit in the various L1-L5 points, here is some low-level technical data from NASA, including To graphics: one showing where the Lagrange points are located, and the other showing the application of putting the James Webb Space Telescope in L2. It is important to know this. If your game says it's putting a solar observatory in L2, many will know you put these in L1 instead. This will give you just an overview on accuracy.

Some Technical Details: It is easy for an object (like a spacecraft) at one of these five points to stay in place relative to the other two bodies (e.g., the Sun and the Earth). In fact, L4 and L5 are stable in that objects there will orbit L4 and L5 with no assistance. Some small asteroids are known to be orbiting the Sun-Earth L4 and L5 points. However, L1, L2, and L3 are metastable so objects around these points slowly drift away into their own orbits around the Sun unless they maintain their positions, for example by using small periodic rocket thrust. This is why L1, L2, and L3 don't "collect" objects like L4 and L5 do.
The James Webb Space Telescope will not be in orbit around the Earth, like the Hubble Space Telescope is - it will actually orbit the Sun, 1.5 million kilometers (1 million miles) away from the Earth at what is called the second Lagrange point or L2. What is special about this orbit is that it lets the telescope to stay in line with the Earth as it moves around the Sun. This allows the satellite's large sunshield to protect the telescope from the light and heat of the Sun and Earth (and Moon).
What Is L2?
Joseph-Louis Lagrange was an 18th century mathematician who found the solution to what is called the “three-body problem.” That is, is there any stable configuration, in which three bodies could orbit each other, yet stay in the same position relative to each other? As it turns out, there are five solutions to this problem - and they are called the five Lagrange points, after their discoverer. The L1, L2, and L3 points are all in line with each other - and L4 and L5 are at the points of equilateral triangles.
Attachments
IMG_0394.JPG
IMG_0394.JPG (23.93 KiB) Viewed 1955 times
George McGinn
Computer Scientist/Cosmologist/Writer/Photographer
Member: IEEE, IEEE Computer Society
IEEE Sensors Council & IoT Technical Community
American Association for the Advancement of Science (AAAS)

Post Reply