Friday, October 25, 2024

Revisiting Numerical Analysis: Applications To The Kepler Equation Of Celestial Mechanics

 Given the importance of numerical analysis in solving a host of otherwise intractable problems, i.e. in plasma physics as well as  celestial mechanics, it is worth taking another look at it.  We can start with a relatively basic celestial mechanics application.

 The Problem in a nutshell: Given Kepler's equation:  

M =   E - e sin E 

Where M is the mean anomaly of an orbiting body, and E is the eccentric anomaly, with e the eccentricity of the orbit, how does one obtain E when M is known?  The simplest method is based on a formula:

tan E  = sin M/ (cos M - e)

But only works provided e   »  0 

Therefore more 'onerous' methods are needed for the more realistic case e > 0.

One possible approach is to use a Taylor series expansion,

An Introduction To Numerical Analysis - Part 1 (Preliminaries)

f(x) = f(x ) + f ’(x o) (x – x o) + f ” (x o) (x – x o)2 / 2!  +……

Applied after reconfiguring the Kepler equation to show:

( E  -  e sin E )  =  M  

Representing a first iteration for starting values of  and  o.

f(E) = f(E ) + f ’(E o) (E – E o) + f ” (E o) (E – E o)2 / 2!  +……

Or,  substituting relevant values:

f(E) = ( E  -  e sin E ) +  (1 – e cos  E o) (E – E o) +  (E – E o)2  (e sin E )/ 2

As a first approximation, let   o,  then:

f(1) = ( o -  e sin o ) +  (1 – e cos o) (E – o) + (E – o)2  (e sin ME0)/ 2

And:  E 1  is computed iteratively so that:  

| 1 – E  | <   | E  –  E o| 

But  The Newton-Raphson method is a more straightforward approach, e.g.

An Introduction To Numerical Analysis (3): Newton's Approximation Method

Which we will apply here. Begin by letting the eccentricity e = 0.5 and the mean anomaly: M' = 24o .35 =  o

 Then convert to radians:  M' (rad) =  p(24o .35)/ 180= 0.425

We use the basic Kepler eqn. to start:  M' = E - e sin E  

and rewrite it in the form:

E = M' +  e sin E

For starting the iterative process:   M' = 0.425    and the goal is to get the differences in E from iteration down to 0.001 rad or less.

 Write for 1

1=  M'  +  e sin o   = 0.425  +  0.5 sin (0.425)

= 0.631

And:  1 - o  = 0.206

The aim then is to keep the iterations going until the difference in E values if less than one thousandth of a radian.

So next we have:

2=  M'  +  e sin 1   = 0.425  +  0.5 sin (0.631)

2=  0.720

And: 2 - 1  = 0.089

Still  not low enough so go to the next iteration:

3 =  M'  +  e sin 2   = 0.425  +  0.5 sin (0.720)

=  0.755

And: 3 - 2  = 0.035

Still  not low enough so go to the next iteration:

4 =  M'  +  e sin 3   = 0.425  +  0.5 sin (0.755)

4 =  0.768

4 - 3  = 0.013

Still too large so another:

5 =  M'  +  e sin 4   = 0.425  +  0.5 sin (0.768)

5 =   0.772

5 - 4  = 0.004

Within 'striking' range but not there yet:

6 =  M'  +  e sin 5   = 0.425  +  0.5 sin (0.772)

6 = 0.774

6 - 5  = 0.002

Once more:

7 =  M'  +  e sin 6   = 0.425  +  0.5 sin (0.774)

7 - 6  = 5.968 x 10-4

Which meets the threshold.

Suggested Problems:

1) The mean anomaly  (M) of Mercury is 179 .796  and its eccentricity is e= 0.20563.

Use this to do a Newton-Raphson iteration to obtain the eccentric anomaly. 

2) The mean anomaly of Earth is 198 .115  and its eccentricity is e = 0.01672.

a) Use this data  to obtain the eccentric anomaly E using a Newton-Raphson iteration.

b) Compare the preceding result with E obtained from applying the Lagrange  expansion theorem, viz.

E = M + e sin M +  ( e 2/ 2) sin 2M +  e3/3! 2e 2  [3 2sin 3M - 3 sin 3M) + ...

c) How do the above results compare with E obtained from the simplest formula:

tan E  = sin M/ (cos M - e)

No comments: