diff --git a/benchmarks/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.jmd b/benchmarks/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.jmd index c39664b80..56268b595 100644 --- a/benchmarks/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.jmd +++ b/benchmarks/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.jmd @@ -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] @@ -64,11 +64,11 @@ 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) ``` @@ -76,8 +76,8 @@ 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 @@ -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 @@ -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]) -``` +``` \ No newline at end of file diff --git a/benchmarks/DynamicalODE/Manifest.toml b/benchmarks/DynamicalODE/Manifest.toml index 8149756f4..6dfb07988 100644 --- a/benchmarks/DynamicalODE/Manifest.toml +++ b/benchmarks/DynamicalODE/Manifest.toml @@ -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" \ No newline at end of file diff --git a/benchmarks/DynamicalODE/Project.toml b/benchmarks/DynamicalODE/Project.toml index 937e84059..09c185ffd 100644 --- a/benchmarks/DynamicalODE/Project.toml +++ b/benchmarks/DynamicalODE/Project.toml @@ -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" +PyPlot = "2.11.5" +SciMLBenchmarks = "0.1.3" +StaticArrays = "1.9.7" +Statistics = "1.11.1" +TaylorIntegration = "0.16.1" \ No newline at end of file diff --git a/benchmarks/DynamicalODE/single_pendulums.jmd b/benchmarks/DynamicalODE/single_pendulums.jmd index 318c5d9eb..fe6302af9 100644 --- a/benchmarks/DynamicalODE/single_pendulums.jmd +++ b/benchmarks/DynamicalODE/single_pendulums.jmd @@ -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"^[^.]*\." => "") @@ -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() @@ -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)) @@ -257,4 +271,4 @@ singlependulum(k, SymplecticEuler(), Δt) ```julia, echo = false using SciMLBenchmarks SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file]) -``` +``` \ No newline at end of file