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: Taylor series polynomial evaluation
t: 19608-element Vector{Float64}:
    0.0
    0.13196071872716866
    0.27014676027774454
    0.42393639891281987
    0.6040489984083115
    0.7470888302887118
    0.8537154042859872
    0.9339017026636446
    0.9948376044666346
    1.041388851851968
    ⋮
 1999.5042725273079
 1999.5431825892883
 1999.585321353172
 1999.6333069204968
 1999.6900342956144
 1999.7583339712685
 1999.841241219841
 1999.9430374708056
 2000.0
u: 19608-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.17496376257492635, -0.17162024831531136, 1.3853825963776936, -1.9910183910679593]
 [-0.12304118303795919, -0.23591024736030314, 1.664519586011635, -1.6072083780763986]
 [-0.06079568084101738, -0.2910297290272401, 1.7937412752512865, -1.201602071416007]
 [0.010706160716735615, -0.33799420974483546, 1.8000388400926512, -0.8222349820550685]
 [0.09036540182670252, -0.3776915547287994, 1.7168539805662775, -0.4948839982771776]
 [0.1759176142431495, -0.4109734546288778, 1.5764514213719614, -0.2305013931376877]
 [0.2641294989793357, -0.4393446214896214, 1.4058003373126944, -0.02725826011770835]
 [0.3515484829667911, -0.46520091680777775, 1.2228269820636821, 0.12590172434120908]
 [0.39178452468014624, -0.47746330663809045, 1.1343053892122903, 0.1855857869403744]

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

retcode: Success
Interpolation: specialized 9th order lazy interpolation
t: 109969-element Vector{Float64}:
    0.0
    0.014218880460979355
    0.029254985035048104
    0.043100405390433874
    0.05734634963317922
    0.07358281996048438
    0.08944879825305627
    0.10587664107366627
    0.1226246832421739
    0.1435160589250072
    ⋮
 1999.8835158404124
 1999.8979213519083
 1999.9116672471735
 1999.9261589673808
 1999.941671239634
 1999.9563042419318
 1999.9717788474554
 1999.987747535878
 2000.0
u: 109969-element Vector{Vector{Float64}}:
 [-0.8, 0.0, 0.0, -0.6276410653920694]
 [-0.7999214879562981, 0.0024498152778736912, 0.00859413588182428, -0.6277599916692193]
 [-0.7996675843508191, 0.005034208510608978, 0.017696409719266754, -0.6281442610049114]
 [-0.7992782930574767, 0.007402685037193327, 0.026103506293378813, -0.6287324566601369]
 [-0.7987218689531319, 0.00982280232842162, 0.034792021624444995, -0.6295711582485729]
 [-0.797894470875534, 0.012553097170943344, 0.04475750864054082, -0.63081399178736]
 [-0.7968864732598321, 0.015184680375190484, 0.0545768793929235, -0.6323213310978815]
 [-0.7956341732152393, 0.017863254033671463, 0.06484608541850159, -0.634183992679585]
 [-0.7941379505860766, 0.02053692742751024, 0.07544018310964262, -0.636395532591486]
 [-0.7919590653442087, 0.023777693641981894, 0.08885831570084259, -0.6395904165292945]
 ⋮
 [-0.10036556390021033, -0.40881948168063087, 1.3263848652497512, -1.0339244230210929]
 [-0.08720892823246354, -0.4219534753803598, 1.3306307677028557, -0.9775508905357453]
 [-0.07478774834940205, -0.43392235178080524, 1.332218341540104, -0.9262738327498371]
 [-0.0618595087150646, -0.4459796381322718, 1.33164440577699, -0.8748026893703007]
 [-0.048236112287607984, -0.4582871679437302, 1.328841777142106, -0.8225447583250936]
 [-0.03560926760707582, -0.46936547175800486, 1.3244363154112246, -0.7758223093021686]
 [-0.022510911045455013, -0.480556024781558, 1.318210211953628, -0.7289957477163415]
 [-0.009284187203683182, -0.4915754146264585, 1.3103724961844636, -0.6833061494007471]
 [0.0006569134087038978, -0.49968918682462987, 1.303542664918598, -0.6499623338440916]

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

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}$ throughout the integration. While the Vern9() solution at the end of the integration time has reached a similar value, it displays a larger Jacobi constant error earlier in time.

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)

bT # TaylorMethod(25) benchmark

BenchmarkTools.Trial: 50 samples with 1 evaluation.
 Range (min … max):   96.911 ms … 145.101 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     100.137 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   101.179 ms ±   7.023 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃   █                                                          
  █▄▃▇█▄▁▃▁▄▃▁▃▃▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃ ▁
  96.9 ms          Histogram: frequency by time          145 ms <

 Memory estimate: 8.66 MiB, allocs estimate: 117851.

bV # Vern9 benchmark

BenchmarkTools.Trial: 20 samples with 1 evaluation.
 Range (min … max):  214.008 ms … 679.944 ms  ┊ GC (min … max):  0.00% … 67.88%
 Time  (median):     235.338 ms               ┊ GC (median):     0.00%
 Time  (mean ± σ):   273.811 ms ± 104.240 ms  ┊ GC (mean ± σ):  16.78% ± 18.28%

    █                                                            
  ▅▁█▃▃▃▃▁▁▁▃▁▁▁▁▃▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃ ▁
  214 ms           Histogram: frequency by time          680 ms <

 Memory estimate: 130.11 MiB, allocs estimate: 1319694.

We notice in this setup, where the TaylorMethod(25) and the Vern9() integrations perform similarly in terms of accuracy, the former performs better in terms of runtime.

We can tune the abstol and reltol for the Vern9() method we so that performance is similar. Such situation has an accuracy cost, which then makes TaylorIntegration a sensible alternative for high-accuracy integrations of non-stiff ODEs 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