# Single Pedulum Comparison

##### Gen Kuroki (黒木玄), Chris Rackauckas

<p><div class="lev1 toc-item"><a href="#Solving-single-pendulums-by-DifferentialEquations.jl" data-toc-modified-id="Solving-single-pendulums-by-DifferentialEquations.jl-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Solving single pendulums by DifferentialEquations.jl</a></div><div class="lev2 toc-item"><a href="#Tests" data-toc-modified-id="Tests-11"><span class="toc-item-num">1.1&nbsp;&nbsp;</span>Tests</a></div><div class="lev2 toc-item"><a href="#Comparison-of-symplectic-Integrators" data-toc-modified-id="Comparison-of-symplectic-Integrators-12"><span class="toc-item-num">1.2&nbsp;&nbsp;</span>Comparison of symplectic Integrators</a></div>

# Solving single pendulums by DifferentialEquations.jl

In this notebook, we shall solve the single pendulum equation:

$\ddot q = -\sin q,$

where $q$ means the angle.

Hamiltonian:

$H(q,p) = \frac{1}{2}p^2 - \cos q + 1.$

Canonical equation:

$\dot q = p, \quad \dot p = - \sin q.$

Initial condition:

$q(0) = 0, \quad p(0) = 2k.$

Exact solution:

$q(t) = 2\arcsin(k\,\mathrm{sn}(t,k)).$

Maximum of $q(t)$:

$\sin(q_{\max}/2) = k, \quad q_{\max} = \max\{q(t)\}.$

Define $y(t)$ by

$y(t) = \sin(q(t)/2) = k\,\mathrm{sn}(t,k), \quad y_{\max} = k.$

# Single pendulums shall be solved numerically.
#
using OrdinaryDiffEq, Elliptic, Printf, DiffEqPhysics, Statistics

sol2q(sol) = [sol.u[i][j] for i in 1:length(sol.u), j in 1:length(sol.u[1])÷2]
sol2p(sol) = [sol.u[i][j] for i in 1:length(sol.u), j in length(sol.u[1])÷2+1:length(sol.u[1])]
sol2tqp(sol) = (sol.t, sol2q(sol), sol2p(sol))

# The exact solutions of single pendulums can be expressed by the Jacobian elliptic functions.
#
sn(u, k) = Jacobi.sn(u, k^2) # the Jacobian sn function

# Use PyPlot.
#
using PyPlot

colorlist = [
"#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
"#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf",
]
cc(k) = colorlist[mod1(k, length(colorlist))]

# plot the sulution of a Hamiltonian problem
#
function plotsol(sol::ODESolution)
local t, q, p
t, q, p = sol2tqp(sol)
local d = size(q)[2]
for j in 1:d
j_str = d > 1 ? "[$j]" : "" plot(t, q[:,j], color=cc(2j-1), label="q$(j_str)", lw=1)
plot(t, p[:,j], color=cc(2j),   label="p$(j_str)", lw=1, ls="--") end grid(ls=":") xlabel("t") legend() end # plot the solution of a Hamiltonian problem on the 2D phase space # function plotsol2(sol::ODESolution) local t, q, p t, q, p = sol2tqp(sol) local d = size(q)[2] for j in 1:d j_str = d > 1 ? "[$j]" : ""
plot(q[:,j], p[:,j], color=cc(j), label="(q$(j_str),p$(j_str))", lw=1)
end
grid(ls=":")
xlabel("q")
ylabel("p")
legend()
end

# plot the energy of a Hamiltonian problem
#
function plotenergy(H, sol::ODESolution)
local t, q, p
t, q, p = sol2tqp(sol)
local energy = [H(q[i,:], p[i,:], nothing) for i in 1:size(q)[1]]
plot(t, energy, label="energy", color="red", lw=1)
grid(ls=":")
xlabel("t")
legend()
local stdenergy_str = @sprintf("%.3e", std(energy))
title("                    std(energy) = $stdenergy_str", fontsize=10) end # plot the numerical and exact solutions of a single pendulum # # Warning: Assume q(0) = 0, p(0) = 2k. (for the sake of laziness) # function plotcomparison(k, sol::ODESolution) local t, q, p t, q, p = sol2tqp(sol) local y = sin.(q/2) local y_exact = k*sn.(t, k) # the exact solution plot(t, y, label="numerical", lw=1) plot(t, y_exact, label="exact", lw=1, ls="--") grid(ls=":") xlabel("t") ylabel("y = sin(q(t)/2)") legend() local error_str = @sprintf("%.3e", maximum(abs.(y - y_exact))) title("maximum(abs(numerical - exact)) =$error_str", fontsize=10)
end

# plot solution and energy
#
function plotsolenergy(H, integrator, Δt, sol::ODESolution)
local integrator_str = replace("$integrator", r"^[^.]*\." => "") figure(figsize=(10,8)) subplot2grid((21,20), ( 1, 0), rowspan=10, colspan=10) plotsol(sol) subplot2grid((21,20), ( 1,10), rowspan=10, colspan=10) plotsol2(sol) subplot2grid((21,20), (11, 0), rowspan=10, colspan=10) plotenergy(H, sol) suptitle("=====$integrator_str,   Δt = $Δt =====") end # Solve a single pendulum # function singlependulum(k, integrator, Δt; t0 = 0.0, t1 = 100.0) local H(p,q,params) = p[1]^2/2 - cos(q[1]) + 1 local q0 = [0.0] local p0 = [2k] local prob = HamiltonianProblem(H, p0, q0, (t0, t1)) local integrator_str = replace("$integrator", r"^[^.]*\." => "")
@printf("%-25s", "$integrator_str:") sol = solve(prob, integrator, dt=Δt) @time local sol = solve(prob, integrator, dt=Δt) sleep(0.1) figure(figsize=(10,8)) subplot2grid((21,20), ( 1, 0), rowspan=10, colspan=10) plotsol(sol) subplot2grid((21,20), ( 1,10), rowspan=10, colspan=10) plotsol2(sol) subplot2grid((21,20), (11, 0), rowspan=10, colspan=10) plotenergy(H, sol) subplot2grid((21,20), (11,10), rowspan=10, colspan=10) plotcomparison(k, sol) suptitle("=====$integrator_str,   Δt = \$Δt    =====")
end

singlependulum (generic function with 1 method)


## Tests

# Single pendulum

k = rand()
integrator = VelocityVerlet()
Δt = 0.1
singlependulum(k, integrator, Δt, t0=-20.0, t1=20.0)

VelocityVerlet():          0.000400 seconds (6.91 k allocations: 359.141 Ki
B)
PyObject Text(0.5, 0.98, '=====    VelocityVerlet(),   Δt = 0.1    =====')

# Two single pendulums

H(q,p,param) = sum(p.^2/2 .- cos.(q) .+ 1)
q0 = pi*rand(2)
p0 = zeros(2)
t0, t1 = -20.0, 20.0
prob = HamiltonianProblem(H, q0, p0, (t0, t1))

integrator = VelocityVerlet()
Δt = 0.1
sol = solve(prob, integrator, dt=Δt)
@time sol = solve(prob, integrator, dt=Δt)

0.001188 seconds (10.13 k allocations: 685.313 KiB)

sleep(0.1)
plotsolenergy(H, integrator, Δt, sol)

PyObject Text(0.5, 0.98, '=====    VelocityVerlet(),   Δt = 0.1    =====')


## Comparison of symplectic Integrators

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

k = 0.999
Δt = 0.1
for integrator in SymplecticIntegrators
singlependulum(k, integrator, Δt)
end

SymplecticEuler():         0.000690 seconds (17.13 k allocations: 936.922 K
iB)
VelocityVerlet():          0.000718 seconds (17.13 k allocations: 936.938 K
iB)
VerletLeapfrog():          0.000827 seconds (18.13 k allocations: 952.688 K
iB)
PseudoVerletLeapfrog():    0.000750 seconds (18.13 k allocations: 952.688 K
iB)
McAte2():                  0.000805 seconds (18.14 k allocations: 952.828 K
iB)
Ruth3():                   0.000989 seconds (20.14 k allocations: 984.188 K
iB)
McAte3():                  0.000857 seconds (20.14 k allocations: 984.109 K
iB)
CandyRoz4():               0.001014 seconds (22.17 k allocations: 1015.906
KiB)
McAte4():                  0.000953 seconds (22.13 k allocations: 1015.172
KiB)
CalvoSanz4():              0.001100 seconds (24.15 k allocations: 1.022 MiB
)
McAte42():                 0.001063 seconds (24.14 k allocations: 1.022 MiB
)
McAte5():                  0.001359 seconds (26.14 k allocations: 1.053 MiB
)
Yoshida6():                0.001230 seconds (30.16 k allocations: 1.114 MiB
)
KahanLi6():                0.001338 seconds (34.14 k allocations: 1.175 MiB
)
McAte8():                  0.001725 seconds (46.16 k allocations: 1.358 MiB
)
KahanLi8():                0.001923 seconds (50.16 k allocations: 1.420 MiB
)
SofSpa10():                0.003019 seconds (86.20 k allocations: 1.970 MiB
)

k = 0.999
Δt = 0.01
for integrator in SymplecticIntegrators[1:4]
singlependulum(k, integrator, Δt)
end

SymplecticEuler():         0.007145 seconds (170.11 k allocations: 8.627 Mi
B)
VelocityVerlet():          0.006243 seconds (170.11 k allocations: 8.627 Mi
B)
VerletLeapfrog():          0.006941 seconds (180.11 k allocations: 8.780 Mi
B)
PseudoVerletLeapfrog():    0.006435 seconds (180.11 k allocations: 8.780 Mi
B)

k = 0.999
Δt = 0.001
singlependulum(k, SymplecticEuler(), Δt)

SymplecticEuler():         0.548139 seconds (1.70 M allocations: 86.218 MiB
, 82.64% gc time)
PyObject Text(0.5, 0.98, '=====    SymplecticEuler(),   Δt = 0.001    =====
')

k = 0.999
Δt = 0.0001
singlependulum(k, SymplecticEuler(), Δt)

SymplecticEuler():         1.583127 seconds (17.00 M allocations: 862.127 M
iB, 39.87% gc time)
PyObject Text(0.5, 0.98, '=====    SymplecticEuler(),   Δt = 0.0001    ====
=')

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","single_pendulums.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.8.0
[78ddff82-25fc-5f2b-89aa-309469cbf16f] DiffEqMonteCarlo 0.14.0
[77a26b50-5914-5dd7-bc55-306e6241c503] DiffEqNoiseProcess 3.2.0
[055956cb-9e8b-5191-98cc-73ae4a59e68a] DiffEqPhysics 3.1.0
[a077e3f3-b75c-5d7f-a0c6-6bc4c8ec64a9] DiffEqProblemLibrary 4.1.0
[41bf760c-e81c-5289-8e54-58b1f1f8abe2] DiffEqSensitivity 3.2.2
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.3.0
[b305315f-e792-5b7a-8f41-49f472929428] Elliptic 0.5.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.2.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
[e88e6eb3-aa80-5325-afca-941959d7151f] Zygote 0.3.0