Skip to content

Commit

Permalink
Updates
Browse files Browse the repository at this point in the history
  • Loading branch information
@ParamThakkar123
ParamThakkar123 committed Oct 20, 2024
2 parents b67b885 + 43bb152 commit ed748c2
Show file tree
Hide file tree
Showing 6 changed files with 284 additions and 234 deletions.
26 changes: 12 additions & 14 deletions benchmarks/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.jmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ end
const oop_q0 = @SVector [0.1, 0.]
const oop_p0 = @SVector [0., 0.5]

function hamilton(du,u,p,t)
function hamilton(du, u, p, t)
dq, q = @views u[3:4], du[3:4]
dp, p = @views u[1:2], du[1:2]

Expand All @@ -64,20 +64,20 @@ function hamilton(du,u,p,t)
end

function g(resid, u, p)
resid[1] = H([u[1],u[2]], [u[3],u[4]], nothing) - E
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 cb = ManifoldProjection(g, nlopts=Dict(:ftol => 1e-13))

const E = H(iip_p0, iip_q0, nothing)
```

For the comparison we will use the following function

```julia
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)]
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)
if mode == :inplace
Expand All @@ -100,23 +100,21 @@ function compare(mode=:inplace, all=true, plt=nothing; tmax=1e2)
GC.gc()
(mode == :inplace && all) && @time sol6 = solve(prob_linear, 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.")
(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.")
(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)
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
Expand Down Expand Up @@ -205,4 +203,4 @@ The benchmarks were performed on a machine with
```julia, echo = false
using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
```
```
2 changes: 1 addition & 1 deletion benchmarks/DynamicalODE/Manifest.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2345,4 +2345,4 @@ version = "3.5.0+0"
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "9c304562909ab2bab0262639bd4f444d7bc2be37"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "1.4.1+1"
version = "1.4.1+1"
21 changes: 11 additions & 10 deletions benchmarks/DynamicalODE/Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,13 +14,14 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"

[compat]
DiffEqCallbacks = "3"
DiffEqPhysics = "3.2"
Elliptic = "1.0"
OrdinaryDiffEq = "6"
ParameterizedFunctions = "5.3"
Plots = "1.4"
PyPlot = "2.9"
SciMLBenchmarks = "0.1"
StaticArrays = "1.0"
Statistics = "1"
DiffEqCallbacks = "2.13, 3, 4"
DiffEqPhysics = "3.15.0"
Elliptic = "1.0.1"
OrdinaryDiffEq = "6.89.0"
ParameterizedFunctions = "5.17.0"
Plots = "2.0.0, 1"
PyPlot = "2.11.5"
SciMLBenchmarks = "0.1.3"
StaticArrays = "1.9.7"
Statistics = "1.11.1"
TaylorIntegration = "0.16.1"
168 changes: 91 additions & 77 deletions benchmarks/DynamicalODE/single_pendulums.jmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,107 +48,116 @@ sol2tqp(sol) = (sol.t, sol2q(sol), sol2p(sol))
#
sn(u, k) = Jacobi.sn(u, k^2) # the Jacobian sn function

# Use PyPlot.
# Use Plot.
#
using PyPlot
using Plots

# Define the color list
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 to plot the solution of a Hamiltonian problem
function plotsol(sol::ODESolution)
local t, q, p
t, q, p = sol2tqp(sol)
local d = size(q)[2]
p1 = plot(title="Solution", xlabel="t", grid=:on)
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="--")
plot!(p1, t, q[:, j], color=cc(2j-1), label="q$(j_str)", linewidth=1)
plot!(p1, t, p[:, j], color=cc(2j), label="p$(j_str)", linewidth=1, linestyle=:dash)
end
grid(ls=":")
xlabel("t")
legend()
return p1
end

# plot the solution of a Hamiltonian problem on the 2D phase space
#
# Function to plot the solution on the 2D phase space
function plotsol2(sol::ODESolution)
local t, q, p
t, q, p = sol2tqp(sol)
local d = size(q)[2]
p2 = plot(title="Phase Space", xlabel="q", ylabel="p", grid=:on)
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)
plot!(p2, q[:, j], p[:, j], color=cc(j), label="(q$(j_str), p$(j_str))", linewidth=1)
end
grid(ls=":")
xlabel("q")
ylabel("p")
legend()
return p2
end

# plot the energy of a Hamiltonian problem
#
# Function to 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)
p3 = plot(t, energy, label="energy", color="red", linewidth=1, xlabel="t", title="Energy", grid=:on)
stdenergy_str = @sprintf("%.3e", std(energy))
title!("std(energy) = $stdenergy_str", fontsize=10)
return p3
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 to compare the numerical and exact solutions of a single pendulum
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)
local y = sin.(q / 2)
local y_exact = k * sn.(t, k) # the exact solution

p4 = plot(t, y, label="numerical", linewidth=1, grid=:on, xlabel="t", ylabel="y = sin(q(t)/2)", title="Comparison")
plot!(p4, t, y_exact, label="exact", linewidth=1, linestyle=:dash)
error_str = @sprintf("%.3e", maximum(abs.(y - y_exact)))
title!("maximum(abs(numerical - exact)) = $error_str", fontsize=10)
return p4
end

# plot solution and energy
#
# Plot solution and energy
function plotsolenergy(H, integrator, Δt, sol::ODESolution)
local integrator_str = replace("$integrator", r"^[^.]*\." => "")
integrator_str = replace("$integrator", r"^[^.]*\\." => "")

figure(figsize=(10,8))
p1 = plotsol(sol)
p2 = plotsol2(sol)
p3 = plotenergy(H, sol)

subplot2grid((21,20), ( 1, 0), rowspan=10, colspan=10)
plotsol(sol)
suptitle = "===== $integrator_str, Δt = $Δt ====="
plot(p1, p2, p3, layout=(3, 1), title=suptitle)
end

subplot2grid((21,20), ( 1,10), rowspan=10, colspan=10)
plotsol2(sol)
# Solve a single pendulum
function singlependulum(k, integrator, Δt; t0 = 0.0, t1 = 100.0)
H(p, q, params) = p[1]^2 / 2 - cos(q[1]) + 1
q0 = [0.0]
p0 = [2k]
prob = HamiltonianProblem(H, p0, q0, (t0, t1))

subplot2grid((21,20), (11, 0), rowspan=10, colspan=10)
plotenergy(H, sol)
integrator_str = replace("$integrator", r"^[^.]*\\." => "")
@printf("%-25s", "$integrator_str:")
sol = solve(prob, integrator, dt=Δt)
@time sol = solve(prob, integrator, dt=Δt)

suptitle("===== $integrator_str, Δt = $Δt =====")
sleep(0.1)

p1 = plotsol(sol)
p2 = plotsol2(sol)
p3 = plotenergy(H, sol)
p4 = plotcomparison(k, sol)

suptitle = "===== $integrator_str, Δt = $Δt ====="
plot(p1, p2, p3, p4, layout=(2, 2), title=suptitle)
end
```

## Tests

```julia
# Single pendulum
using Plots
using OrdinaryDiffEq

# 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 H(p, q, params) = p[1]^2 / 2 - cos(q[1]) + 1 # Hamiltonian for single pendulum
local q0 = [0.0] # Initial position
local p0 = [2k] # Initial momentum
local prob = HamiltonianProblem(H, p0, q0, (t0, t1))

local integrator_str = replace("$integrator", r"^[^.]*\." => "")
Expand All @@ -157,28 +166,17 @@ function singlependulum(k, integrator, Δt; t0 = 0.0, t1 = 100.0)
@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)
# Create plots using Plots.jl
p1 = plot(sol.t, map(x -> x[1], sol.u), label="Position", title="Position vs Time")
p2 = plot(sol.t, map(x -> x[2], sol.u), label="Momentum", title="Momentum vs Time")
p3 = plot(sol.t, map(p -> H(p[1], p[2], nothing), sol.u), label="Energy", title="Energy vs Time")
p4 = plot(sol.t, map(x -> x[1] - k, sol.u), label="Comparison", title="Comparison Plot")

suptitle("===== $integrator_str, Δt = $Δt =====")
# Combine all plots in a layout
layout = @layout [a b; c d]
plot(p1, p2, p3, p4, layout=layout, size=(1000, 800), title="===== $integrator_str, Δt = $Δt =====")
end
```

## Tests

```julia
# Single pendulum

k = rand()
integrator = VelocityVerlet()
Expand All @@ -189,8 +187,24 @@ 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)
using Plots
using OrdinaryDiffEq

function plotsolenergy(H, integrator, Δt, sol::ODESolution)
local integrator_str = replace("$integrator", r"^[^.]*\." => "")

# Create plots using Plots.jl
p1 = plot(sol, label="Position", title="Position vs Time")
p2 = plot(sol, label="Momentum", title="Momentum vs Time")
p3 = plot(sol.t, map(p -> H(p[2], p[3], nothing), sol.u), label="Energy", title="Energy vs Time")

# Combine all plots in a layout
layout = @layout [a b; c]
plot(p1, p2, p3, layout=layout, size=(1000, 800), title="===== $integrator_str, Δt = $Δt =====")
end

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))
Expand Down Expand Up @@ -257,4 +271,4 @@ singlependulum(k, SymplecticEuler(), Δt)
```julia, echo = false
using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
```
```
Loading

0 comments on commit ed748c2

Please sign in to comment.