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 OrdinaryDiffEq.jl together with TaylorIntegration.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.

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

\[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), \tag{1}\]

where

\[V(x, y) = - \frac{1-\mu}{\sqrt{(x-\mu)^2+y^2}} - \frac{\mu}{\sqrt{(x+1-\mu)^2+y^2}}.\tag{2}\]

is the gravitational potential associated to the primaries. The RHS of Eq. (1) 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{aligned} \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{aligned}\]

We define this system of ODEs in a way that allows the use of the @taylorize macro from TaylorIntegration.jl, which for the present example allows important speed-ups. For more details about the specifics of the use of @taylorize, see this section.

using TaylorIntegration
@taylorize function pcr3bp!(dq, q, param, t)
    local μ = param[1]
    local onemμ = 1 - μ
    x1 = q[1]-μ
    x1sq = x1^2
    y = q[2]
    ysq = y^2
    r1_1p5 = (x1sq+ysq)^1.5
    x2 = q[1]+onemμ
    x2sq = x2^2
    r2_1p5 = (x2sq+ysq)^1.5
    dq[1] = q[3] + q[2]
    dq[2] = q[4] - q[1]
    dq[3] = (-((onemμ*x1)/r1_1p5) - ((μ*x2)/r2_1p5)) + q[4]
    dq[4] = (-((onemμ*y )/r1_1p5) - ((μ*y )/r2_1p5)) - q[3]
    return nothing
end

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, fmt = :png)
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)")

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 Vector{Float64}:
 -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 can be used via its common interface bindings with OrdinaryDiffEq.jl; both packages need to be loaded explicitly.

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

using OrdinaryDiffEq
prob = ODEProblem(pcr3bp!, q0, tspan, p)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 2000.0)
u0: 4-element Vector{Float64}:
 -0.8
  0.0
  0.0
 -0.6276410653920694

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

solT = solve(prob, TaylorMethod(25), abstol=1e-15);

retcode: Success
Interpolation: 3rd order Hermite
t: 15343-element Vector{Float64}:
    0.0
    0.13196071872716866
    0.27014676027774454
    0.42393639891281987
    0.6040489984083115
    0.7470888302887118
    0.8537154042859872
    0.9339017026636446
    0.9948376044666346
    1.041388851851968
    ⋮
 1998.110314252107
 1998.2589546650302
 1998.4384851263433
 1998.6592370959424
 1998.9309210910958
 1999.269325515496
 1999.5855925750302
 1999.823076680848
 2000.0
u: 15343-element Vector{Vector{Float64}}:
 [-0.8, 0.0, 0.0, -0.6276410653920694]
 [-0.7932072459999212, 0.021999013383924005, 0.08140695643428739, -0.6377638461316928]
 [-0.7711573074953806, 0.04021587810521527, 0.17675113719300067, -0.6688864600381992]
 [-0.7274049294817609, 0.04873797433924279, 0.30569454759447606, -0.7263176192504541]
 [-0.6471962297961017, 0.0339035917847401, 0.5097162066272879, -0.8242465081730088]
 [-0.555741961444557, -0.004539329245841953, 0.7525448101567436, -0.9233378179280098]
 [-0.4645151896861849, -0.05237933172838659, 1.0331768994048929, -0.99547157889477]
 [-0.37557401376068666, -0.09959841168708136, 1.3584539352908696, -1.016191194976285]
 [-0.2895689539623897, -0.13989162192605029, 1.7261635478522546, -0.9558013178409663]
 [-0.20780076710609213, -0.1693782791151218, 2.1143731596742112, -0.7823572471371022]
 ⋮
 [0.534118578164382, -0.14693555710440903, 0.7859164023751788, 1.1579509194031394]
 [0.6149984984553903, -0.060989727578516506, 0.5212511901406414, 1.1430001114721466]
 [0.6834754715470739, 0.02229888834659935, 0.28920706861122025, 1.084615544384861]
 [0.736803466533942, 0.09573289098244057, 0.08117021523093437, 1.0074366885093562]
 [0.765094033015664, 0.15394471391571177, -0.1217985096106119, 0.9345717267235158]
 [0.7442743989643611, 0.2034319568377962, -0.3636008271512402, 0.8849113345785443]
 [0.6597563679169343, 0.25547057129471346, -0.6408398352694411, 0.8597038019456145]
 [0.54200468670874, 0.31122171045132707, -0.9309778166067714, 0.8099257946089322]
 [0.41204958706401396, 0.3607674400585044, -1.2214470621671558, 0.6920710776913027]

As mentioned above, we will now solve the same problem prob with the Vern9 method from OrdinaryDiffEq, which the DifferentialEquations documentation recommends for high-accuracy (i.e., very low tolerance) integrations of non-stiff problems. Note that, besides setting an absolute tolerance abstol=1e-15, we're setting a relative tolerance reltol=1e-15 [2]. We have found that for the current problem this is a good balance between speed and accuracy for the Vern9 method, i.e., the Vern9 integration becomes noticeably slower (although more accurate) if either abstol or reltol are set to lower values.

using OrdinaryDiffEq

solV = solve(prob, Vern9(), abstol=1e-15, reltol=1e-15); #solve `prob` with the `Vern9` method

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

plot(solT, vars=(1, 2), linewidth=1, fmt = :png)
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")

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, fmt = :png)
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")

We note that both orbits display the same qualitative features, and also some differences. For instance, the TaylorMethod(25) solution gets closer to the primary than that the Vern9(). 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
nothing # hide

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

plot(solT.t, H.(solT.u), label="TaylorMethod(25)", fmt = :png, yformatter = :plain)
plot!(solV.t, H.(solV.u), label="Vern9()")
xlabel!("t")
ylabel!("H")

In the scale shown we observe that, while both solutions display a preservation of the Jacobi constant to a certain degree, the Vern9() solution suffers sudden jumps during the integration.

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="TaylorMethod(25)", legend=:topleft, fmt = :png, yformatter = :plain)
plot!(solV.t, abs.(δEV), label="Vern9()")
ylims!(10^-16, 10^-10)
xlabel!("t")
ylabel!("dE")

We notice that the Jacobi constant absolute error for the TaylorMethod(25) solution remains bounded below $10^{-13}$. The Vern9() solution is, at the end of the integration time, two orders of magnitude above the former for the same quantity.

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

using BenchmarkTools
bT = @benchmark solve($prob, $(TaylorMethod(25)), abstol=1e-15)
bV = @benchmark solve($prob, $(Vern9()), abstol=1e-15, reltol=1e-15)
nothing # hide

bT # TaylorMethod(25) benchmark

BenchmarkTools.Trial: 8 samples with 1 evaluation.
 Range (min … max):  580.989 ms … 800.839 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     658.954 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   671.574 ms ±  68.512 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █          █  █  █       █   █          █                   █  
  █▁▁▁▁▁▁▁▁▁▁█▁▁█▁▁█▁▁▁▁▁▁▁█▁▁▁█▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  581 ms           Histogram: frequency by time          801 ms <

 Memory estimate: 6.33 MiB, allocs estimate: 61572.

bV # Vern9 benchmark

We notice that the TaylorMethod(25) and the Vern9() integrations perform similarly. Yet, as shown above, the former preserves the Jacobi constant to a higher accuracy by two orders of magnitude.

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. For example, for the Vern9() method we can set abstol=1e-20, reltol=1e-15 to get better accuracies. Such situation has a performance cost, which then makes TaylorIntegration comparable or even faster in some cases; see [2].

Finally, as mentioned above, a crucial way in which TaylorIntegration provides high accuracy at competitive speeds is through the use of the @taylorize macro; see this section for details. Currently, TaylorIntegration supports the use of @taylorize via the common interface with DifferentialEquations only for in-place ODEProblem.

References and notes

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

[2] SciMLBenchmarks.jl/DynamicalODE