How to solve Kepler's equation $M=E-\varepsilon \sin E$ for $E$?

Question

I'm trying to create a program to solve a set of Kepler's Equation and I cannot isolate the single variable to use the expression in my program.

The Kepler Equation is $$M = E - \varepsilon \sin(E)$$

I will enter the values of $M$ and $\varepsilon$ and wish to find the value of $E$.

The website Wolfram Alpha could find a solution for this input $30(\frac{\pi}{180})=x-0.3 \sin(x)$

How can I find $E$?

EDIT:

I would like to propose this algorithm (Javascript) to solve the equation. It might require some adjustments depending on the programming language used. I did some basic tests with it, and would like feedback on it.

The following algorithm doesn't have an user input method and is treating the expression:

$0 = -x + \sin(-x)$

which is related to:

$x = \sin(-x)$

I used this algorithm to solve $E$ in Kepler's Equation rewriting it as (And swapping M and e for numerical constants, that in my case will be user inputs):

$0 = -M + E - e*sin(E)$

<script>
var start = Number.MIN_SAFE_INTEGER; end = Number.MAX_VALUE;

var result = null;

while(result === null) {
    var resultStart = Math.abs(-start + Math.sin(-start));
    var resultEnd = Math.abs(-end + Math.sin(-end));
    if(resultStart === 0) {
        result = start;
    } else if (resultEnd === 0) {
        result = end;
    } else {
        if (resultStart > resultEnd) {
            var startOld = start;
            start = (start+end)/2;
            if(startOld === start) { // underflow
                result = start;
            }
        } else {
            var endOld = end;
            end = (start+end)/2;
            if(endOld === end) { // underflow
                result = end;
            }
        }
    }
}

console.log(result);

</script>

Answer 1

There is no known closed form inverse for Kepler's equation. Newton's Method works well numerically, except in the case of near parabolic orbits ($\varepsilon\approx1$). It says to iterate $$ E_{n+1}=\frac{M+\varepsilon\sin(E_n)-\varepsilon E_n\cos(E_n)}{1-\varepsilon\cos(E_n)} $$ until it converges.

Answer 2

An analytical and elegant solution is the following Kapteyn / Fourier series:

$$E(M)=M+\sum_{n=1}\frac{ 2 J_n(n\epsilon) }{n} \sin(n M)$$

where $J_n(x)$ are the Bessel functions.

Other analytical solutions based on series expansion are discussed in:

"Solving Kepler's Equation Over Three Centuries", Peter Colwell, 1993

For implementation of Bessel functions in Javascript, several codes are available. For example:

http://www.mhtl.uwaterloo.ca/old/courses/me3532/js/bessel.html

https://github.com/SheetJS/bessel

Answer 3

$\def\L{\operatorname L} \def\J{\operatorname J}$

@giorgiomugnaini gave a Fourier series solution and the following integral representations was not found elsewhere online, but here are other solutions to the Kepler equation

Laplace Transform:

We apply the Laplace inversion theorem and one could expand $\sin(\phi)$, like in the Fourier series derivation, or:

$$f(\phi)=\sin(\phi)+a\phi=x\implies \phi=\L^{-1}_x(\L_s(\phi(t))$$

We define $\phi(x)=0\not<x\not<\pi a$, substitute $\phi(t)\to t$, and integrate by parts:

$$\L_s(\phi(t))=\int_0^{\pi a}e^{-st}\phi(t)dt=\int_0^\pi e^{-sf(t)}t\ df(t)=-\frac\pi s e^{-\pi a s}+\int_0^\pi e^{-s\sin(t)-a s t}dt$$

$-\frac\pi s e^{-\pi a s}$’s inverse Laplace transform is a Heaviside theta function which is zero for $x<\pi a$. Alternatively, if $f(x)$ is monotonic, then there is global inverse, so another setup is:

$$\L_s(h(t))=\int_0^\infty e^{-st}dt=\int_0^\infty t e^{-sf(t)}df(t)=\frac1s\int_0^\infty e^{-s\sin(t)-a s t}dt$$

Finally, expand $e^{c\sin(x)}$ via Maclaurin series, integrate, and sum over to obtain an Anger J. Therefore:

$$\bbox[3px,border: solid blue 3px]{\sin(\phi)+a\phi=x\implies \phi=\int_{ic-\infty}^{ic+\infty}\frac{e^{isx}\J_{a s}(s)}{is(e^{\pi i a s}+1)}ds,\int_{ic+\infty}^{ic-\infty}\frac{e^{isx}}{2s}\csc(\pi a s)\J_{as}(s)ds}$$

The first integral represents a period of $\phi(x)$ while the second is the global inverse. shown here:

Fourier Transform:

The Fourier inversion theorem uses Anger J:

$$\bbox[3px,border: solid blue 3px]{\sin(\phi)-a\phi=x\implies \phi=\int_0^\infty\frac2t\sin(tx)\operatorname J_{at}(t)dt}$$

shown here:

Power Series:

$\frac{\sin(\phi)+a\phi}{a+1}=x$, $\phi(x)\approx x$ satisfies $$a\phi^2\phi’+2(a+1)\Phi\phi’-2(a+1)x\phi\phi’+2a\phi’=2(a+1),y(0)=0,y’(0)=1$$ where $\Phi=\int_0^x\phi(t)dt$. One notices and substitutes $\phi=\sum\limits_{n=0}^\infty a_nx^{2n+1}$. Then, gathering powers of $x$ to find a recurrence relation for $a_n$ and letting $x\to\frac x{a+1}$ gives:

$$\sin(\phi)+a\phi=x\implies \phi=\sum_{n=0}^\infty a_n\left(\frac x{a+1}\right)^{2n+1}\\2(a+1)(2n+1)a_{n+1}=\sum_{k=0}^{n-1}(a+1)a_ka_{n-k-1}(2(n-k)-1)\frac{2k+1}{k+1}-\sum_{m=0}^{n-k-1}\sum_{k=0}^n aa_ka_ma_{n-k-m-1}(2m+1)\\a_n=\left\{1,\frac1{6(a+1)},\frac{9-a}{120(a+1)^2 },\frac{225-54a+a^2 }{5040(a+1)^3},\frac{11025-4131a+243a^2-a^3}{362880(a+1)^4},\frac{893025-457200a+50166a^2-1008a^3+a^4}{39916800(a+1)^5},\frac{108056025-70301925a+11708154a^2-520218a^3+4077a^4-a^5}{6227020800(a+1)^6},\dots\right\}$$

This matches the series reversion of $\sin(x)+ax$ and the recurrence relation can be tested with -1/(2(a+1)(2n+1))Sum[Sum[aa[k]a[n-k-j-1] a[j] (2j+1),{j,0,n-k-1}],{k,0,n}]-Sum[(a+1) a[k] a[n-k-1](2(n-k-1)+1)(2k+1)/(k+1),{k,0,n-1}] It seems that there is no single sum recurrence relation except for $a=1$ to invert $\sin^{-1}(x)+x\sqrt{1-x^2}$. The series can take $a,x\in\Bbb C$ as shown here:

Answer 4

For $E\neq 0$, Kepler's equation can be solved by Hyper Lambert W:

$$E-\varepsilon\sin(E)=M$$ $$E\ e^{\ln(1-\varepsilon\frac{\sin(E)}{E})}=M$$ $$E=HW\left(\left\{\ln(1-\varepsilon\frac{\sin(x)}{x})\right\}_1;M\right)$$

So we have a closed form for $E$, and the series representations of Hyper Lambert W give some hints for calculating $E$.

[Galidakis 2005] Galidakis , I. N.: On solving the p-th complex auxiliary equation $f^{(p)}(z)=z$. Complex Variables 50 (2005) (13) 977-997

Galidakis, I. N.: On some applications of the generalized hyper-Lambert functions. Complex Variables and Elliptic Equations 52 (2007) (12) 1101-1119

[Dubinov/Galidakis 2007] Dubinov, A.; Galidakis, Y.: Explicit solution of the Kepler equation. Physics of Particles and Nuclei Letters 4 (2007) 213-216

[Masson] Masson, Paul: Analytic Physics - An Exact Solution of the Complex Kepler Equation

Answer 5

Let $M>0$, and let's take the change of variable $E\rightarrow\left(1-\varepsilon\right)\log\left(1-1/z\right)$ then the equation becomes

\begin{equation} \left(1-\varepsilon\right)\log\left(1-1/z\right)-\varepsilon\sin\left(\left(1-\varepsilon\right)\log\left(1-1/z\right)\right)-M=0 \end{equation}

Multiplying the equation by z, we define the function

\begin{equation} f(z)=z\left(1-\varepsilon\right)\log\left(1-1/z\right)-\varepsilon z\sin\left(\left(1-\varepsilon\right)\log\left(1-1/z\right)\right)-Mz \end{equation}

and applying the Burniston-Siewert method to solve transcendental equations to this equation we obtain

\begin{equation} E=\log\left(1-1/z_{0}\right),\quad z_{0}=\frac{M-2\left(1-\varepsilon\right)^{2}}{2M}+m \end{equation}

where

\begin{equation} m=\frac{1}{\pi}{\displaystyle \int\limits _{0}^{1}\arctan\left(\frac{\left(1-\varepsilon\right)\log\left(\frac{1}{t}-1\right)-M-\varepsilon R\left(t\right)}{\pi\left(1-\varepsilon\right)-\varepsilon I\left(t\right)}\right)\,dt} \end{equation}

and

\begin{equation} R\left(t\right)=\textrm{Re}\left[\sin\left(\left(1-\varepsilon\right)\log\left(1/t-1\right)+i\left(1-\varepsilon\right)\pi\right)\right] \end{equation}

\begin{equation} I\left(t\right)=\textrm{Im}\left[\sin\left(\left(1-\varepsilon\right)\log\left(1/t-1\right)+i\left(1-\varepsilon\right)\pi\right)\right] \end{equation}

The expanded form of these functions is

\begin{equation} R\left(t\right)=\sin\left[\left(1-\varepsilon\right)\log\left(1/t-1\right)\right]\cosh\left[\left(1-\varepsilon\right)\pi\right] \end{equation}

\begin{equation} I\left(t\right)=\cos\left[\left(1-\varepsilon\right)\log\left(1/t-1\right)\right]\sinh\left[\left(1-\varepsilon\right)\pi\right] \end{equation}

Putting everything into a single expression, we get

where the integral involved is completely real.

Or

Denoting by

\begin{equation} K\left(E\right)=E-\varepsilon\sin\left(E\right) \end{equation}

An inverse graph for $\varepsilon=1/2$ will be

and taking $K^{-1}(K(E))$

Furthermore, the expression found gives us an approximation when $\left|M\right|\rightarrow0$

$$ \begin{align} E & \sim \left(1-\varepsilon\right)\log\left[1-\left[\frac{M-2\left(1-\varepsilon\right)^{2}}{2M}+\frac{1}{\pi}{\displaystyle \int\limits _{0}^{1}\frac{\left(1-\varepsilon\right)\log\left(\frac{1}{t}-1\right)-M-\varepsilon R\left(t\right)}{\pi\left(1-\varepsilon\right)-\varepsilon I\left(t\right)}\,dt}\right]^{-1}\right] \\ & \sim\left(1-\varepsilon\right)\log\left[1-\left[\frac{M-2\left(1-\varepsilon\right)^{2}}{2M}+\frac{1}{\pi}{\displaystyle \int\limits _{0}^{1}\frac{-M}{\pi\left(1-\varepsilon\right)-\varepsilon I\left(t\right)}\,dt}\right]^{-1}\right] \\ & \sim\left(1-\varepsilon\right)\log\left[1-\left[\frac{M-2\left(1-\varepsilon\right)^{2}}{2M}-\frac{M}{\pi}c\right]^{-1}\right] \\ & \sim\left(1-\varepsilon\right)\log\left(1-\frac{2\pi M}{M-2\left(1-\varepsilon\right)^{2}-2cM^{2}}\right) \end{align} $$

then

\begin{equation} E\sim\left(1-\varepsilon\right)\log\left(\frac{2cM^{2}+2\pi\left(1-\varepsilon\right)^{2}+\pi x}{2cM^{2}+2\pi\left(1-\varepsilon\right)^{2}-\pi x}\right){,\quad\left|M\right|\rightarrow0} \end{equation}

Ref:

-E. E. Burniston, C.E. Siewert. The use of Riemann problems in solving a class of transcendental equations. Mathematical Proceedings of the Cambridge Philosophical Society. 73. 111 - 118. (1973).

-Henrici P., Applied and computational complex analysis, Volume III. Wiley, New York. 183-191 (1986).

Answer 6

$$M=E-\varepsilon\sin(E)=(1-\varepsilon)E-\varepsilon\sum_{n=1}^\infty (-1)^n\,\frac{E^{2 n+1}}{(2 n+1)!}\tag 1$$ Using $\alpha=(1-\varepsilon)$ and power series reversion $$E=\frac M \alpha-(\alpha-1)\sum_{n=1}^\infty\frac{M^{2n+1}}{(2 n+1)!\,\,\alpha^{3n+1} }\,P_n(\alpha)\tag 2$$ where the first polynomials are $$\left( \begin{array}{cc} n & P_n(\alpha) \\ 1 & 1 \\ 2 & 9 \alpha -10 \\ 3 & 225 \alpha ^2-504 \alpha +280 \\ 4 & 11025 \alpha ^3-37206 \alpha ^2+41580 \alpha -15400 \\ 5 & 893025 \alpha ^4-4029300 \alpha ^3+6779916 \alpha ^2-5045040 \alpha +1401400 \\ \end{array} \right)$$

Notice that, knowing the coefficients of $(1)$, we also know all coefficients of $(2)$.

Repeating @Tyma Gaidash' example $\left(M=2,\varepsilon=-\frac {10}{23}\right)$ and only the above terms, the result is $E=1.565231$ while the solution is $1.565224$.