Keplers Equation: Finding Eccentric Anomaly E

[NoobixCube]

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?

 

[D H]

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 E_0=M[/tex] is OK, but you can do much better. For example, E_0=M+e\sin M[/tex]. You can do even better, but this is simple and very good.<br /> <br /> <br /> [*]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.

 

[NoobixCube]

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;
______________________________________

 

[D H]

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

x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}

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

E_{n+1} = E_n - \frac{E_n-e\sin E_n - M}{1-e\cos E_n}

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.

 

[NoobixCube]

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.

 

[Helios]

Ya know, I once saw an online calculator for this sort of thing. Just plug in the numbers and click. But now I can't find it. I bet it's still online somewhere.

Yea, like here
http://www.geocities.com/SiliconValley/2902/kepler.htm

 

[NoobixCube]

I have implemented my new code D H. The run time is far superior than it was before

 

[NoobixCube]

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.

Could you be a little more open about these protection issues D H ?

 

[D H]

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.

 

[NoobixCube]

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 E_n and E_n-1 start at? Could you maybe point me towards a text etc?

 

[NoobixCube]

I have just read that E₀ and E₁ 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)?

 

[NoobixCube]

I have just read that E₀ and E₁ 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)?

Got it sorted.
I Set:
E₀=M
E₁=M+e*sin(M)


Hopefully this holds for varing P and t

 

[D H]

I tend to use Newton's method. Problems arise only for very large eccentricities (~0.9 or greater) and bad initial guesses.

 

[NoobixCube]

I tend to use Newton's method. Problems arise only for very large eccentricities (~0.9 or greater) and bad initial guesses.

Do you use it over the Secant method because of the run time involved?

 

[D H]

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.

 

[DuncanM]

There is an online calculator here too:

http://www.akiti.ca/KeplerSolver.html"

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.

 

[NoobixCube]

There is an online calculator here too:

http://www.akiti.ca/KeplerSolver.html"

I have many data points so to use an online calculator would be impractical...

 

[NoobixCube]

Has anyone used the 'fzero' function in matlab?

 

[D H]

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?

 

[NoobixCube]

MSc, at the moment. Thanks for the info