Interoperability with DifferentialEquations.jl

Interoperability with DifferentialEquations.jl

Here, we show an example of interoperability between TaylorIntegration.jl and DifferentialEquations.jl, i.e., how to use TaylorIntegration.jl from the DifferentialEquations ecosystem. The basic requirement is to load DiffEqBase.jl, which sets-up the common interface. Below, we shall also use OrdinaryDiffEq.jl to compare the accuracy of TaylorIntegration.jl with respect to high-accuracy methods for non-stiff problems (Vern9 method). While DifferentialEquations offers many macros to simplify certain aspects, we do not rely on them simply because using properly @taylorize improves the performance.

Note

Currently, the only keyword argument supported by DiffEqBase.solve that is implemented in TaylorIntegration.jl is :saveat. The keyword argument :parse_eqs is available in order to control the use of methods defined via @taylorize.

The problem we will integrate in this example is the planar circular restricted three-body problem (PCR3BP, also capitalized as PCRTBP). The PCR3BP describes the motion of a body with negligible mass $m_3$ under the gravitational influence of two bodies with masses $m_1$ and $m_2$, such that $m_1 \ge m_2$. It is assumed that $m_3$ is much smaller than the other two masses so it does not influence their motion, and therefore it is simply considered as a massless test particle. The body with the greater mass $m_1$ is referred as the primary, and $m_2$ as the secondary. These bodies are together called the primaries and are assumed to describe Keplerian circular orbits about their center of mass, which is placed at the origin of the reference frame. It is further assumed that the orbit of the third body takes place in the orbital plane of the primaries. A full treatment of the PCR3BP may be found in [1].

The ratio $\mu = m_2/(m_1+m_2)$ is known as the mass parameter. Using mass units such that $m_1+m_2=1$, we have $m_1=1-\mu$ and $m_2=\mu$. In this example, we assume the mass parameter to have a value $\mu=0.01$.

using Plots

const μ = 0.01

The Hamiltonian for the PCR3BP in the synodic frame (i.e., a frame which rotates such that the primaries are at rest on the $x$ axis) is

\[\begin{equation} \label{eq-pcr3bp-hamiltonian} H(x, y, p_x, p_y) = \frac{1}{2}(p_x^2+p_y^2) - (x p_y - y p_x) + V(x, y) \end{equation}\]

where

\[\begin{equation} \label{eq-pcr3bp-potential} V(x, y) = - \frac{1-\mu}{\sqrt{(x-\mu)^2+y^2}} - \frac{\mu}{\sqrt{(x+1-\mu)^2+y^2}}. \end{equation}\]

is the gravitational potential associated to the primaries. The RHS of Eq. (\ref{eq-pcr3bp-hamiltonian}) is also known as the Jacobi constant, since it is a preserved quantity of motion in the PCR3BP. We will use this property to check the accuracy of the solutions computed.

V(x, y) = - (1-μ)/sqrt((x-μ)^2+y^2) - μ/sqrt((x+1-μ)^2+y^2)
H(x, y, px, py) = (px^2+py^2)/2 - (x*py-y*px) + V(x, y)
H(x) = H(x...)

The equations of motion for the PCR3BP are

\[\begin{eqnarray} \label{eqs-motion-pcr3bp} \dot{x} &=& p_x + y \\ \dot{y} &=& p_y - x \\ \dot{p_x} &=& - \frac{(1-\mu)(x-\mu)}{((x-\mu)^2+y^2)^{3/2}} - \frac{\mu(x+1-\mu)}{((x+1-\mu)^2+y^2)^{3/2}} + p_y \\ \dot{p_y} &=& - \frac{(1-\mu)y }{((x-\mu)^2+y^2)^{3/2}} - \frac{\mu y }{((x+1-\mu)^2+y^2)^{3/2}} - p_x. \end{eqnarray}\]

We define this system of ODEs using the most naive approach

function f(dq, q, param, t)
    local μ = param[1]
    x, y, px, py = q
    dq[1] = px + y
    dq[2] = py - x
    dq[3] = - (1-μ)*(x-μ)*((x-μ)^2+y^2)^-1.5 - μ*(x+1-μ)*((x+1-μ)^2+y^2)^-1.5 + py
    dq[4] = - (1-μ)*y    *((x-μ)^2+y^2)^-1.5 - μ*y      *((x+1-μ)^2+y^2)^-1.5 - px
    return nothing
end

Note that DifferentialEquations offers interesting alternatives to write these equations of motion in a simpler and more convenient way, for example, using the macro @ode_def, see ParameterizedFunctions.jl. We have not used that flexibility here because TaylorIntegration.jl has @taylorize, which under certain circumstances allows to important speed-ups.

We shall define the initial conditions $q_0 = (x_0, y_0, p_{x,0}, p_{y,0})$ such that $H(q_0) = J_0$, where $J_0$ is a prescribed value. In order to do this, we select $y_0 = p_{x,0} = 0$ and compute the value of $p_{y,0}$ for which $H(q_0) = J_0$ holds.

We consider a value for $J_0$ such that the test particle is able to display close encounters with both primaries, but cannot escape to infinity. We may obtain a first approximation to the desired value of $J_0$ if we plot the projection of the zero-velocity curves on the $x$-axis.

ZVC(x) =  -x^2/2 + V(x, zero(x)) # projection of the zero-velocity curves on the x-axis

plot(ZVC, -2:0.001:2, label="zero-vel. curve", legend=:topleft)
plot!([-2, 2], [-1.58, -1.58], label="J0 = -1.58")
ylims!(-1.7, -1.45)
xlabel!("x")
ylabel!("J")
title!("Zero-velocity curves (x-axis projection)")
-2 -1 0 1 2 -1.70 -1.65 -1.60 -1.55 -1.50 -1.45 Zero-velocity curves (x-axis projection) x J zero-vel. curve J0 = -1.58

Notice that the maxima in the plot correspond to the Lagrangian points $L_1$, $L_2$ and $L_3$; below we shall concentrate in the value $J_0 = -1.58$.

J0 = -1.58

We define a function py!, which depends on the initial condition $q_0 = (x_0, 0, 0, p_{y,0})$ and the Jacobi constant value $J_0$, such that it computes an adequate value $p_{y,0}$ for which we have $H(q_0)=J_0$ and updates (in-place) the initial condition accordingly.

function py!(q0, J0)
    @assert iszero(q0[2]) && iszero(q0[3]) # q0[2] and q0[3] have to be equal to zero
    q0[4] = q0[1] + sqrt( q0[1]^2-2( V(q0[1], q0[2])-J0 ) )
    nothing
end

We are now ready to generate an appropriate initial condition.

q0 = [-0.8, 0.0, 0.0, 0.0]
py!(q0, J0)
q0
4-element Array{Float64,1}:
 -0.8
  0.0
  0.0
 -0.6276410653920694

We note that the value of q0 has been updated. We can check that the value of the Hamiltonian evaluated at the initial condition is indeed equal to J0.

H(q0) == J0
true

Following the DifferentialEquations.jl tutorial, we define an ODEProblem for the integration; TaylorIntegration.jl can be used via its common interface bindings with DiffEqBase.jl; both packages need to be loaded explicitly.

tspan = (0.0, 1000.0)
p = [μ]

using TaylorIntegration, DiffEqBase
prob = ODEProblem(f, q0, tspan, p)
WARNING: using DataStructures.update! in module TaylorIntegration conflicts with an existing identifier.
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 1000.0)
u0: [-0.8, 0.0, 0.0, -0.627641]

We solve prob using a 25-th order Taylor method, with a local absolute tolerance $\epsilon_\mathrm{tol} = 10^{-20}$.

solT = solve(prob, TaylorMethod(25), abstol=1e-20);
retcode: Success
Interpolation: 1st order linear
t: 13782-element view(::Array{Float64,1}, 1:13782) with eltype Float64:
    0.0
    0.08167917432844125
    0.1648062016344043
    0.25194621452252675
    0.34542063977828563
    0.44783200479034113
    0.5596562661907026
    0.6551542940610247
    0.7361173484102514
    0.8038909079882202
    ⋮
  999.5057885626638
  999.6134222959215
  999.7037150677216
  999.780461965901
  999.8461067752269
  999.902790056689
  999.9521071897881
  999.9949632045395
 1000.0
u: 13782-element Array{Array{Float64,1},1}:
 [-0.8, 0.0, 0.0, -0.627641]
 [-0.797405, 0.013901, 0.0497572, -0.631547]
 [-0.789379, 0.0269558, 0.102796, -0.643337]
 [-0.774966, 0.0382713, 0.16329, -0.663658]
 [-0.75238, 0.0463273, 0.236125, -0.694037]
 [-0.718664, 0.0485296, 0.328746, -0.737382]
 [-0.67018, 0.0406328, 0.45187, -0.797005]
 [-0.617795, 0.0232668, 0.584923, -0.85796]
 [-0.563845, -0.00059466, 0.72977, -0.915361]
 [-0.510224, -0.027846, 0.886601, -0.963953]
 ⋮
 [-0.365033, -0.547722, 1.03508, -0.158738]
 [-0.30456, -0.523352, 1.16406, -0.0555769]
 [-0.240076, -0.498595, 1.28972, 0.0636653]
 [-0.173914, -0.472613, 1.40809, 0.205528]
 [-0.108177, -0.444571, 1.51278, 0.3743]
 [-0.0443663, -0.413652, 1.59576, 0.572798]
 [0.0161301, -0.379313, 1.64685, 0.800576]
 [0.0715682, -0.34171, 1.65495, 1.05024]
 [0.0781893, -0.336715, 1.65243, 1.08298]

As mentioned above, we load OrdinaryDiffEq in order to solve the same problem prob now with the Vern9 method, which the DifferentialEquations.jl documentation recommends for high-accuracy (i.e., very low tolerance) integrations of non-stiff problems.

using OrdinaryDiffEq

solV = solve(prob, Vern9(), abstol=1e-20); #solve `prob` with the `Vern9` method
retcode: Success
Interpolation: specialized 9th order lazy interpolation
t: 8204-element Array{Float64,1}:
    0.0
    1.0e-6
    6.221690001226793e-6
    3.2087951572557424e-5
    0.0001686533590778439
    0.0008624119118352291
    0.004356627833659284
    0.024714826801768876
    0.1356423149353504
    0.3957733870751676
    ⋮
  999.4261512009862
  999.6191020393967
  999.71466447773
  999.7819314238017
  999.8217522388055
  999.8554129476396
  999.8958265527937
  999.9451341818717
 1000.0
u: 8204-element Array{Array{Float64,1},1}:
 [-0.8, 0.0, 0.0, -0.627641]
 [-0.8, 1.72359e-7, 6.04267e-7, -0.627641]
 [-0.8, 1.07236e-6, 3.75956e-6, -0.627641]
 [-0.8, 5.53065e-6, 1.93897e-5, -0.627641]
 [-0.8, 2.90689e-5, 0.000101912, -0.627641]
 [-0.8, 0.000148644, 0.000521127, -0.627642]
 [-0.799993, 0.000750877, 0.00263263, -0.627652]
 [-0.799763, 0.00425491, 0.0149456, -0.628]
 [-0.792821, 0.0225695, 0.0837728, -0.63833]
 [-0.737015, 0.0483953, 0.279717, -0.714024]
 ⋮
 [0.374578, -0.211395, 0.365016, 0.878803]
 [0.331974, -0.0663988, -0.620518, 1.36192]
 [0.225905, 0.0448727, -1.73933, 1.4738]
 [0.0584876, 0.108359, -3.58527, 0.125023]
 [-0.0775195, 0.049754, -2.56822, -3.03248]
 [-0.122762, -0.0561194, -0.37664, -3.11734]
 [-0.119589, -0.160743, 0.569221, -2.32457]
 [-0.0917824, -0.251539, 0.906529, -1.62064]
 [-0.0545129, -0.321111, 1.00516, -1.09382]

We plot in the $x-y$ synodic plane the solution obtained with TaylorIntegration.jl:

plot(solT, vars=(1, 2), linewidth=1)
scatter!([μ, -1+μ], [0,0], leg=false) # positions of the primaries
xlims!(-1+μ-0.2, 1+μ+0.2)
ylims!(-0.8, 0.8)
xlabel!("x")
ylabel!("y")
-1.0 -0.5 0.0 0.5 1.0 -0.6 -0.3 0.0 0.3 0.6 x y

Note that the orbit obtained displays the expected dynamics: the test particle explores the regions surrounding both primaries, located at the red dots, without escaping to infinity. For comparison, we now plot the orbit corresponding to the solution obtained with the Vern9() integration; note that the scales are identical.

plot(solV, vars=(1, 2), linewidth=1)
scatter!([μ, -1+μ], [0,0], leg=false) # positions of the primaries
xlims!(-1+μ-0.2, 1+μ+0.2)
ylims!(-0.8, 0.8)
xlabel!("x")
ylabel!("y")
-1.0 -0.5 0.0 0.5 1.0 -0.6 -0.3 0.0 0.3 0.6 x y

We note that the orbits do not display the same qualitative features. In particular, the Vern9() integration displays an orbit which does not visit the secondary, as it was the case in the integration using Taylor's method, and stays far enough from $m_1$. The question is which integration should we trust?

We can obtain a quantitative comparison of the validity of both integrations through the preservation of the Jacobi constant:

ET = H.(solT.u)
EV = H.(solV.u)
δET = ET .- J0
δEV = EV .- J0

We plot first the value of the Jacobi constant as function of time.

plot(solT.t, H.(solT.u), label="TaylorIntegration.jl")
plot!(solV.t, H.(solV.u), label="Vern9()")
xlabel!("t")
ylabel!("H")
0 250 500 750 1000 -2.2 -2.0 -1.8 -1.6 t H TaylorIntegration.jl Vern9()

Clearly, the integration with Vern9() does not conserve the Jacobi constant; actually, the fact that its value is strongly reduced leads to the artificial trapping displayed above around $m_1$. We notice that the loss of conservation of the Jacobi constant is actually not related to a close approach with $m_1$.

We now plot, in log scale, the abs of the absolute error in the Jacobi constant as a function of time, for both solutions:

plot(solT.t, abs.(δET), yscale=:log10, label="TaylorIntegration.jl")
plot!(solV.t, abs.(δEV), label="Vern9()")
ylims!(10^-18, 10^4)
xlabel!("t")
ylabel!("dE")
0 250 500 750 1000 10 - 15 10 - 10 10 - 5 10 0 t dE TaylorIntegration.jl Vern9()

We notice that the Jacobi constant absolute error for the TaylorIntegration.jl solution remains bounded below $5\times 10^{-14}$, despite of the fact that the solution displays many close approaches with $m_2$.

Finally, we comment on the time spent by each integration.

@elapsed solve(prob, TaylorMethod(25), abstol=1e-20);
4.44343092
@elapsed solve(prob, Vern9(), abstol=1e-20);
0.083511567

The integration with TaylorMethod() takes much longer than that using Vern9(). Yet, as shown above, the former preserves the Jacobi constant to a high accuracy, whereas the latter solution loses accuracy in the sense of not conserving the Jacobi constant, which is an important property to trust the result of the integration. A fairer comparison is obtained by pushing the native methods of DiffEqs to reach similar accuracy for the integral of motion, as the one obtained by TaylorIntegration.jl. Such comparable situation has a performance cost, which then makes TaylorIntegration.jl comparable or even faster in some cases; see [2].

Finally, as mentioned above, a way to improve the integration time in TaylorIntegration is using the macro @taylorize; see this section for details. Under certain circumstances it is possible to improve the performance, also with the common interface with DifferentialEquations, which restricts some of the great flexibility that DifferentialEquations allows when writing the function containing the differential equations.

References

[1] Murray, Carl D., Stanley F. Dermott. Solar System dynamics. Cambridge University Press, 1999.

[2] DiffEqBenchmarks.jl/DynamicalODE