Simple gravity pendulum¶

Simple gravity pendulum:

$\frac{d^2\theta}{dt^2} = -\frac{g}{l}\sin\theta.$

Initial condition for $-1 < y_{\mathrm{max}} < 1$

$\theta(0) = 0, \quad \dot\theta(0) = \sqrt{\frac{g}{l}}\;2y_{\mathrm{max}}$

Exact solution:

$\theta(t) = 2\arcsin\left( y_{max}\mathrm{sn}\left(\sqrt{\frac{g}{l}}\;t\;\Bigg|\;k=y_{\mathrm{max}}\right) \right)$

@show VERSION

second(a) = a[begin+1]

@time using Printf
@time using DifferentialEquations
@time using StaticArrays
@time using Parameters
@time using Elliptic

@time using Plots
@time pyplot(fmt=:svg)

@time using Unitful
Base.zero(::Type{Quantity{T}}) where T = zero(T)
const m = u"m"
const s = u"s"
val(x) = x
val(x::Unitful.AbstractQuantity) = x.val

@time using PhysicalConstants
g = PhysicalConstants.CODATA2018.StandardAccelerationOfGravitation

VERSION = v"1.7.0-DEV.1129"
  0.001262 seconds (1.61 k allocations: 117.172 KiB, 135.22% compilation time)
  8.646942 seconds (21.28 M allocations: 1.542 GiB, 4.93% gc time)
  0.021959 seconds (33.37 k allocations: 1.902 MiB, 97.53% compilation time)
  0.000396 seconds (251 allocations: 16.391 KiB)
  0.003836 seconds (3.21 k allocations: 282.844 KiB, 28.24% compilation time)
  3.418940 seconds (6.87 M allocations: 508.058 MiB, 7.96% gc time)
  4.658772 seconds (12.27 M allocations: 735.629 MiB, 3.57% gc time, 0.35% compilation time)
  0.022045 seconds (33.37 k allocations: 1.902 MiB, 97.57% compilation time)
  0.312106 seconds (459.72 k allocations: 30.415 MiB)

WARNING: Wrapping `Vararg` directly in UnionAll is deprecated (wrap the tuple instead).

Out[1]:

Standard acceleration of gravitation (g_n)
Value                         = 9.80665 m s^-2
Standard uncertainty          = (exact)
Relative standard uncertainty = (exact)
Reference                     = CODATA 2018

In [2]:

# parameters with dimensions
g = PhysicalConstants.CODATA2018.StandardAccelerationOfGravitation # = 9.80665 m/s²
l = 1.0m
θ0 = 0.0
ymax = 0.9999 # this setting is very stiff!
dθ0 = √(g/l)*2ymax
t0 = 0.0s
t1 = 30.0s

# dimensionless parameters
p = val.((g, l)) 
tspan = val.((t0, t1))
v0 = val.([dθ0])
u0 = val.([θ0])
vu0 = val.([dθ0, θ0])

(; p, ymax, vu0)

Out[2]:

(p = (9.80665, 1.0), ymax = 0.9999, vu0 = [6.262487929909805, 0.0])

In [3]:

# exact solution for θ0 = 0, dθ(0)/dt = √(g/l)*2y_max
sn(u, k) = Jacobi.sn(u, k^2)
θ_exact_sol(p, t, ymax) = ((g, l) = p; 2asin(ymax*sn(√(g/l)*t, ymax)))

t = range(tspan...; length=1001)
plot(t, θ_exact_sol.(Ref(p), t, ymax); label="exact solution", ylim=(-4, 5.5))
plot!(; size=(400, 250))

Out[3]:

┌ Warning: `vendor()` is deprecated, use `BLAS.get_config()` and inspect the output instead
│   caller = npyinitialize() at numpy.jl:67
└ @ PyCall C:\Users\genkuroki\.julia\packages\PyCall\BD546\src\numpy.jl:67

In [4]:

# total energy
E(v, u, p, t) = ((g, l) = p; v[1]^2/2 - (g/l)*(cos(u[1]) - 1))

Out[4]:

E (generic function with 1 method)

In [5]:

# plot function
function plot_pendulum(sol, ymax; title="solution: θ(t)")
    @unpack tspan, p = sol.prob
    t = range(tspan..., length=501)
    dθ = (t -> sol(t)[1]).(t)
    θ = (t -> sol(t)[2]).(t)
    θ_exact = θ_exact_sol.(Ref(p), t, ymax)
    
    P = plot(t, [θ θ_exact]; label=["numerical" "exact"], ls=:auto, ylim=(-4, 6.5))
    plot!(xtick=0:2:30, ytick=-10:10)
    title!(title; titlefontsize=12)
    
    v0, u0 = first.((dθ, θ))
    E_exact = fill(E(v0, u0, p, 0.0), length(t))
    E_numerical = E.(dθ, θ, Ref(p), t)
    Q = plot(t, E_numerical - E_exact; label="")
    title!("total energy error"; titlefontsize=12)
    plot!(xtick=0:2:30)
    
    plot(P, Q; size=(800, 250), layout=(1, 2))
end

Out[5]:

plot_pendulum (generic function with 1 method)

ODEProblem¶

ODEProblem: in-place version¶

In [6]:

# u = [dθ/dt, θ]
# du = [d²θ/dt², dθ/dt]
function f!(du, u, p, t)
    g, l = p
    du[1] = -(g/l)*sin(u[2])
    du[2] = u[1]
end

ode_prob = ODEProblem(f!, vu0, tspan, p)
@time sol = solve(ode_prob)
@time sol = solve(ode_prob)
@time sol = solve(ode_prob)
@show length(sol.t)

plot_pendulum(sol, ymax)

 24.767610 seconds (57.54 M allocations: 3.365 GiB, 3.73% gc time)
  0.000175 seconds (1.16 k allocations: 122.750 KiB)
  0.000151 seconds (1.16 k allocations: 122.750 KiB)
length(sol.t) = 113

Out[6]:

In [7]:

# u = [dθ/dt, θ]
# du = [d²θ/dt², dθ/dt]
function f!(du, u, p, t)
    g, l = p
    du[1] = -(g/l)*sin(u[2])
    du[2] = u[1]
end

ode_prob = ODEProblem(f!, vu0, tspan, p)
@time sol = solve(ode_prob, RKO65(); dt=0.05)
@time sol = solve(ode_prob, RKO65(); dt=0.05)
@time sol = solve(ode_prob, RKO65(); dt=0.05)
@show length(sol.t)

plot_pendulum(sol, ymax)

  2.895777 seconds (6.72 M allocations: 386.049 MiB, 3.10% gc time)
  0.000254 seconds (4.84 k allocations: 487.719 KiB)
  0.000253 seconds (4.84 k allocations: 487.719 KiB)
length(sol.t) = 601

Out[7]:

ODEProblem: SVector version¶

In [8]:

# u = [dθ/dt, θ]
# du = [d²θ/dt², dθ/dt]
function f(u, p, t)
    g, l = p
    du = SVector(
        -(g/l)*sin(u[2]),
        u[1]
    )
end

ode_prob = ODEProblem(f, SVector{2}(vu0), tspan, p)
@time sol = solve(ode_prob)
@time sol = solve(ode_prob)
@time sol = solve(ode_prob)
@show length(sol.t)

plot_pendulum(sol, ymax)

  3.441609 seconds (10.25 M allocations: 628.632 MiB, 4.34% gc time)
  0.000074 seconds (187 allocations: 40.328 KiB)
  0.000068 seconds (187 allocations: 40.328 KiB)
length(sol.t) = 114

Out[8]:

In [9]:

# u = [dθ/dt, θ]
# du = [d²θ/dt², dθ/dt]
function f(u, p, t)
    g, l = p
    du = SVector(
        -(g/l)*sin(u[2]),
        u[1]
    )
end

ode_prob = ODEProblem(f, SVector{2}(vu0), tspan, p)
@time sol = solve(ode_prob, RKO65(); dt=0.05)
@time sol = solve(ode_prob, RKO65(); dt=0.05)
@time sol = solve(ode_prob, RKO65(); dt=0.05)
@show length(sol.t)

plot_pendulum(sol, ymax)

  0.651616 seconds (2.57 M allocations: 153.325 MiB, 5.61% gc time)
  0.000243 seconds (622 allocations: 124.797 KiB)
  0.000112 seconds (622 allocations: 124.797 KiB)
length(sol.t) = 601

Out[9]:

SecondOrderODEProblem¶

In [10]:

# u = [θ], v = [dθ/dt]
# du = [dθ/dt], dv = [d²θ/dt²]
f2!(dv, v, u, p, t) = ((g, l) = p; dv[1] = -(g/l)*sin(u[1]))

soode_prob = SecondOrderODEProblem(f2!, v0, u0, tspan, p)
@time sol = solve(soode_prob)
@time sol = solve(soode_prob)
@time sol = solve(soode_prob)
@show length(sol.t)

plot_pendulum(sol, ymax)

  5.784741 seconds (26.73 M allocations: 1.312 GiB, 6.42% gc time)
  0.000185 seconds (2.01 k allocations: 210.750 KiB)
  0.000179 seconds (2.01 k allocations: 210.750 KiB)
length(sol.t) = 113

Out[10]:

In [11]:

# u = [θ], v = [dθ/dt]
# du = [dθ/dt], dv = [d²θ/dt²]
f2!(dv, v, u, p, t) = ((g, l) = p; dv[1] = -(g/l)*sin(u[1]))

soode_prob = SecondOrderODEProblem(f2!, v0, u0, tspan, p)
@time sol = solve(soode_prob, Yoshida6(); dt=0.1)
@time sol = solve(soode_prob, Yoshida6(); dt=0.1)
@time sol = solve(soode_prob, Yoshida6(); dt=0.1)
@show length(sol.t)

plot_pendulum(sol, ymax)

  1.609789 seconds (6.28 M allocations: 325.435 MiB, 4.86% gc time)
  0.000181 seconds (2.14 k allocations: 215.203 KiB)
  0.000211 seconds (2.14 k allocations: 215.203 KiB)
length(sol.t) = 301

Out[11]:

DynamicalODEProblem¶

In [12]:

# u = [θ], v = [dθ/dt]
# du = [dθ/dt], dv = [d²θ/dt²]
f1!(du, v, u, p, t) = @. du = v
f2!(dv, v, u, p, t) = ((g, l) = p; dv[1] = -(g/l)*sin(u[1]))

dynode_prob = DynamicalODEProblem(f2!, f1!, v0, u0, tspan, p)
@time sol = solve(dynode_prob)
@time sol = solve(dynode_prob)
@time sol = solve(dynode_prob)
@show length(sol.t)

plot_pendulum(sol, ymax)

  3.786276 seconds (7.95 M allocations: 421.076 MiB, 11.20% gc time)
  0.000202 seconds (2.01 k allocations: 210.750 KiB)
  0.000185 seconds (2.01 k allocations: 210.750 KiB)
length(sol.t) = 113

Out[12]:

In [13]:

# u = [θ], v = [dθ/dt]
# du = [dθ/dt], dv = [d²θ/dt²]
f1!(du, v, u, p, t) = @. du = v
f2!(dv, v, u, p, t) = ((g, l) = p; dv[1] = -(g/l)*sin(u[1]))

dynode_prob = DynamicalODEProblem(f2!, f1!, v0, u0, tspan, p)
@time sol = solve(dynode_prob, Yoshida6(); dt=0.1)
@time sol = solve(dynode_prob, Yoshida6(); dt=0.1)
@time sol = solve(dynode_prob, Yoshida6(); dt=0.1)
@show length(sol.t)

plot_pendulum(sol, ymax)

  1.197207 seconds (2.85 M allocations: 153.838 MiB, 1.93% gc time)
  0.000305 seconds (2.14 k allocations: 215.203 KiB)
  0.000161 seconds (2.14 k allocations: 215.203 KiB)
length(sol.t) = 301

Out[13]:

HamiltonianProblem¶

In [14]:

# u = [θ], v = [dθ/dt]
H(v, u, p, t) = ((g, l) = p; v[1]^2/2 - (g/l)*(cos(u[1]) + 1))

hamil_prob = HamiltonianProblem(H, v0, u0, tspan, p)
@time sol = solve(hamil_prob)
@time sol = solve(hamil_prob)
@time sol = solve(hamil_prob)
@show length(sol.t)

plot_pendulum(sol, ymax)

  3.876122 seconds (9.92 M allocations: 534.090 MiB, 3.83% gc time)
  0.000228 seconds (2.01 k allocations: 211.141 KiB)
  0.000214 seconds (2.01 k allocations: 211.141 KiB)
length(sol.t) = 113

Out[14]:

In [15]:

# u = [θ], v = [dθ/dt]
H(v, u, p, t) = ((g, l) = p; v[1]^2/2 - (g/l)*(cos(u[1]) + 1))

hamil_prob = HamiltonianProblem(H, v0, u0, tspan, p)
@time sol = solve(hamil_prob, Yoshida6(); dt=0.1)
@time sol = solve(hamil_prob, Yoshida6(); dt=0.1)
@time sol = solve(hamil_prob, Yoshida6(); dt=0.1)
@show length(sol.t)

plot_pendulum(sol, ymax)

  1.235402 seconds (2.95 M allocations: 159.639 MiB, 2.50% gc time)
  0.000278 seconds (2.14 k allocations: 215.516 KiB)
  0.000295 seconds (2.14 k allocations: 215.516 KiB)
length(sol.t) = 301

Out[15]:

Simple Runge-Kutta Order 4¶

In [16]:

function f!(du, u, p, t)
    g, l = p
    du[1] = -(g/l)*sin(u[2])
    du[2] = u[1]
end

p = (g, l) = (9.80665, 1.0)
ymax = 0.9999
v0 = [√(g/l)*2ymax]
u0 = [0.0]
tspan = (0.0, 30.0)

Out[16]:

(0.0, 30.0)

In [17]:

function my_simple_RK4(f!, u0, tspan, p; dt=0.1)
    t = range(tspan...; step=dt)
    n = length(t)
    u = similar(typeof(u0)[], n)
    u[1] = u0
    du1, du2, du3, du4 = similar(u0), similar(u0), similar(u0), similar(u0)
    dt_half = dt/2
    for k in 2:n
        u_prev, t_prev = u[k-1], t[k-1]
        f!(du1, u_prev              , p, t_prev          )
        f!(du2, u_prev + du1*dt_half, p, t_prev + dt_half)
        f!(du3, u_prev + du2*dt_half, p, t_prev + dt_half)
        f!(du4, u_prev + du3*dt     , p, t_prev + dt     )
        u_next = @. u_prev + (du1 + 2du2 + 2du3 + du4)*dt/6
        u[k] = u_next
    end
    (; t, u)
end

Out[17]:

my_simple_RK4 (generic function with 1 method)

In [18]:

using SimpleDiffEq

dt = 0.1
sol = solve(ODEProblem(f!, vu0, tspan, p), SimpleRK4(); dt)
my_sol = my_simple_RK4(f!, vu0, tspan, p; dt)

@show sol.t == my_sol.t
@show sol.u ≈ my_sol.u

plot_pendulum(sol, ymax; title="simple Ringe-Kutta order 4 (dt = $dt)")

sol.t == my_sol.t = true
sol.u ≈ my_sol.u = true

Out[18]:

In [19]:

using SimpleDiffEq

dt = 0.05
sol = solve(ODEProblem(f!, vu0, tspan, p), SimpleRK4(); dt)
my_sol = my_simple_RK4(f!, vu0, tspan, p; dt)

@show sol.t == my_sol.t
@show sol.u ≈ my_sol.u

plot_pendulum(sol, ymax; title="simple Ringe-Kutta order 4 (dt = $dt)")

sol.t == my_sol.t = true
sol.u ≈ my_sol.u = true

Out[19]:

In [20]:

using SimpleDiffEq

dt = 0.1
sol_simpleRK4 = solve(ODEProblem(f!, vu0, tspan, p), SimpleRK4(); dt)
sol_RK4 = solve(ODEProblem(f!, vu0, tspan, p), RK4(); dtmax=1.1dt)

@show length(sol_simpleRK4.t)
@show length(sol_RK4.t)

t = range(tspan..., length=1001)
dθ_simpleRK4 = (t -> sol_simpleRK4(t)[1]).(t)
θ_simpleRK4 = (t -> sol_simpleRK4(t)[2]).(t)
dθ_RK4 = (t -> sol_RK4(t)[1]).(t)
θ_RK4 = (t -> sol_RK4(t)[2]).(t)
θ_exact = θ_exact_sol.(Ref(p), t, ymax)

P = plot(t, [θ_simpleRK4 θ_RK4 θ_exact];
    label=["SimpleRK4" "RK4" "exact"], ls=[:solid  :dash :dashdot], ylim=(-3.5, 7))
plot!(xtick=0:2:30, ytick=-10:10)
title!("Simple and adaptive stepping order-4 Runge-Kutta's (dt = $dt)"; titlefontsize=8)

E_exact = fill(E(v0, u0, p, 0.0), length(t))
E_simpleRK4 = E.(dθ_simpleRK4, θ_simpleRK4, Ref(p), t)
E_RK4 = E.(dθ_RK4, θ_RK4, Ref(p), t)

Q1 = plot(t, E_simpleRK4 - E_exact; label="") 
title!("total energy error of SimpleRK4"; titlefontsize=8)
plot!(xtick=0:2:30)
plot!(ylim=(-0.05, 0.005))

Q2 = plot(t, E_RK4 - E_exact; label="") 
title!("total energy error of RK4"; titlefontsize=8)
plot!(xtick=0:2:30)
plot!(ylim=(-0.05, 0.005))

plot(P, Q1, Q2; size=(800, 250), layout=@layout[a [b; c]], tickfontsize=7)

length(sol_simpleRK4.t) = 301
length(sol_RK4.t) = 299

Out[20]:

In [21]:

using SimpleDiffEq

dt = 0.05
sol_simpleRK4 = solve(ODEProblem(f!, vu0, tspan, p), SimpleRK4(); dt)
sol_RK4 = solve(ODEProblem(f!, vu0, tspan, p), RK4(); dtmax=1.01dt)

@show length(sol_simpleRK4.t)
@show length(sol_RK4.t)

t = range(tspan..., length=1001)
dθ_simpleRK4 = (t -> sol_simpleRK4(t)[1]).(t)
θ_simpleRK4 = (t -> sol_simpleRK4(t)[2]).(t)
dθ_RK4 = (t -> sol_RK4(t)[1]).(t)
θ_RK4 = (t -> sol_RK4(t)[2]).(t)
θ_exact = θ_exact_sol.(Ref(p), t, ymax)

P = plot(t, [θ_simpleRK4 θ_RK4 θ_exact];
    label=["SimpleRK4" "RK4" "exact"], ls=[:solid  :dash :dashdot], ylim=(-3.5, 7))
plot!(xtick=0:2:30, ytick=-10:10)
title!("Simple and adaptive stepping order-4 Runge-Kutta's (dt = $dt)"; titlefontsize=8)

E_exact = fill(E(v0, u0, p, 0.0), length(t))
E_simpleRK4 = E.(dθ_simpleRK4, θ_simpleRK4, Ref(p), t)
E_RK4 = E.(dθ_RK4, θ_RK4, Ref(p), t)

Q1 = plot(t, E_simpleRK4 - E_exact; label="") 
title!("total energy error of SimpleRK4"; titlefontsize=8)
plot!(xtick=0:2:30)
plot!(ylim=(-0.001, 0.0001))

Q2 = plot(t, E_RK4 - E_exact; label="") 
title!("total energy error of RK4"; titlefontsize=8)
plot!(xtick=0:2:30)
plot!(ylim=(-0.001, 0.0001))

plot(P, Q1, Q2; size=(800, 250), layout=@layout[a [b; c]], tickfontsize=7)

length(sol_simpleRK4.t) = 601
length(sol_RK4.t) = 600

Out[21]:

Benchmark¶

In [22]:

f2!(dv, v, u, p, t) = ((g, l) = p; dv[1] = -(g/l)*sin(u[1]))

p = (g, l) = (9.80665, 1.0)
ymax = 0.9999
v0 = [√(g/l)*2ymax]
u0 = [0.0]
tspan = (0.0, 30.0)

dt = 0.1

Out[22]:

0.1

In [23]:

Explicit_RungeKutta = [
    # Euler
    Midpoint
    Heun
    Ralston
    RK4
    BS3
    OwrenZen3
    OwrenZen4
    OwrenZen5
    DP5
    Tsit5
    # Anas5
    FRK65
    PFRK87
    # RKO65
    TanYam7
    DP8
    TsitPap8
    Feagin10
    Feagin12
    Feagin14
    BS5
    Vern6
    Vern7
    Vern8
    Vern9
]

Out[23]:

23-element Vector{Type}:
 Midpoint
 Heun
 Ralston
 RK4
 BS3
 OwrenZen3
 OwrenZen4
 OwrenZen5
 DP5
 Tsit5
 FRK65
 PFRK87
 TanYam7
 DP8
 TsitPap8
 Feagin10
 Feagin12
 Feagin14
 BS5
 Vern6
 Vern7
 Vern8
 Vern9

In [24]:

for integrator in Explicit_RungeKutta
    soode_prob = SecondOrderODEProblem(f2!, v0, u0, tspan, p)
    sol = solve(soode_prob, integrator())
    @printf "%28s:" integrator
    @time sol = solve(soode_prob, integrator())
    plot_pendulum(sol, ymax; title="$integrator()") |> display
end

                    Midpoint:  0.001479 seconds (16.75 k allocations: 1.762 MiB)
                        Heun:  0.001125 seconds (15.79 k allocations: 1.672 MiB)
                     Ralston:  0.001233 seconds (16.77 k allocations: 1.764 MiB)
                         RK4:  0.000428 seconds (2.06 k allocations: 219.203 KiB)
                         BS3:  0.000290 seconds (2.60 k allocations: 271.969 KiB)
                   OwrenZen3:  0.000636 seconds (7.94 k allocations: 825.234 KiB)
                   OwrenZen4:  0.000323 seconds (3.73 k allocations: 392.547 KiB)
                   OwrenZen5:  0.000611 seconds (3.17 k allocations: 327.203 KiB)
                         DP5:  0.000168 seconds (1.23 k allocations: 130.766 KiB)
                       Tsit5:  0.000226 seconds (1.98 k allocations: 208.984 KiB)
                       FRK65:  7.280693 seconds (21.00 M allocations: 2.050 GiB, 69.38% gc time)
                      PFRK87:  2.286901 seconds (7.00 M allocations: 710.817 MiB)
                     TanYam7:  0.000217 seconds (713 allocations: 74.141 KiB)
                         DP8:  0.000419 seconds (1.16 k allocations: 118.516 KiB)
                    TsitPap8:  0.000159 seconds (540 allocations: 58.547 KiB)
                    Feagin10:  0.000236 seconds (612 allocations: 66.859 KiB)
                    Feagin12:  0.000178 seconds (664 allocations: 75.797 KiB)
                    Feagin14:  0.000355 seconds (750 allocations: 90.844 KiB)
                         BS5:  0.000158 seconds (1.50 k allocations: 156.688 KiB)
                       Vern6:  0.000221 seconds (2.04 k allocations: 211.516 KiB)
                       Vern7:  0.000236 seconds (2.05 k allocations: 215.609 KiB)
                       Vern8:  0.000193 seconds (1.70 k allocations: 184.609 KiB)
                       Vern9:  0.000329 seconds (1.58 k allocations: 176.250 KiB)

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ SciMLBase C:\Users\genkuroki\.julia\packages\SciMLBase\grNUR\src\integrator_interface.jl:331
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ SciMLBase C:\Users\genkuroki\.julia\packages\SciMLBase\grNUR\src\integrator_interface.jl:331
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ SciMLBase C:\Users\genkuroki\.julia\packages\SciMLBase\grNUR\src\integrator_interface.jl:331
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ SciMLBase C:\Users\genkuroki\.julia\packages\SciMLBase\grNUR\src\integrator_interface.jl:331

In [25]:

RKN_integrator = [
    Nystrom4
    IRKN3
    IRKN4
    ERKN4
    ERKN5
    Nystrom4VelocityIndependent
    Nystrom5VelocityIndependent
    DPRKN6
    DPRKN8
    DPRKN12
]

Out[25]:

10-element Vector{DataType}:
 Nystrom4
 IRKN3
 IRKN4
 ERKN4
 ERKN5
 Nystrom4VelocityIndependent
 Nystrom5VelocityIndependent
 DPRKN6
 DPRKN8
 DPRKN12

See https://diffeq.sciml.ai/stable/solvers/dynamical_solve/

In [26]:

for integrator in RKN_integrator
    soode_prob = SecondOrderODEProblem(f2!, v0, u0, tspan, p)
    sol = solve(soode_prob, integrator(); dt=0.1)
    @printf "%28s:" integrator
    @time sol = solve(soode_prob, integrator(); dt)
    plot_pendulum(sol, ymax; title="$integrator(dt=$dt)") |> display
end

                    Nystrom4:  0.000207 seconds (2.15 k allocations: 215.594 KiB)
                       IRKN3:  0.000135 seconds (2.15 k allocations: 215.906 KiB)
                       IRKN4:  0.000183 seconds (2.17 k allocations: 217.141 KiB)
                       ERKN4:  0.000061 seconds (558 allocations: 57.750 KiB)
                       ERKN5:  0.000080 seconds (642 allocations: 65.703 KiB)
 Nystrom4VelocityIndependent:  0.000173 seconds (2.15 k allocations: 215.422 KiB)
 Nystrom5VelocityIndependent:  0.000144 seconds (2.16 k allocations: 216.266 KiB)
                      DPRKN6:  0.000124 seconds (1.34 k allocations: 143.172 KiB)
                      DPRKN8:  0.000074 seconds (599 allocations: 63.734 KiB)
                     DPRKN12:  0.000073 seconds (345 allocations: 41.312 KiB)

In [27]:

symplectic_integrator = [
    SymplecticEuler
    VelocityVerlet
    VerletLeapfrog
    PseudoVerletLeapfrog
    McAte2
    Ruth3
    McAte3
    CandyRoz4
    McAte4
    CalvoSanz4
    McAte42
    McAte5
    Yoshida6
    KahanLi6
    McAte8
    KahanLi8
    SofSpa10
]

Out[27]:

17-element Vector{DataType}:
 SymplecticEuler
 VelocityVerlet
 VerletLeapfrog
 PseudoVerletLeapfrog
 McAte2
 Ruth3
 McAte3
 CandyRoz4
 McAte4
 CalvoSanz4
 McAte42
 McAte5
 Yoshida6
 KahanLi6
 McAte8
 KahanLi8
 SofSpa10