A little bit of Mathematics.

[GWS]

Hi folks,

It crossed my mind that some visitors to this site, might be wondering if CBasic could be useful for mathematical studies.

Well of course it can  .. and here is a little example relating to Kepler's Equation of Elliptical Motion.

A good description of this problem can be found here:

http://www.akiti.ca/KeplerEquation.html

- and there is much more discussion on the Web.

The problem is a general one in mathematics - to find the root of a nonlinear equation.

If f(x) = ax^2 + b*x + c = 0, we seek the value(s) of x which satisfy this equation - these are the 'root' values.

One method of solution is Linear Iteration.  This can have it's little problems, but for this example these will not arise.

Given f(x) = 0, we first express this as x = g(x).
There are often several ways of writing this, and not all of them are suitable.

For example with an equation:   f(x) = x^2 - x - 2 = 0
we can write:

(a) x = x^2 - 2
(b) x = sqrt(2 + x)
(c) x = 1 + 2/x

What we are looking for, is a form where we can generate successive approximations to the required root by the recursion formula:
  x _i+1  = g(x _i).

(This means we get the next approximation to x from the expression g(x) using the current value of x).

In the example above, option (a) will not work for this method, but option (b) will ..

Using x = 0 as a starting value we get:
x 1 = sqrt(2) = 1.414
x 2 = sqrt(3.414) = 1.847
x 3 = sqrt(3.847) = 1.961
x 4 = sqrt(3.961) = 1.990
x 5 = sqrt(3.990) = 1.9975

... and we are obviously converging on the result 2.0 , which is the correct answer.

So with that preliminary discussion, let's have a look at Kepler's equation of Elliptical Motion.

f(E) = M - E + e * sin(E) = 0

where E = the angle of Eccentric Anomaly (radians)
M = the angle of Mean Anomaly (radians)
e = the orbit eccentricity

We require to find the value of E - the Eccentric Anomaly, to at leat 6 figure accuracy.
To date, no explicit solution to this equation has been found - so we will attempt a Linear Iteration method.

To proceed, we re-write the equation to get g(E) = M + e * sin(E), and with a starting estimate for E = M.
(Note we calculate in radians since CBasic trigonometry functions work in radians).

The program uses a Console window, since this is the most direct approach for an infrequently used program.
If we intended to do lots of these orbit calculations it would be nicer to use a full Window user interface, with edit boxes for input values, and a result area.

To test the program, try using 28.0 degrees as the Mean Anomaly, 0.2 as the Eccentricity, and 1.0e-10 as the acceptable Convergence difference limit.
You should get the result 38.48866974 degrees in 13 iterations.

Here's the program ..

Code Select Expand



' Example of Linear Iteration to solve Kepler's Eccentric Anomaly equation
' GWS Jan 2009
openconsole
cls

def MeanAnomaly,Eccentricity,EccentricAnomaly,difference,Delta,x,pi:double
def a$,b$,c$:string

input "Input Mean Anomaly (Degrees):  ",MeanAnomaly
input "Input Eccentricity:  ",Eccentricity
' express this in the form (eg. 0.00001 or 1.0e-5) whichever is convenient.
input "Iteration Convergence Difference:  ",Delta

setprecision 8 :' for output purposes

pi = 4 * atan(1)
d2r = pi/180 :' degrees to radians conversion factor

' convert MeanAnomaly to radians ..
MeanAnomaly = MeanAnomaly * d2r

' Kepler's equation is:  f(E) = MeanAnomaly - E + Eccentricity * sin(E) = 0
' from which:  E = M + e * sin(E)
' (where E is the Eccentric Anomaly)

cls

color 10,0
print "Solution of Kepler's Eccentric Anomaly":print
color 7,0

' try up to 20 iterations to converge ..
for i = 1 to 20

if(i = 1)
' start with the first estimate eqaul to the Mean Anomaly ..
x = MeanAnomaly
else
x = EccentricAnomaly
endif

EccentricAnomaly = MeanAnomaly + Eccentricity * sin(x)
difference = abs(EccentricAnomaly - x)

a$ = using("###",i)
b$ = using("   #.########",EccentricAnomaly)
c$ = using("   #.##########",difference)


print "Iteration " + a$ + "    Eccentric" + b$ + "    Difference" + c$
if (difference < Delta)
print:print:color 12,0
print "Eccentric Anomaly (Degrees)  =  ", x/d2r
i = 20
endif

next i

do:until inkey$<>""
closeconsole
end



all the best,

Graham

[LarryMc]

Nice discussion and example Graham.

Larry