From a4570e07af79d44642a3cf2d4312cc0d605e78d1 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Wed, 25 Sep 2024 22:52:46 +0530 Subject: [PATCH 1/9] Changes --- benchmarks/DynamicalODE/Project.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/benchmarks/DynamicalODE/Project.toml b/benchmarks/DynamicalODE/Project.toml index 937e84059..9f60b4f57 100644 --- a/benchmarks/DynamicalODE/Project.toml +++ b/benchmarks/DynamicalODE/Project.toml @@ -14,7 +14,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42" [compat] -DiffEqCallbacks = "3" +DiffEqCallbacks = "2.13, 3" DiffEqPhysics = "3.2" Elliptic = "1.0" OrdinaryDiffEq = "6" @@ -24,3 +24,4 @@ PyPlot = "2.9" SciMLBenchmarks = "0.1" StaticArrays = "1.0" Statistics = "1" +TaylorIntegration = "0.16" \ No newline at end of file From dfedecebb3ed117fa783a750fbf4a27f44f074fe Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 15 Oct 2024 01:14:43 +0530 Subject: [PATCH 2/9] Updates --- benchmarks/DynamicalODE/Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/DynamicalODE/Project.toml b/benchmarks/DynamicalODE/Project.toml index 9f60b4f57..a31c389a5 100644 --- a/benchmarks/DynamicalODE/Project.toml +++ b/benchmarks/DynamicalODE/Project.toml @@ -14,7 +14,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42" [compat] -DiffEqCallbacks = "2.13, 3" +DiffEqCallbacks = "4.0.0" DiffEqPhysics = "3.2" Elliptic = "1.0" OrdinaryDiffEq = "6" @@ -24,4 +24,4 @@ PyPlot = "2.9" SciMLBenchmarks = "0.1" StaticArrays = "1.0" Statistics = "1" -TaylorIntegration = "0.16" \ No newline at end of file +TaylorIntegration = "0.16.1" \ No newline at end of file From 7be8b4a5390ea04dcd6a2652cbaae181e927df01 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Wed, 16 Oct 2024 00:19:05 +0530 Subject: [PATCH 3/9] Fixes --- ...n-Heiles_energy_conservation_benchmark.jmd | 26 +++++++++---------- benchmarks/DynamicalODE/single_pendulums.jmd | 2 +- 2 files changed, 13 insertions(+), 15 deletions(-) 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/single_pendulums.jmd b/benchmarks/DynamicalODE/single_pendulums.jmd index 318c5d9eb..de4c2df86 100644 --- a/benchmarks/DynamicalODE/single_pendulums.jmd +++ b/benchmarks/DynamicalODE/single_pendulums.jmd @@ -257,4 +257,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 From fb171fcbae0cbb2b2017e55865a891407fa9b341 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Wed, 16 Oct 2024 00:24:23 +0530 Subject: [PATCH 4/9] updates --- benchmarks/DynamicalODE/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/DynamicalODE/Project.toml b/benchmarks/DynamicalODE/Project.toml index a31c389a5..de3d715b4 100644 --- a/benchmarks/DynamicalODE/Project.toml +++ b/benchmarks/DynamicalODE/Project.toml @@ -14,7 +14,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42" [compat] -DiffEqCallbacks = "4.0.0" +DiffEqCallbacks = "2.13, 3, 4" DiffEqPhysics = "3.2" Elliptic = "1.0" OrdinaryDiffEq = "6" From 22886280519f34ae4b9cff0383e6f657a05df776 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Fri, 18 Oct 2024 13:13:55 +0530 Subject: [PATCH 5/9] Updates --- benchmarks/DynamicalODE/single_pendulums.jmd | 45 ++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/benchmarks/DynamicalODE/single_pendulums.jmd b/benchmarks/DynamicalODE/single_pendulums.jmd index de4c2df86..e94d25b4b 100644 --- a/benchmarks/DynamicalODE/single_pendulums.jmd +++ b/benchmarks/DynamicalODE/single_pendulums.jmd @@ -179,6 +179,34 @@ end ```julia # 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 k = rand() integrator = VelocityVerlet() @@ -189,6 +217,23 @@ singlependulum(k, integrator, Δt, t0=-20.0, t1=20.0) ```julia # Two single pendulums +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 + H(q,p,param) = sum(p.^2/2 .- cos.(q) .+ 1) q0 = pi*rand(2) p0 = zeros(2) From 61c475fcfa3cd9cbe399998dffd53f9831efd220 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Fri, 18 Oct 2024 17:44:34 +0530 Subject: [PATCH 6/9] Updates --- benchmarks/DynamicalODE/Project.toml | 3 +-- benchmarks/DynamicalODE/single_pendulums.jmd | 4 ++-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/benchmarks/DynamicalODE/Project.toml b/benchmarks/DynamicalODE/Project.toml index de3d715b4..8190c62de 100644 --- a/benchmarks/DynamicalODE/Project.toml +++ b/benchmarks/DynamicalODE/Project.toml @@ -19,8 +19,7 @@ DiffEqPhysics = "3.2" Elliptic = "1.0" OrdinaryDiffEq = "6" ParameterizedFunctions = "5.3" -Plots = "1.4" -PyPlot = "2.9" +Plots = "2.0.0" SciMLBenchmarks = "0.1" StaticArrays = "1.0" Statistics = "1" diff --git a/benchmarks/DynamicalODE/single_pendulums.jmd b/benchmarks/DynamicalODE/single_pendulums.jmd index e94d25b4b..cf8ecb2a1 100644 --- a/benchmarks/DynamicalODE/single_pendulums.jmd +++ b/benchmarks/DynamicalODE/single_pendulums.jmd @@ -48,9 +48,9 @@ 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 Plot colorlist = [ "#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", From cd7aae2553e9c7791e2ae3020f1fa3b97c09ada1 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Fri, 18 Oct 2024 17:46:13 +0530 Subject: [PATCH 7/9] Updates --- benchmarks/DynamicalODE/Project.toml | 3 ++- benchmarks/DynamicalODE/single_pendulums.jmd | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/benchmarks/DynamicalODE/Project.toml b/benchmarks/DynamicalODE/Project.toml index 8190c62de..de3d715b4 100644 --- a/benchmarks/DynamicalODE/Project.toml +++ b/benchmarks/DynamicalODE/Project.toml @@ -19,7 +19,8 @@ DiffEqPhysics = "3.2" Elliptic = "1.0" OrdinaryDiffEq = "6" ParameterizedFunctions = "5.3" -Plots = "2.0.0" +Plots = "1.4" +PyPlot = "2.9" SciMLBenchmarks = "0.1" StaticArrays = "1.0" Statistics = "1" diff --git a/benchmarks/DynamicalODE/single_pendulums.jmd b/benchmarks/DynamicalODE/single_pendulums.jmd index cf8ecb2a1..875ace816 100644 --- a/benchmarks/DynamicalODE/single_pendulums.jmd +++ b/benchmarks/DynamicalODE/single_pendulums.jmd @@ -50,7 +50,7 @@ sn(u, k) = Jacobi.sn(u, k^2) # the Jacobian sn function # Use Plot. # -using Plot +using PyPlot colorlist = [ "#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", @@ -179,6 +179,7 @@ end ```julia # Single pendulum +using PyPlot 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] From a146fdb5a7042f8feb1e756ef38ab5568d4b6c3a Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Sat, 19 Oct 2024 19:39:04 +0530 Subject: [PATCH 8/9] Changes Made --- benchmarks/DynamicalODE/Manifest.toml | 2 +- benchmarks/DynamicalODE/Project.toml | 18 +- benchmarks/DynamicalODE/single_pendulums.jmd | 172 ++++++++----------- 3 files changed, 80 insertions(+), 112 deletions(-) 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 de3d715b4..09c185ffd 100644 --- a/benchmarks/DynamicalODE/Project.toml +++ b/benchmarks/DynamicalODE/Project.toml @@ -15,13 +15,13 @@ TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42" [compat] DiffEqCallbacks = "2.13, 3, 4" -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" +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 875ace816..87c70b965 100644 --- a/benchmarks/DynamicalODE/single_pendulums.jmd +++ b/benchmarks/DynamicalODE/single_pendulums.jmd @@ -50,128 +50,100 @@ sn(u, k) = Jacobi.sn(u, k^2) # the Jacobian sn function # 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"^[^.]*\." => "") - - figure(figsize=(10,8)) + integrator_str = replace("$integrator", r"^[^.]*\\." => "") - subplot2grid((21,20), ( 1, 0), rowspan=10, colspan=10) - plotsol(sol) + p1 = plotsol(sol) + p2 = plotsol2(sol) + p3 = plotenergy(H, 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 =====") + suptitle = "===== $integrator_str, Δt = $Δt =====" + plot(p1, p2, p3, layout=(3, 1), title=suptitle) 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)) + H(p, q, params) = p[1]^2 / 2 - cos(q[1]) + 1 + q0 = [0.0] + p0 = [2k] + prob = HamiltonianProblem(H, p0, q0, (t0, t1)) - local integrator_str = replace("$integrator", r"^[^.]*\." => "") + integrator_str = replace("$integrator", r"^[^.]*\\." => "") @printf("%-25s", "$integrator_str:") sol = solve(prob, integrator, dt=Δt) - @time local sol = solve(prob, integrator, dt=Δt) + @time 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 =====") + + 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 ``` @@ -179,9 +151,11 @@ end ```julia # Single pendulum -using PyPlot +using Plots +using OrdinaryDiffEq + 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 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)) @@ -192,21 +166,16 @@ 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) + # 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") + p4 = plot(sol.t, sol.u[1] .- k, label="Comparison", title="Comparison Plot") - 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 =====") + # 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 k = rand() @@ -218,25 +187,24 @@ singlependulum(k, integrator, Δt, t0=-20.0, t1=20.0) ```julia # Two single pendulums +using Plots +using OrdinaryDiffEq + 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) + # 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") - suptitle("===== $integrator_str, Δt = $Δt =====") + # 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) +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)) From 5535ed2aee09e84323a0f84359ea0787270a9f21 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Sat, 19 Oct 2024 21:18:18 +0530 Subject: [PATCH 9/9] Updates --- benchmarks/DynamicalODE/single_pendulums.jmd | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/benchmarks/DynamicalODE/single_pendulums.jmd b/benchmarks/DynamicalODE/single_pendulums.jmd index 87c70b965..fe6302af9 100644 --- a/benchmarks/DynamicalODE/single_pendulums.jmd +++ b/benchmarks/DynamicalODE/single_pendulums.jmd @@ -155,9 +155,9 @@ using Plots using OrdinaryDiffEq 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"^[^.]*\." => "") @@ -168,10 +168,10 @@ function singlependulum(k, integrator, Δt; t0 = 0.0, t1 = 100.0) sleep(0.1) # 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") - p4 = plot(sol.t, sol.u[1] .- k, label="Comparison", title="Comparison Plot") + 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") # Combine all plots in a layout layout = @layout [a b; c d]