Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update single_pendulums.jmd #890

Closed
wants to merge 6 commits into from
Closed

Update single_pendulums.jmd #890

Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b2fb938
Update single_pendulums.jmd
Spinachboul Apr 2, 2024
445b32a
Update single_pendulums.jmd
Spinachboul Apr 5, 2024
127fa8d
Update single_pendulums.jmd
Spinachboul Apr 11, 2024
58dc8a3
Update single_pendulums.jmd
Spinachboul Apr 11, 2024
700f80f
Update single_pendulums.jmd
Spinachboul Apr 11, 2024
ef2d917
Update single_pendulums.jmd
Spinachboul Apr 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
237 changes: 49 additions & 188 deletions benchmarks/DynamicalODE/single_pendulums.jmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,221 +37,82 @@ $$y(t) = \sin(q(t)/2) = k\,\mathrm{sn}(t,k), \quad y_{\max} = k.$$

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

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

# Define a color list for plotting
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
```

## Tests

```julia
# Single pendulum

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

```julia
# 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 = McAte4()
Δt = 0.1
sol = solve(prob, integrator, dt=Δt)
@time sol = solve(prob, integrator, dt=Δt)

sleep(0.1)
plotsolenergy(H, integrator, Δt, sol)
cc(k) = colorlist[mod1(k, length(colorist))]
```

## Comparison of symplectic Integrators

```julia
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)
# Plot the solution of a Hamiltonian problem
t, q, p = sol2tqp(sol)
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=:dash, alpha=0.8)
end
plot!(grid=true) # Add grid with appropriate linestyle
xlabel!("t") # Update x-axis label
ylabel!("q, p") # Update y-axis label
plot!(legend=:bottomleft)
```

```julia
k = 0.999
Δt = 0.01
for integrator in SymplecticIntegrators[1:4]
singlependulum(k, integrator, Δt)
# Plot the solution of a Hamiltonian problem on the 2D phase space
t, q, p = sol2tqp(sol)
println("t: ", t)
println("q: ", q)
println("p: ", p)
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
plot!(grid=true, linestyle=:dot) # Add grid with appropriate linestyle
xlabel!("q") # Update x-axis label
ylabel!("p") # Update y-axis label
plot!(legend=:bottomleft)
```

```julia
k = 0.999
Δt = 0.001
singlependulum(k, SymplecticEuler(), Δt)
# Plot the energy of a Hamiltonian problem
t, q, p = sol2tqp(sol)
energy = [H(q[i,:], p[i,:], nothing) for i in 1:size(q)[1]]
plot(t, energy, label="energy", color="red", lw=1)
plot!(grid=true, linestyle=:dot) # Add grid with appropriate linestyle
xlabel!("t") # Update x-axis label
ylabel!("Energy") # Update y-axis label
plot!(legend=:bottomleft)
stdenergy_str = @sprintf("%.3e", std(energy))
title!("std(energy) = $stdenergy_str", fontsize=10) # Update title
```

```julia
k = 0.999
Δt = 0.0001
singlependulum(k, SymplecticEuler(), Δt)
# Plot the numerical and exact solutions of a single pendulum
t, q, p = sol2tqp(sol)
y = sin.(q/2)
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=:dash, alpha=0.8, color="red")
plot!(grid=true) # Add grid with appropriate linestyle
xlabel!("t") # Update x-axis label
ylabel!("y = sin(q(t)/2)") # Update y-axis label
plot!(legend=:bottomleft)
error_str = @sprintf("%.3e", maximum(abs.(y - y_exact)))
title!("maximum(abs(numerical - exact)) = $error_str", fontsize=10)
```

```julia, echo = false
Expand Down