Solving Kepler's equation using iterative technique

[dgalvan123]

I am working on using a Runge Kutta 4th order code (in the C programming language) to predict a Keplerian orbit in 2 dimensions. My code seems to be working ok, producing a reasonable elliptical orbit given initial conditions x0, y0, Vx0 and Vy0 (initial position and velocity), however I'd like to compare the code output to an "exact" solution to Kepler's equation, giving radius and true anomaly (r, theta) as a function of time.

The problem here is that, in order to get r and theta as a function of time, I need to solve for the eccentric anomaly (psi) at a given time. Kepler's equation is:

psi = w*t + e*Sin(psi).

Hence, I'm trying to invert and solve for psi given w (average angular velocity = 2pi/period), t (time) and e (eccentricity). Once I get psi as a function of t, I can use that to get position in 2 D.

My question is: how can I solve the Kepler equation? I know that some iterative technique will probably be needed. . . perhaps one where I let:

psi(0) = 0
and
psi(n+1) = w*t + e*Sin(psi(n))

(keeping time constant and just continuously advancing "n")

. . . but how will I know when to stop? How will I know when I've iterated enough such that psi(N) is a good approximation to the exact solution?

Thanks for your time!

 

[pkleinod]

Hi and welcome to PF!

An iterative approach should not be necessary.
I assume that the centre of force is at the focus of the ellipse on
the positive x-axis at a distance ae from the origin, where a is the
semi-major axis and e the eccentricity of the ellipse. Your numerical
solution then yields a position vector x(t) and a corresponding
velocity vector v(t) as functions of the time t. At a particular
time t, you know x(t) and therefore the vector r(t):
<br /> r(t) = x(t) - ae i,<br />
where i is the unit vector along the x-axis. So you know the
distance of the point x from the centre of force and you know the
angle theta that r makes with the x-axis. From the equation for
the (elliptical) orbit and the angle theta, compute the
exact magnitude of the vector r and compare with the value from
your code:

<br /> r_{ex} = \frac{a(1-e^2)}{1 + e \cos{\theta}}<br />

The values of the semi-major axis and eccentricity are obtained
from the total energy (E=K/r + L²/(2mr²))
and the angular momentum L, which are in turn specified by your
initial conditions.

 

[dgalvan123]

Thanks very much for your help!

 

[D H]

I am working on using a Runge Kutta 4th order code (in the C programming language) to predict a Keplerian orbit in 2 dimensions. My code seems to be working ok, producing a reasonable elliptical orbit given initial conditions x0, y0, Vx0 and Vy0 (initial position and velocity), however I'd like to compare the code output to an "exact" solution to Kepler's equation, giving radius and true anomaly (r, theta) as a function of time.

The problem here is that, in order to get r and theta as a function of time, I need to solve for the eccentric anomaly (psi) at a given time. Kepler's equation is:

psi = w*t + e*Sin(psi).

What you are calling "w" here is typically called the mean motion, Mdot or \dot M. "w" or \omega is typically reserved for angular velocity, that is, the time derivative of the true anomaly. The mean motion is a constant in the two body problem. Angular velocity is not constant. To get the true anomaly as a function of time, you do indeed need to solve Kepler's equation for the eccentric anomaly.

Newton's method is particularly quick even for highly eccentric orbits. Use the mean anomaly as an initial guess and iterate via Newton's method. Stop when the error reaches some sufficiently small value, such as 1e-12 radians.