Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Solving Kepler's equation using iterative technique

  1. Mar 12, 2008 #1
    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
    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!
  2. jcsd
  3. Mar 13, 2008 #2
    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):
    r(t) = x(t) - ae i,
    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:

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

    The values of the semi-major axis and eccentricity are obtained
    from the total energy (E=K/r + L2/(2mr2))
    and the angular momentum L, which are in turn specified by your
    initial conditions.
  4. Mar 13, 2008 #3
    Thanks very much for your help!
  5. Mar 13, 2008 #4

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    What you are calling "w" here is typically called the mean motion, Mdot or [itex]\dot M[/itex]. "w" or [itex]\omega[/itex] 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.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook