@@ -155,9 +155,9 @@ using Plots
155155using OrdinaryDiffEq
156156
157157function singlependulum(k, integrator, Δt; t0 = 0.0, t1 = 100.0)
158- local H(p, q, params) = p[1]^2 / 2 - cos(q[1]) + 1
159- local q0 = [0.0]
160- local p0 = [2k]
158+ local H(p, q, params) = p[1]^2 / 2 - cos(q[1]) + 1 # Hamiltonian for single pendulum
159+ local q0 = [0.0] # Initial position
160+ local p0 = [2k] # Initial momentum
161161 local prob = HamiltonianProblem(H, p0, q0, (t0, t1))
162162
163163 local integrator_str = replace("$integrator", r"^[^.]*\." => "")
@@ -168,10 +168,10 @@ function singlependulum(k, integrator, Δt; t0 = 0.0, t1 = 100.0)
168168 sleep(0.1)
169169
170170 # Create plots using Plots.jl
171- p1 = plot(sol, label="Position", title="Position vs Time")
172- p2 = plot(sol, label="Momentum", title="Momentum vs Time")
173- p3 = plot(sol.t, map(p -> H(p[2 ], p[3 ], nothing), sol.u), label="Energy", title="Energy vs Time")
174- p4 = plot(sol.t, sol.u [1] . - k, label="Comparison", title="Comparison Plot")
171+ p1 = plot(sol.t, map(x -> x[1], sol.u) , label="Position", title="Position vs Time")
172+ p2 = plot(sol.t, map(x -> x[2], sol.u) , label="Momentum", title="Momentum vs Time")
173+ p3 = plot(sol.t, map(p -> H(p[1 ], p[2 ], nothing), sol.u), label="Energy", title="Energy vs Time")
174+ p4 = plot(sol.t, map(x -> x [1] - k, sol.u) , label="Comparison", title="Comparison Plot")
175175
176176 # Combine all plots in a layout
177177 layout = @layout [a b; c d]
0 commit comments