Merge pull request #1067 from ParamThakkar123/DynamicalODE

DynamicalODE Changes

Author: ChrisRackauckas

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,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
@@ -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])
-```
+```

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"

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"

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])
-```
+```