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 o ) + f ’(x o) (x – x o) + f ” (x o) (x – x o)2 / 2! +……
Applied after reconfiguring the Kepler equation to show:
( E o - e sin E o ) = M o
Representing a first iteration for starting values of E o and M o.
f(E) = f(E o ) + f ’(E o) (E
– E o) + f ” (E o) (E – E o)2 /
2! +……
Or, substituting relevant values:
f(E) = ( E o - e sin E o ) + (1 – e cos E o) (E – E o) + (E – E o)2 (e sin E o )/ 2
As a first approximation, let E o = M o, then:
f(E 1) = ( M o - e sin M o ) + (1 – e cos M o) (E – M o) + (E – M o)2 (e sin ME0)/ 2
And: E 1 is computed iteratively so that:
| E 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 = E o
Then convert to radians: M' (rad) = p(24o .35)/ 180o = 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' = E o = 0.425 and the goal is to get the differences in E from iteration down to 0.001 rad or less.
Write for E 1
E 1= M' + e sin E o = 0.425 + 0.5 sin (0.425)
E 1 = 0.631
And: E 1 - E 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:
E 2= M' + e sin E 1 = 0.425 + 0.5 sin (0.631)
E 2= 0.720
And: E 2 - E 1 = 0.089
Still not low enough so go to the next iteration:
E 3 = M' + e sin E 2 = 0.425 + 0.5 sin (0.720)
E 3 = 0.755
And: E 3 - E 2 = 0.035
Still not low enough so go to the next iteration:
E 4 = M' + e sin E 3 = 0.425 + 0.5 sin (0.755)
E 4 = 0.768
E 4 - E 3 = 0.013
Still too large so another:
E 5 = M' + e sin E 4 = 0.425 + 0.5 sin (0.768)
E 5 = 0.772
E 5 - E 4 = 0.004
Within 'striking' range but not there yet:
E 6 = M' + e sin E 5 = 0.425 + 0.5 sin (0.772)
E 6 = 0.774
E 6 - E 5 = 0.002
Once more:
E 7 = M' + e sin E 6 = 0.425 + 0.5 sin (0.774)
E 7 - E 6 = 5.968 x 10-4
Which meets the threshold.
Suggested Problems:
1) The mean anomaly (M) of Mercury is 179 o .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 o .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:
Post a Comment