Hénon–Heiles Energy Conservation

Sebastian Micluța-Câmpeanu, Chris Rackauckas

In this notebook we will study the energy conservation properties of several high-order methods for the Hénon–Heiles system. We will se how the energy error behaves at very thight tolerances and how different techniques such as using symplectic solvers or manifold projections benchmark against each other. The Hamiltonian for this system is given by:

\[ \mathcal{H}=\frac{1}{2}(p_1^2 + p_2^2) + \frac{1}{2}\left(q_1^2 + q_2^2 + 2q_1^2 q_2 - \frac{2}{3}q_2^3\right) \]

We will also compare the in place apporach with the out of place approach by using Arrays (for the in place version) and StaticArrays (for out of place versions). In order to separate these two, we will define the relevant functions and initial conditions in the InPlace and OutofPlace modules. In this way the rest of the code will work for both.

using OrdinaryDiffEq, Plots, DiffEqCallbacks
using TaylorIntegration, LinearAlgebra

T(p) = 1//2 * norm(p)^2
V(q) = 1//2 * (q[1]^2 + q[2]^2 + 2q[1]^2 * q[2]- 2//3 * q[2]^3)
H(p,q, params) = T(p) + V(q)

module InPlace
    using ParameterizedFunctions

    function (dq,p,q,params,t)
        dq[1] = p[1]
        dq[2] = p[2]
    end

    function (dp,p,q,params,t)
        dp[1] = -q[1] * (1 + 2q[2])
        dp[2] = -q[2] - (q[1]^2 - q[2]^2)
    end

    const q0 = [0.1, 0.]
    const p0 = [0., 0.5]
    const u0 = vcat(p0,q0)

    henon = @ode_def HamiltonEqs begin
        dp1 = -q1 * (1 + 2q2)
        dp2 = -q2 - (q1^2 - q2^2)
        dq1 = p1
        dq2 = p2
    end

end

module OutOfPlace
    using StaticArrays

    function (p, q, params, t)
        p
    end

    function (p, q, params, t)
        dp1 = -q[1] * (1 + 2q[2])
        dp2 = -q[2] - (q[1]^2 - q[2]^2)
        @SVector [dp1, dp2]
    end

    const q0 = @SVector [0.1, 0.]
    const p0 = @SVector [0., 0.5]
    const u0 = vcat(p0,q0)

    henon(z, p, t) = SVector(
        -z[3] * (1 + 2z[4]),
        -z[4] - (z[3]^2 - z[4]^2),
        z[1],
        z[2]
    )

end

function g(resid, u)
    resid[1] = H([u[1],u[2]], [u[3],u[4]], nothing) - E
    resid[2:4] .= 0
end

const cb = ManifoldProjection(g, nlopts=Dict(:ftol=>1e-13))

const E = H(InPlace.p0, InPlace.q0, nothing)
0.13

For the comparison we will use the following function

energy_err(sol) = map(i->H([sol[1,i], sol[2,i]], [sol[3,i], sol[4,i]], nothing)-E, 1:length(sol.u))
abs_energy_err(sol) = [abs.(H([sol[1,j], sol[2,j]], [sol[3,j], sol[4,j]], nothing) - E) for j=1:length(sol.u)]

function compare(mode=InPlace, all=true, plt=nothing; tmax=1e2)
    prob1 = DynamicalODEProblem(mode., mode., mode.p0, mode.q0, (0., tmax))
    prob2 = ODEProblem(mode.henon, mode.u0, (0., tmax))

    GC.gc()
    (mode == InPlace  && all) && @time sol1 = solve(prob2, Vern9(), callback=cb, abstol=1e-14, reltol=1e-14)
    GC.gc()
    @time sol2 = solve(prob1, KahanLi8(), dt=1e-2, maxiters=1e10)
    GC.gc()
    @time sol3 = solve(prob1, SofSpa10(), dt=1e-2, maxiters=1e8)
    GC.gc()
    @time sol4 = solve(prob2, Vern9(), abstol=1e-14, reltol=1e-14)
    GC.gc()
    @time sol5 = solve(prob1, DPRKN12(), abstol=1e-14, reltol=1e-14)
    GC.gc()
    (mode == InPlace && all) && @time sol6 = solve(prob2, TaylorMethod(50), abstol=1e-20)

    (mode == InPlace && all) && println("Vern9 + ManifoldProjection max energy error:\t"*
        "$(maximum(abs_energy_err(sol1)))\tin\t$(length(sol1.u))\tsteps.")
    println("KahanLi8 max energy error:\t\t\t$(maximum(abs_energy_err(sol2)))\tin\t$(length(sol2.u))\tsteps.")
    println("SofSpa10 max energy error:\t\t\t$(maximum(abs_energy_err(sol3)))\tin\t$(length(sol3.u))\tsteps.")
    println("Vern9 max energy error:\t\t\t\t$(maximum(abs_energy_err(sol4)))\tin\t$(length(sol4.u))\tsteps.")
    println("DPRKN12 max energy error:\t\t\t$(maximum(abs_energy_err(sol5)))\tin\t$(length(sol5.u))\tsteps.")
    (mode == InPlace && all) && println("TaylorMethod max energy error:\t\t\t$(maximum(abs_energy_err(sol6)))"*
        "\tin\t$(length(sol6.u))\tsteps.")

    if plt == nothing
        plt = plot(xlabel="t", ylabel="Energy error")
    end
    (mode == InPlace && all) && plot!(sol1.t, energy_err(sol1), label="Vern9 + ManifoldProjection")
    plot!(sol2.t, energy_err(sol2), label="KahanLi8", ls=mode==InPlace ? :solid : :dash)
    plot!(sol3.t, energy_err(sol3), label="SofSpa10", ls=mode==InPlace ? :solid : :dash)
    plot!(sol4.t, energy_err(sol4), label="Vern9", ls=mode==InPlace ? :solid : :dash)
    plot!(sol5.t, energy_err(sol5), label="DPRKN12", ls=mode==InPlace ? :solid : :dash)
    (mode == InPlace && all) && plot!(sol6.t, energy_err(sol6), label="TaylorMethod")

    return plt
end
compare (generic function with 4 methods)

The mode argument choses between the in place approach and the out of place one. The all parameter is used to compare only the integrators that support both the in place and the out of place versions (we reffer here only to the 6 high order methods chosen bellow). The plt argument can be used to overlay the results over a previous plot and the tmax keyword determines the simulation time.

Note:

  1. The Vern9 method is used with ODEProblem because of performance issues with ArrayPartition indexing which manifest for DynamicalODEProblem.

  2. The NLsolve call used by ManifoldProjection was modified to use ftol=1e-13 in order to obtain a very low energy error.

Here are the results of the comparisons between the in place methods:

compare(tmax=1e2)
18.761184 seconds (51.98 M allocations: 2.391 GiB, 10.69% gc time)
  4.001114 seconds (11.54 M allocations: 548.627 MiB, 7.63% gc time)
  2.484655 seconds (6.23 M allocations: 278.580 MiB, 5.75% gc time)
  2.367251 seconds (4.42 M allocations: 188.205 MiB, 4.56% gc time)
  2.546640 seconds (5.70 M allocations: 285.209 MiB, 7.11% gc time)
  1.487699 seconds (3.28 M allocations: 212.654 MiB, 6.23% gc time)
Vern9 + ManifoldProjection max energy error:	2.7755575615628914e-16	in	1881
	steps.
KahanLi8 max energy error:			4.9404924595819466e-15	in	10001	steps.
SofSpa10 max energy error:			5.440092820663267e-15	in	10001	steps.
Vern9 max energy error:				2.7755575615628914e-16	in	941	steps.
DPRKN12 max energy error:			2.220446049250313e-16	in	334	steps.
TaylorMethod max energy error:			1.942890293094024e-16	in	80	steps.
compare(tmax=1e3)
0.109845 seconds (625.18 k allocations: 52.859 MiB, 21.15% gc time)
  0.198922 seconds (1.50 M allocations: 87.745 MiB, 56.57% gc time)
  0.294718 seconds (1.50 M allocations: 87.746 MiB, 60.47% gc time)
  0.015506 seconds (121.39 k allocations: 13.125 MiB)
  0.007512 seconds (79.48 k allocations: 3.551 MiB)
  0.146455 seconds (1.57 M allocations: 245.919 MiB, 15.10% gc time)
Vern9 + ManifoldProjection max energy error:	5.440092820663267e-15	in	18659
	steps.
KahanLi8 max energy error:			1.815214645262131e-14	in	100002	steps.
SofSpa10 max energy error:			2.8033131371785203e-14	in	100002	steps.
Vern9 max energy error:				5.440092820663267e-15	in	9330	steps.
DPRKN12 max energy error:			1.582067810090848e-15	in	3245	steps.
TaylorMethod max energy error:			5.551115123125783e-16	in	784	steps.
compare(tmax=1e4)
1.526229 seconds (6.25 M allocations: 525.696 MiB, 38.92% gc time)
  2.232151 seconds (15.00 M allocations: 831.610 MiB, 60.18% gc time)
  4.546507 seconds (15.00 M allocations: 831.611 MiB, 73.96% gc time)
  2.119072 seconds (1.21 M allocations: 129.866 MiB, 93.13% gc time)
  0.372237 seconds (791.58 k allocations: 34.745 MiB, 81.60% gc time)
  3.290178 seconds (15.73 M allocations: 2.065 GiB, 60.74% gc time)
Vern9 + ManifoldProjection max energy error:	4.1300296516055823e-14	in	1864
29	steps.
KahanLi8 max energy error:			3.161360062620133e-14	in	1000001	steps.
SofSpa10 max energy error:			1.136590821460004e-13	in	1000001	steps.
Vern9 max energy error:				4.1300296516055823e-14	in	93215	steps.
DPRKN12 max energy error:			1.2045919817182948e-14	in	32364	steps.
TaylorMethod max energy error:			2.7200464103316335e-15	in	7826	steps.
compare(tmax=5e4)
6.490849 seconds (31.51 M allocations: 2.576 GiB, 20.89% gc time)
 11.779758 seconds (75.00 M allocations: 4.061 GiB, 61.56% gc time)
 16.448389 seconds (75.00 M allocations: 4.061 GiB, 63.96% gc time)
  0.758374 seconds (6.06 M allocations: 649.249 MiB)
  8.827864 seconds (3.96 M allocations: 171.275 MiB, 95.95% gc time)
 19.844553 seconds (78.64 M allocations: 10.175 GiB, 66.39% gc time)
Vern9 + ManifoldProjection max energy error:	1.0000333894311098e-13	in	9320
97	steps.
KahanLi8 max energy error:			1.2331802246023926e-13	in	5000001	steps.
SofSpa10 max energy error:			1.5035195310986182e-13	in	5000001	steps.
Vern9 max energy error:				2.114697306154767e-13	in	466050	steps.
DPRKN12 max energy error:			6.164513344231182e-14	in	161766	steps.
TaylorMethod max energy error:			8.770761894538737e-15	in	39124	steps.

We can see that as the simulation time increases, the energy error increases. For this particular example the energy error for all the methods is comparable. For relatively short simulation times, if a highly accurate solution is required, the symplectic method is not recommended as its energy error fluctuations are larger than for other methods. An other thing to notice is the fact that the two versions of Vern9 behave identically, as expected, untill the energy error set by ftol is reached.

We will now compare the in place with the out of place versions. In the plots bellow we will use a dashed line for the out of place versions.

function in_vs_out(;all=false, tmax=1e2)
    println("In place versions:")
    plt = compare(InPlace, all, tmax=tmax)
    println("\nOut of place versions:")
    plt = compare(OutOfPlace, false, plt; tmax=tmax)
end
in_vs_out (generic function with 1 method)

First, here is a summary of all the available methods for tmax = 1e3:

in_vs_out(all=true, tmax=1e3)
In place versions:
  0.106317 seconds (625.18 k allocations: 52.859 MiB, 17.85% gc time)
  0.243425 seconds (1.50 M allocations: 87.745 MiB, 63.65% gc time)
  0.256927 seconds (1.50 M allocations: 87.746 MiB, 49.43% gc time)
  0.018548 seconds (121.39 k allocations: 13.125 MiB)
  0.009028 seconds (79.48 k allocations: 3.551 MiB)
  0.153892 seconds (1.57 M allocations: 245.919 MiB, 13.76% gc time)
Vern9 + ManifoldProjection max energy error:	5.440092820663267e-15	in	18659
	steps.
KahanLi8 max energy error:			1.815214645262131e-14	in	100002	steps.
SofSpa10 max energy error:			2.8033131371785203e-14	in	100002	steps.
Vern9 max energy error:				5.440092820663267e-15	in	9330	steps.
DPRKN12 max energy error:			1.582067810090848e-15	in	3245	steps.
TaylorMethod max energy error:			5.551115123125783e-16	in	784	steps.

Out of place versions:
  3.695491 seconds (8.18 M allocations: 432.014 MiB, 6.74% gc time)
  1.594871 seconds (3.75 M allocations: 199.755 MiB, 6.12% gc time)
  5.157518 seconds (9.11 M allocations: 457.788 MiB, 5.06% gc time)
  2.080734 seconds (4.42 M allocations: 218.469 MiB, 5.13% gc time)
KahanLi8 max energy error:			1.2684298056342413e-14	in	100002	steps.
SofSpa10 max energy error:			1.9345636204093353e-14	in	100002	steps.
Vern9 max energy error:				3.635980405647388e-15	in	9330	steps.
DPRKN12 max energy error:			1.3600232051658168e-15	in	3238	steps.

Now we will compare the in place and the out of place versions, but only for the integrators that are compatible with StaticArrays

in_vs_out(tmax=1e2)
In place versions:
  0.008173 seconds (150.11 k allocations: 8.322 MiB)
  0.010488 seconds (150.11 k allocations: 8.323 MiB)
  0.001433 seconds (12.33 k allocations: 1.348 MiB)
  0.001311 seconds (8.29 k allocations: 370.516 KiB)
KahanLi8 max energy error:			4.9404924595819466e-15	in	10001	steps.
SofSpa10 max energy error:			5.440092820663267e-15	in	10001	steps.
Vern9 max energy error:				2.7755575615628914e-16	in	941	steps.
DPRKN12 max energy error:			2.220446049250313e-16	in	334	steps.

Out of place versions:
  0.003331 seconds (20.07 k allocations: 1.989 MiB)
  0.005749 seconds (20.07 k allocations: 1.990 MiB)
  0.000841 seconds (1.97 k allocations: 547.156 KiB)
  0.000646 seconds (1.77 k allocations: 129.453 KiB)
KahanLi8 max energy error:			2.7478019859472624e-15	in	10001	steps.
SofSpa10 max energy error:			4.6074255521944e-15	in	10001	steps.
Vern9 max energy error:				2.7755575615628914e-16	in	940	steps.
DPRKN12 max energy error:			1.6653345369377348e-16	in	329	steps.
in_vs_out(tmax=1e3)
In place versions:
  0.218601 seconds (1.50 M allocations: 87.745 MiB, 52.16% gc time)
  0.315477 seconds (1.50 M allocations: 87.746 MiB, 57.17% gc time)
  0.014030 seconds (121.39 k allocations: 13.125 MiB)
  0.008943 seconds (79.48 k allocations: 3.551 MiB)
KahanLi8 max energy error:			1.815214645262131e-14	in	100002	steps.
SofSpa10 max energy error:			2.8033131371785203e-14	in	100002	steps.
Vern9 max energy error:				5.440092820663267e-15	in	9330	steps.
DPRKN12 max energy error:			1.582067810090848e-15	in	3245	steps.

Out of place versions:
  0.044473 seconds (200.08 k allocations: 25.946 MiB, 11.32% gc time)
  0.064648 seconds (200.08 k allocations: 25.947 MiB, 10.99% gc time)
  0.008354 seconds (18.76 k allocations: 4.889 MiB)
  0.005371 seconds (16.70 k allocations: 1.398 MiB)
KahanLi8 max energy error:			1.2684298056342413e-14	in	100002	steps.
SofSpa10 max energy error:			1.9345636204093353e-14	in	100002	steps.
Vern9 max energy error:				3.635980405647388e-15	in	9330	steps.
DPRKN12 max energy error:			1.3600232051658168e-15	in	3238	steps.
in_vs_out(tmax=1e4)
In place versions:
  2.109015 seconds (15.00 M allocations: 831.610 MiB, 49.93% gc time)
  4.400183 seconds (15.00 M allocations: 831.611 MiB, 69.05% gc time)
  2.000854 seconds (1.21 M allocations: 129.866 MiB, 91.33% gc time)
  0.366431 seconds (791.58 k allocations: 34.745 MiB, 76.15% gc time)
KahanLi8 max energy error:			3.161360062620133e-14	in	1000001	steps.
SofSpa10 max energy error:			1.136590821460004e-13	in	1000001	steps.
Vern9 max energy error:				4.1300296516055823e-14	in	93215	steps.
DPRKN12 max energy error:			1.2045919817182948e-14	in	32364	steps.

Out of place versions:
  0.609828 seconds (2.00 M allocations: 198.370 MiB, 37.44% gc time)
  1.091297 seconds (2.00 M allocations: 198.371 MiB, 49.46% gc time)
  0.107677 seconds (186.54 k allocations: 44.810 MiB, 24.70% gc time)
  0.051721 seconds (166.17 k allocations: 12.727 MiB)
KahanLi8 max energy error:			5.5261351050717167e-14	in	1000001	steps.
SofSpa10 max energy error:			2.3342439092743916e-14	in	1000001	steps.
Vern9 max energy error:				4.343747583845925e-14	in	93216	steps.
DPRKN12 max energy error:			1.2323475573339238e-14	in	32350	steps.
in_vs_out(tmax=5e4)
In place versions:
  9.043584 seconds (75.00 M allocations: 4.061 GiB, 41.65% gc time)
 16.132004 seconds (75.00 M allocations: 4.061 GiB, 57.78% gc time)
  0.876920 seconds (6.06 M allocations: 649.249 MiB)
  8.384477 seconds (3.96 M allocations: 171.275 MiB, 94.14% gc time)
KahanLi8 max energy error:			1.2331802246023926e-13	in	5000001	steps.
SofSpa10 max energy error:			1.5035195310986182e-13	in	5000001	steps.
Vern9 max energy error:				2.114697306154767e-13	in	466050	steps.
DPRKN12 max energy error:			6.164513344231182e-14	in	161766	steps.

Out of place versions:
  2.595858 seconds (10.00 M allocations: 991.827 MiB, 26.12% gc time)
  4.268870 seconds (10.00 M allocations: 991.828 MiB, 35.35% gc time)
  1.949432 seconds (932.22 k allocations: 227.099 MiB, 78.97% gc time)
  0.281741 seconds (830.48 k allocations: 54.221 MiB, 8.69% gc time)
KahanLi8 max energy error:			9.178768856088482e-14	in	5000001	steps.
SofSpa10 max energy error:			1.1290968160437842e-13	in	5000001	steps.
Vern9 max energy error:				2.2346013928142838e-13	in	466051	steps.
DPRKN12 max energy error:			6.533662499919046e-14	in	161742	steps.

As we see from the above comparisons, the StaticArray versions are significantly faster and use less memory. The speedup provided for the out of place version is more proeminent at larger values for tmax. We can see again that if the simulation time is increased, the energy error of the symplectic methods is less noticeable compared to the rest of the methods.

The benchmarks were performed on a machine with

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

Appendix

These benchmarks are a part of the DiffEqBenchmarks.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqBenchmarks.jl

To locally run this tutorial, do the following commands:

using DiffEqBenchmarks
DiffEqBenchmarks.weave_file("DynamicalODE","Henon-Heiles_energy_conservation_benchmark.jmd")

Computer Information:

Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, haswell)

Package Information:

Status: `/home/crackauckas/.julia/environments/v1.1/Project.toml`
[c52e3926-4ff0-5f6e-af25-54175e0327b1] Atom 0.8.5
[bcd4f6db-9728-5f36-b5f7-82caef46ccdb] DelayDiffEq 5.2.0
[bb2cbb15-79fc-5d1e-9bf1-8ae49c7c1650] DiffEqBenchmarks 0.1.0
[459566f4-90b8-5000-8ac3-15dfb0a30def] DiffEqCallbacks 2.5.2
[f3b72e0c-5b89-59e1-b016-84e28bfd966d] DiffEqDevTools 2.7.2
[78ddff82-25fc-5f2b-89aa-309469cbf16f] DiffEqMonteCarlo 0.14.0
[77a26b50-5914-5dd7-bc55-306e6241c503] DiffEqNoiseProcess 3.1.0
[055956cb-9e8b-5191-98cc-73ae4a59e68a] DiffEqPhysics 3.1.0
[a077e3f3-b75c-5d7f-a0c6-6bc4c8ec64a9] DiffEqProblemLibrary 4.1.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.3.0
[b305315f-e792-5b7a-8f41-49f472929428] Elliptic 0.5.0
[e5e0dc1b-0480-54bc-9374-aad01c23163d] Juno 0.7.0
[7f56f5a3-f504-529b-bc02-0b1fe5e64312] LSODA 0.4.0
[c030b06c-0b6d-57c2-b091-7029874bd033] ODE 2.4.0
[54ca160b-1b9f-5127-a996-1867f4bc2a2c] ODEInterface 0.4.5
[09606e27-ecf5-54fc-bb29-004bd9f985bf] ODEInterfaceDiffEq 3.1.0
[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.6.0
[2dcacdae-9679-587a-88bb-8b444fb7085b] ParallelDataTransfer 0.5.0
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 4.1.1
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 0.24.0
[d330b81b-6aea-500a-939a-2ce795aea3ee] PyPlot 2.8.1
[90137ffa-7385-5640-81b9-e52037218182] StaticArrays 0.10.3
[789caeaf-c7a9-5a7d-9973-96adeb23e2a0] StochasticDiffEq 6.1.1
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 3.4.0
[92b13dbe-c966-51a2-8445-caca9f8a7d42] TaylorIntegration 0.4.1
[44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9] Weave 0.9.0