Keplers Equation

  • Thread starter NoobixCube
  • Start date
  • #1
155
0

Main Question or Discussion Point

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?
 

Answers and Replies

  • #2
D H
Staff Emeritus
Science Advisor
Insights Author
15,393
682
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.
 
  • #3
155
0
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;
______________________________________
 
  • #4
D H
Staff Emeritus
Science Advisor
Insights Author
15,393
682
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.
 
  • #5
155
0
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.
 
  • #6
249
38
Last edited:
  • #7
155
0
I have implemented my new code D H. The run time is far superior than it was before
 
  • #8
155
0
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 ?
 
  • #9
D H
Staff Emeritus
Science Advisor
Insights Author
15,393
682
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.
 
  • #10
155
0
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?
 
  • #11
155
0
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)?
 
  • #12
155
0
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)?
Got it sorted.
I Set:
E0=M
E1=M+e*sin(M)


Hopefully this holds for varing P and t :yuck:
 
  • #13
D H
Staff Emeritus
Science Advisor
Insights Author
15,393
682
I tend to use Newton's method. Problems arise only for very large eccentricities (~0.9 or greater) and bad initial guesses.
 
  • #14
155
0
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?
 
  • #15
D H
Staff Emeritus
Science Advisor
Insights Author
15,393
682
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.
 
  • #16
98
2
There is an online calculator here too:

http://www.akiti.ca/KeplerSolver.html" [Broken]

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.
 
Last edited by a moderator:
  • #17
155
0
There is an online calculator here too:

http://www.akiti.ca/KeplerSolver.html" [Broken]
I have many data points so to use an online calculator would be impractical...
 
Last edited by a moderator:
  • #18
155
0
Has anyone used the 'fzero' function in matlab?
 
  • #19
D H
Staff Emeritus
Science Advisor
Insights Author
15,393
682
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?
 
  • #20
155
0
MSc, at the moment. Thanks for the info
 

Related Threads for: Keplers Equation

  • Last Post
Replies
1
Views
799
  • Last Post
Replies
4
Views
1K
  • Last Post
Replies
9
Views
3K
Replies
2
Views
6K
Replies
17
Views
2K
Replies
4
Views
1K
Replies
2
Views
2K
Top