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

Keplers Equation

  1. Jun 17, 2008 #1
    Hi all,
    I am trying to find the eccentric anomaly E of an orbit of a binary system knowing the mean anomaly M through Newtons solution:
    E=M+e*sin(E)
    where an iterative solution can be found.
    I am using this to fit radial velocity data vs. time. But calculating E is hindering my progress (getting a good fit to data and length of time taken to calculate E). Is there another way to find E other than Newtons method?
     
  2. jcsd
  3. Jun 17, 2008 #2

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    The Newton-Raphson method for solving the reverse Kepler equation converges very quickly with a good initial guess; five steps even for very large (0.99) eccentricities. Some questions:
    • Are you sure you have implemented it correctly?

    • Are you truncating M to some 2*pi range? (0 to 2*pi is OK, -pi to pi is better).

    • What is your initial guess? A constant value (e.g., zero) may not even converge. An initial guess of [itex]E_0=M[/tex] is OK, but you can do much better. For example, [itex]E_0=M+e\sin M[/tex]. You can do even better, but this is simple and very good.

    • If timing is still a problem, return sin(E) and cos(E) as well as E. These are needed to calculate the solution to Kepler's equation and are needed to calculate the radial velocities. Use the returned values instead of calling sin and cos.
     
  4. Jun 17, 2008 #3
    I am using Levenberg-Marquardt algorithm to fit 5 orbital parameters of a binary exoplanet system. So the algorithm alters the period P and the time of periastron T when deducing a fit . I would have to think of how to implement a check when the program is altering those two parameters. let me know of any changes you may think would be wise for finding E.

    This is my code, hope this helps.
    ____________________________________
    function f= kepler(e,M)%where M=(2*pi/P)*(t-T)

    E0=0;
    E=M;
    while(abs(E-E0)>1.0e-6)

    E0=E;
    E= M + e*sin(E);%guess I should change this to : E=E+e*sin(E)
    end
    f=E;
    ______________________________________
     
  5. Jun 17, 2008 #4

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    Question number 1: Are you implementing it correctly? The answer is no. That algorithm is not using Newton-Raphson method. In general, the Newton-Raphson iterator for finding a zero of some function f(x) is

    [tex]x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}[/tex]

    You want to find the zero of [itex]f(E) = E-e\sin E - M[/itex] given some mean anomaly M. The Newton-Raphson iterator for the reverse Kepler problem is

    [tex]E_{n+1} = E_n - \frac{E_n-e\sin E_n - M}{1-e\cos E_n}[/tex]

    Your algorithm is incredibly slow to converge. A Newton-Raphson iterator tends to be much, much faster. Then key problem with any Newton-Raphson iterator is that it can shoot off to never-never land. This can happen here with large eccentricities (nearly 1) and small values of M. A little protection will stop that from happening.
     
  6. Jun 18, 2008 #5
    Thanks for the info D H. I have been reading up now on the Newton-Raphson method, and the methodology is beginning to make a bit more sense. My initial attempt is very wayward.
     
  7. Jun 18, 2008 #6
    Last edited: Jun 18, 2008
  8. Jun 20, 2008 #7
    I have implemented my new code D H. The run time is far superior than it was before
     
  9. Jun 25, 2008 #8
    Could you be a little more open about these protection issues D H ?
     
  10. Jun 25, 2008 #9

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    In this case, its pretty simple. The eccentric and mean anomalies are equal at 0 and pi. If, for example, the mean anomaly is between 0 and pi, an iteration that sends the eccentric anomaly outside [noparse][0,pi]/noparse] is wrong. Similarly, for a mean anomaly between -pi and 0 the eccentric anomaly must lie in [noparse][-pi,0]/noparse]. A step that goes outside the bounds or too many steps means something is wrong; you should revert to something that converges less quickly but is more stable, such as the secant method.
     
  11. Jun 26, 2008 #10
    I decided to check out the secant method you suggested since run time is not an issue at the
    moment. But my problem now lies in how to choose my starting points for the secant method. i.e. what should En and En-1 start at? Could you maybe point me towards a text etc?
     
  12. Jun 26, 2008 #11
    I have just read that E0 and E1 the initial starting points have to be sufficiently close to the root for the secant method to converge. But given my function f(E) how do ensure that generally with varying sets of data (namely differing time domains and orbital periods P)?
     
  13. Jun 27, 2008 #12
    Got it sorted.
    I Set:
    E0=M
    E1=M+e*sin(M)


    Hopefully this holds for varing P and t :yuck:
     
  14. Jun 27, 2008 #13

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    I tend to use Newton's method. Problems arise only for very large eccentricities (~0.9 or greater) and bad initial guesses.
     
  15. Jun 27, 2008 #14
    Do you use it over the Secant method because of the run time involved?
     
  16. Jun 27, 2008 #15

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    The secant method isn't guaranteed to converge, either. If you want something that converges you need use something real slow, like bisection. To find a root of some function, you start by finding a pair of points whose function values have opposite sign. Once you have found such a pair, evaluate the function at the midpoint. Throw out the point whose function value is of the same sign as this new function value. Repeat. Each step cuts the error by half. This is called linear convergence. It will take fifteen or so iterations with a function like Kepler's equation to find a root.
     
  17. Jun 29, 2008 #16
    There is an online calculator here too:

    Kepler Equation Calculator

    I think it uses a combination of bisection and secant method and usually converges to a solution in about 7 or 8 iterations.

    Making a good initial guess helps it converge fast too. For example, we know a solution exists between M and (M + e). If e is small, the search interval is pretty small to begin with. In addition, the starting interval is cut again before the iteration process begins.

    If you wanted to be really meticulous about the speed of program execution, I suppose you might want to consider Pade approximations or making tables.
     
  18. Jun 30, 2008 #17
    I have many data points so to use an online calculator would be impractical...
     
  19. Jun 30, 2008 #18
    Has anyone used the 'fzero' function in matlab?
     
  20. Jun 30, 2008 #19

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    Yes, but why bother? It is built to find the zero of any function. As a general purpose zero-finder it will be much slower than any reasonable special purpose function you can write to solve a function a very well-behaved function such as Kepler's equation.

    What level are you in in school?
     
  21. Jul 2, 2008 #20
    MSc, at the moment. Thanks for the info
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?