# Taylor integration of the Kepler problem

Here, we try to reproduce __exactly__ the [Kepler problem integration example](http://nbviewer.jupyter.org/github/JuliaDiff/TaylorSeries.jl/blob/master/examples/1-KeplerProblem.ipynb) made by Luis Benet in [JuliaDiff/TaylorSeries.jl](https://github.com/JuliaDiff/TaylorSeries.jl).

The Kepler problem is the basis of planetary motion; it describes the motion of a secondary body (e.g., a planet, asteroid, comet, etc.), around a primary body (e.g., the Sun). In cartesian coordinates over the orbital plane, the Hamiltonian for the Kepler problem reads:

$$
\begin{align}
H_{\mathrm{Kepler}} &= \frac{1}{2\mu}(p_x^2+p_y^2)-\frac{\mu}{\sqrt{x^2+y^2}}
\end{align}
$$

where $\mu=G(m_1+m_2)$, $G$ is the gravitational constant, $m_1$ is the mass of the primary body and $m_2$ is the mass of the secondary body. If we write $\vec r_1 = (x_1,y_1)$ and $\vec r_2 = (x_2,y_2)$, respectively, for the position of the primary and secondary body, then the vector $\vec r = \vec r_2-\vec r_1$, with coordinates $\vec r = (x, y)=(x_2-x_1,y_2-y_1)$, represents the position of the secondary body relative to the primary body; i.e., $\vec r$ represents the so-called relative coordinates. In terms of the vector $\vec{r}$, the position $\vec{r}=0$ corresponds to the position of the primary body.

Using Hamilton equations, we can obtain the equations of motion for the Kepler problem:

$$
\begin{align}
\dot x &= u \\
\dot y &= v \\
\dot u &= -\frac{\mu x}{(x^2+y^2)^{3/2}}\\
\dot v &= -\frac{\mu y}{(x^2+y^2)^{3/2}}
\end{align}
$$

Note that the canonical momenta are $p_x$ and $p_y$, while $u$ and $v$ are, respectively, the  $x$ and $y$ components of the velocity; i.e., $p_x = \mu u$ and $p_y = \mu v$. 


First of all, we shall include all relevant packages:


```julia
using TaylorIntegration, Plots, LaTeXStrings
pyplot()
```

    Fontconfig warning: ignoring UTF-8: not a valid region tag





    Plots.PyPlotBackend()



Some parameters necessary for the integration:

+ $\mu$: the gravitational parameter
+ `q0`: the initial condition (we will select an initial condition which corresponds to elliptical motion)
+ `order`: the order of the Taylor expansion
+ `t_max`: the final time of the integration
+ `abs_tol`: the absolute tolerance
+ `n_iter`: the number of time-steps


```julia
const μ = 1.0
const q0 = [0.19999999999999996, 0.0, 0.0, 3.0] # a initial condition for elliptical motion
const order = 28
const t0 = 0.0
const t_max = 10000*(2π) # we are just taking a wild guess about the period ;)
const abs_tol = 1.0E-20
const steps = 500000
```




    500000



As usual, we write down the equations of motion into a `function`, which here we will name `kepler!`. `q` represents the system state; i.e., the set of values of the dynamical variables at a given instant; `dq` represents the time derivatives of the components of `q`.


```julia
#an auxiliary array which helps with optimization:
const r_p3d2 = Array{TaylorSeries.Taylor1{Float64}}(1)

#the equations of motion for the Kepler problem:
function kepler!(t, q, dq)
    r_p3d2[1] = (q[1]^2+q[2]^2)^(3/2)
    
    dq[1] = q[3]
    dq[2] = q[4]
    dq[3] = -μ*q[1]/r_p3d2[1]
    dq[4] = -μ*q[2]/r_p3d2[1]
    
    nothing
end
```




    kepler! (generic function with 1 method)



The Taylor integration:


```julia
t, q = taylorinteg(kepler!, q0, t0, 0.01, order, abs_tol, maxsteps=2); #warm-up lap
@time t, q = taylorinteg(kepler!, q0, t0, t_max, order, abs_tol, maxsteps=steps);
```

     43.060857 seconds (464.94 M allocations: 39.523 GB, 12.87% gc time)


The final state:


```julia
t[end]/(2pi), q[end,:]
```




    (10000.0,[0.2,2.72854e-9,-2.27378e-8,3.0])



Let's extract the values of $x$, $y$, $u$ and $v$ for each time-step:


```julia
x, y, u, v = view(q,:,1), view(q,:,2), view(q,:,3), view(q,:,4);
```

## Orbital motion

The initial conditions we selected correspond to a elliptical orbit with a relatively high eccentricity: $e=0.8$ for the initial condition `q0=[0.19999999999999996, 0.0, 0.0, 3.0]`. How does the motion of the planet/asteroid/comet looks like? Well, let's plot its orbit over the $x-y$ plane (the orbit is shown in blue; the position of the primary body is shown as a yellow dot):


```julia
scatter(

[0.0], [0.0],
label = L"\mathrm{Center\; of\; attraction}",
ms=10

)
scatter!(

x[1:75:end], y[1:75:end],
title = L"\mathrm{Orbital\;\;motion}", 
xaxis = (L"x", (-2.0, 0.5), -2.0:0.5:2.0), 
yaxis = (L"y", (-0.8, 0.8), -2.0:0.5:2.0),
label = L"\mathrm{Keplerian\; ellipse}",
ms=0.5

)
```








## Energy conservation

Below, we write the energy function of the Kepler problem; i.e., the Hamiltonian $H_{\mathrm{Kepler}}$ in terms of $x$, $y$, $u$ and $v$:


```julia
H(x_, y_, u_, v_) = 0.5*(u_^2+v_^2)-μ/sqrt(x_^2+y_^2)
```




    H (generic function with 1 method)



Now, using the function `H` defined above, we calculate the energy during each time-step:


```julia
E = H.(x, y, u, v);
```

$E_0=H_{\mathrm{Kepler}}(x_0,y_0,u_0,v_0)$, the initial value of the energy, is:


```julia
E0=E[1]
```




    -0.5000000000000009



We define $\delta E$ as the relative error in the energy; i.e.:

$$
\begin{align}
\delta E(t) &= \frac{E(t)-E_0}{E_0}
\end{align}
$$


```julia
δE = (E-E0)/(E0)
```




    436153-element Array{Float64,1}:
     -0.0        
      1.77636e-15
      1.77636e-15
      8.88178e-16
      8.88178e-16
      2.66454e-15
      3.55271e-15
      2.66454e-15
      3.55271e-15
      2.66454e-15
      3.9968e-15 
      2.66454e-15
      3.10862e-15
      ⋮          
     -7.12763e-14
     -7.19425e-14
     -7.14984e-14
     -7.10543e-14
     -7.19425e-14
     -7.10543e-14
     -7.10543e-14
     -6.92779e-14
     -7.01661e-14
     -6.92779e-14
     -6.92779e-14
     -7.10543e-14



The analytical solution preserves the energy; i.e., $\delta E(t)=0$ for the analytical solution. Thus, we expect our solution to be close to zero. So, even if the $\delta E$ is not perfectly zero during the whole integration, it is comparable to Julia's `Float64` machine-epsilon.

Below, we plot $\delta E$ in units of Julia's `Float64` machine-epsilon as a function of time, $t$. In Julia, the machine-epsilon for the `Float64` type has a value


```julia
eps(Float64)
```




    2.220446049250313e-16




```julia
plot(

t/(2π), δE/eps(Float64),
title = L"\mathrm{Energy\;\; relative\;\; error\;\; vs\;\; time}", 
xaxis = (L"t\mathrm{\,(orbital\;\;periods)}"),
yaxis = (L"\delta E \mathrm{\;\;(machine\;\;epsilons)}"),
label = L"\mathrm{Energy\;\; relative\;\; error\;\; vs\;\; time}"

)
```








Now, how does the energy error distribute around zero?


```julia
histogram(

δE/eps(Float64),
title = L"\mathrm{Distribution\;\;of\;\;energy\;\;relative\;\;error}", 
xaxis = (L"\delta E"),
yaxis = (L"\mathrm{Number\;\;of\;\;}\delta E\mathrm{\;\;values\;\;within\;\;bin\;\;range}"),
leg = false,
nbins=20

)
```

    [1m[34mINFO: binning = 20
    [0m








The energy error behaves roughly as a random walk, which means that the numerical error in the energy is dominated by rounding errors.

## Angular momentum conservation

Analogously to the energy, we now focus on the angular momentum, which is preserved by the analytical solution too. The value of the angular momentum is given by
$$
\begin{align}
L &= xv-yu
\end{align}
$$
We write the angular momentum function as:


```julia
ang_mom(x_, y_, u_, v_) = x_.*v_-y.*u_
```




    ang_mom (generic function with 1 method)



So the angular momentum during each time-step is:


```julia
L = ang_mom(x, y, u, v);
```

$L_0=x_0v_0-y_0u_0$, the initial value of the angular momentum, is:


```julia
L0 = L[1]
```




    0.5999999999999999



We define $\delta L$ as the relative error in the angular momentum; i.e.:

$$
\begin{align}
\delta L(t) &= \frac{L(t)-L_0}{L_0}
\end{align}
$$


```julia
δL = (L-L0)/L0;
```

Just as we did for the energy, we now plot $\delta L$ in units of `eps(Float64)` vs $t$:


```julia
plot(

t/(2π), δL/eps(Float64),
title = L"\mathrm{Angular\;\;momentum\;\; relative\;\; error\;\; vs\;\; time}", 
xaxis = (L"t\mathrm{\,(orbital\;\;periods)}"),
yaxis = (L"\delta L \mathrm{\;\;(machine\;\;epsilons)}"),
label = L"\mathrm{Angular\;\;momentum\;\; relative\;\; error\;\; vs\;\; time}",
leg=false,
color=:orange

)
```








Again, we see that the angular momentum relative error is comparable to `eps(Float64)`; the maximum variation is ~$150$ machine epsilons. How does the distribution of this error look like?


```julia
histogram(

δL/eps(Float64),
title = L"\mathrm{Distribution\;\;of\;\;angular\;\;momentum\;\;relative\;\;error}", 
xaxis = (L"\delta L"),
yaxis = (L"\mathrm{Number\;\;of\;\;}\delta L\mathrm{\;\;values\;\;within\;\;bin\;\;range}"),
leg = false,
nbins=20,
color=:orange

)


```

    [1m[34mINFO: binning = 20
    [0m








The distribution of $\delta L$ shows a peak near zero and roughly symmetrical around this value. This means that the error in the angular momentum is dominated also by rounding errors of floating-point arithmetic.

Lastly, we reproduce __exactly__ the last plot shown in the original version of this example, authored by Luis Benet and included in `JuliaDiff/TaylorSeries.jl`'s [Kepler problem integration example](http://nbviewer.jupyter.org/github/JuliaDiff/TaylorSeries.jl/blob/master/examples/1-KeplerProblem.ipynb) jupyter notebook.

A $\delta E$, $\delta L$ plot vs $t$:


```julia
plot(

[t/(2π), t/(2π)], [δE/eps(Float64), δL/eps(Float64)],
title = L"\delta E\mathrm{\;\;(blue),\;\;}\delta L\mathrm{\;\;(green)\;\;vs\;\;time}", 
xaxis = (L"t\mathrm{\,(orbital\;\;periods)}"),
yaxis = (L"\delta E\mathrm{,\;\;}\delta L\;\;\mathrm{(machine\;\;epsilons)}"),
leg=false,

)
```









```julia

```

dff4432ce359d1fd69fc00eecf9675eecd7c0378
ipynb
Jupyter Notebook
examples/Kepler-problem.ipynb
SebastianM-C/TaylorIntegration.jl
f3575ee1caba43e21312062d960613ec2ccba325
MIT
2016-09-22T22:32:12.000Z
2022-03-23T13:35:18.000Z
examples/Kepler-problem.ipynb
SebastianM-C/TaylorIntegration.jl
f3575ee1caba43e21312062d960613ec2ccba325
MIT
2016-09-21T05:43:08.000Z
2022-03-15T02:55:17.000Z
examples/Kepler-problem.ipynb
SebastianM-C/TaylorIntegration.jl
f3575ee1caba43e21312062d960613ec2ccba325
MIT
2016-09-24T04:37:11.000Z
2022-03-25T13:48:07.000Z
Qwen/Qwen-72B
1. YES
2. YES