33Here, we show an example of interoperability between ` TaylorIntegration.jl ` and
44[ ` DifferentialEquations.jl ` ] ( https://github.com/JuliaDiffEq/DifferentialEquations.jl ) , i.e.,
55how to use ` TaylorIntegration.jl ` from the ` DifferentialEquations `
6- ecosystem. the basic requirement is to load ` DiffEqBase.jl ` , which
6+ ecosystem. The basic requirement is to load ` DiffEqBase.jl ` , which
77sets-up the common interface.
8- Below, we shall also use ` ParameterizedFunctions.jl ` to define the appropriate
9- system of ODEs, and use ` OrdinaryDiffEq.jl ` to compare
8+ Below, we shall also use ` OrdinaryDiffEq.jl ` to compare
109the accuracy of ` TaylorIntegration.jl ` with respect to
1110high-accuracy methods for non-stiff problems (` Vern9 ` method).
11+ While ` DifferentialEquations ` offers many macros to simplify certain
12+ aspects, we do not rely on them simply because using properly ` @taylorize `
13+ improves the performance.
1214
1315The problem we will integrate in this example is the planar circular restricted
1416three-body problem (PCR3BP, also capitalized as PCRTBP). The PCR3BP describes
1517the motion of a body with negligible mass `` m_3 `` under the gravitational influence of two
1618bodies with masses `` m_1 `` and `` m_2 `` , such that `` m_1 \ge m_2 `` . It is
17- assumed that `` m_3 `` is much smaller than the other two masses, and therefore it
19+ assumed that `` m_3 `` is much smaller than the other two masses so it
20+ does not influence their motion, and therefore it
1821is simply considered as a massless test particle.
1922The body with the greater mass `` m_1 `` is referred as the * primary* , and
2023`` m_2 `` as the * secondary* . These bodies are together
@@ -33,6 +36,7 @@ using Plots
3336const μ = 0.01
3437nothing # hide
3538```
39+
3640The Hamiltonian for the PCR3BP in the synodic frame (i.e., a frame which rotates
3741such that the primaries are at rest on the `` x `` axis) is
3842``` math
@@ -51,7 +55,7 @@ V(x, y) = - \frac{1-\mu}{\sqrt{(x-\mu)^2+y^2}} - \frac{\mu}{\sqrt{(x+1-\mu)^2+y^
5155is the gravitational potential associated to the primaries. The RHS of Eq.
5256(\ref{eq-pcr3bp-hamiltonian}) is also known as the * Jacobi constant* , since it is a
5357preserved quantity of motion in the PCR3BP. We will use this property to check
54- the accuracy of the obtained solution .
58+ the accuracy of the solutions computed .
5559``` @example common
5660V(x, y) = - (1-μ)/sqrt((x-μ)^2+y^2) - μ/sqrt((x+1-μ)^2+y^2)
5761H(x, y, px, py) = (px^2+py^2)/2 - (x*py-y*px) + V(x, y)
@@ -69,26 +73,35 @@ The equations of motion for the PCR3BP are
6973 \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.
7074\end{eqnarray}
7175```
72- We define this system of ODEs with ` ParameterizedFunctions.jl `
73- ``` @example common
74- using ParameterizedFunctions
7576
76- f = @ode_def PCR3BP begin
77- dx = px + y
78- dy = py - x
79- dpx = - (1-μ)*(x-μ)*((x-μ)^2+y^2)^-1.5 - μ*(x+1-μ)*((x+1-μ)^2+y^2)^-1.5 + py
80- dpy = - (1-μ)*y *((x-μ)^2+y^2)^-1.5 - μ*y *((x+1-μ)^2+y^2)^-1.5 - px
81- end μ
77+ We define this system of ODEs using the most naive approach
78+ ``` @example common
79+ function f(dq, q, param, t)
80+ local μ = param[1]
81+ x, y, px, py = q
82+ dq[1] = px + y
83+ dq[2] = py - x
84+ dq[3] = - (1-μ)*(x-μ)*((x-μ)^2+y^2)^-1.5 - μ*(x+1-μ)*((x+1-μ)^2+y^2)^-1.5 + py
85+ dq[4] = - (1-μ)*y *((x-μ)^2+y^2)^-1.5 - μ*y *((x+1-μ)^2+y^2)^-1.5 - px
86+ return nothing
87+ end
8288nothing # hide
8389```
90+ Note that ` DifferentialEquations ` offers interesting alternatives to write
91+ these equations of motion in a simpler and more convenient way,
92+ for example, using the macro ` @ode_def ` ,
93+ see [ ` ParameterizedFunctions.jl ` ] ( https://github.com/JuliaDiffEq/ParameterizedFunctions.jl ) .
94+ We have not used that flexibility here because ` TaylorIntegration.jl ` has
95+ [ ` @taylorize ` ] ( @ref ) , which under certain circumstances allows to
96+ important speed-ups.
8497
8598We shall define the initial conditions `` q_0 = (x_0, y_0, p_{x,0}, p_{y,0}) `` such that
8699`` H(q_0) = J_0 `` , where `` J_0 `` is a prescribed value. In order to do this,
87100we select `` y_0 = p_{x,0} = 0 `` and compute the value of `` p_{y,0} `` for which
88101`` H(q_0) = J_0 `` holds.
89102
90103We consider a value for `` J_0 `` such that the test particle is
91- able to display close encounters with both primaries, but cannot escape to infinity.
104+ able to display close encounters with * both* primaries, but cannot escape to infinity.
92105We may obtain a first
93106approximation to the desired value of `` J_0 `` if we plot the projection of the
94107zero-velocity curves on the `` x `` -axis.
@@ -135,8 +148,9 @@ Hamiltonian evaluated at the initial condition is indeed equal to `J0`.
135148H(q0) == J0
136149```
137150
138- Following [ the tutorial] ( http://docs.juliadiffeq.org/latest/tutorials/ode_example.html )
139- of ` DifferentialEquations.jl ` , we define an ` ODEProblem ` for the integration;
151+ Following the ` DifferentialEquations.jl `
152+ [ tutorial] ( http://docs.juliadiffeq.org/latest/tutorials/ode_example.html ) ,
153+ we define an ` ODEProblem ` for the integration;
140154` TaylorIntegration.jl ` can be used via its common interface bindings with ` DiffEqBase.jl ` ; both packages need to be loaded explicitly.
141155``` @example common
142156tspan = (0.0, 1000.0)
@@ -233,18 +247,28 @@ Finally, we comment on the time spent by each integration.
233247``` @example common
234248@elapsed solve(prob, Vern9(), abstol=1e-20);
235249```
236- Clearly, the integration with ` TaylorMethod() ` takes * much longer* than that using
250+ The integration with ` TaylorMethod() ` takes * much longer* than that using
237251` Vern9() ` . Yet, as shown above, the former preserves the Jacobi constant
238- to a high accuracy while displaying the correct dynamics; whereas the latter
252+ to a high accuracy, whereas the latter
239253solution loses accuracy in the sense of not conserving the Jacobi constant,
240254which is an important property to trust the result of the integration.
255+ A fairer comparison is obtained by pushing the native methods of ` DiffEqs `
256+ to reach similar accuracy for the integral of motion, as the one
257+ obtained by ` TaylorIntegration.jl ` . Such comparable situation has
258+ a performance cost, which then makes ` TaylorIntegration.jl ` comparable
259+ or even faster in some cases; see [[ 2]] (@ref refsPCR3BP).
241260
242- Finally, we shall mention here that a way to improve the integration time in
261+ Finally, as mentioned above, a way to improve the integration time in
243262` TaylorIntegration ` is using the macro [ ` @taylorize ` ] ( @ref ) ; see
244- [ this section] (@ref taylorize) for details. Yet, using the macro is not
245- currently compatible with the common interface with ` DifferentialEquations ` .
263+ [ this section] (@ref taylorize) for details. Under certain circumstances
264+ it is possible to improve the performance, also with the common interface
265+ with ` DifferentialEquations ` , which restricts some of the great
266+ flexibility that ` DifferentialEquations ` allows when writing the function
267+ containing the differential equations.
246268
247269
248270### [ References] (@id refsPCR3BP)
249271
250272[ 1] Murray, Carl D., Stanley F. Dermott. Solar System dynamics. Cambridge University Press, 1999.
273+
274+ [ 2] [ DiffEqBenchmarks.jl/DynamicalODE] ( https://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqBenchmarks.jl/blob/master/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.ipynb )
0 commit comments