finish problem 3 - almost done with problem 4

1 files changed, 175 insertions, 0 deletions

diff --git a/hw3/Hamiltonian.jl b/hw3/Hamiltonian.jl
new file mode 100644
index 0000000..bc8f31a
--- /dev/null
+++ b/hw3/Hamiltonian.jl
@@ -0,0 +1,175 @@
+using DifferentialEquations
+using Plots
+using Polyhedra # for finding the volume of the phase-space ensemble
+import CDDLib
+lib = CDDLib.Library()
+
+g = 9.8
+sim_time = 100.0
+
+function tendency!(du, u, p, t)
+    q_1, q_2, p_1, p_2 = u
+    g = p
+
+    denominator = 1.0 + (sin(q_1 - q_2)*sin(q_1 - q_2))
+
+    dq_1 = (p_1 - p_2 * cos(q_1 - q_2)) / denominator
+    dq_2 = (-p_1*cos(q_1 - q_2) + 2*p_2) / denominator
+
+    k_1 = p_1*p_2*sin(q_1 - q_2) / denominator
+    k_2 = sin(2*(q_1 - q_2)) * (p_1*p_1 + 2*p_2*p_2 - 2*p_1*p_2*cos(q_1 - q_2)) / (2*denominator*denominator)
+
+    dp_1 = -2 * sin(q_1) - k_1 + k_2
+    dp_2 = - sin(q_2) + k_1 - k_2 
+    
+    du[1] = dq_1
+    du[2] = dq_2
+    du[3] = dp_1
+    du[4] = dp_2
+end
+
+function energy(u, p)
+    q_1, q_2, p_1, p_2 = u
+    g = p
+
+    T = 0.5 * (p_1*p_1 + p_2*p_2)
+    V = - 2 * cos(q_1) - cos(q_2)
+    return T + V
+end
+
+function build_ensemble(q_start, q_end, p_start, p_end, n)
+    q_1s = []
+    q_2s = []
+    p_1s = []
+    p_2s = []
+    for i in 1:n
+        q_1 = q_start + (q_end - q_start) * rand()
+        q_2 = q_start + (q_end - q_start) * rand()
+        p_1 = p_start + (p_end - p_start) * rand()
+        p_2 = p_start + (p_end - p_start) * rand()
+        push!(q_1s, q_1)
+        push!(q_2s, q_2)
+        push!(p_1s, p_1)
+        push!(p_2s, p_2)
+    end
+
+    return [q_1s, q_2s, p_1s, p_2s]
+end
+
+# take a list and reduce theta to the interval [-π, π]
+function clean_θ(θ::Vector{Float64})
+    rθ = []
+    for i in 1:length(θ)
+        tmp = θ[i] % (2 * π)
+        if tmp > π
+            tmp = tmp - 2 * π
+        elseif tmp < -π
+            tmp = tmp + 2 * π
+        end
+        push!(rθ, tmp)
+    end
+    return rθ
+end
+
+function run_simulation(u_0, tspan, p)
+    prob = ODEProblem(tendency!, u_0, tspan, p)
+    sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # control simulation
+    sol[1, :] = clean_θ(sol[1, :])
+    sol[2, :] = clean_θ(sol[2, :])
+    return sol
+end
+
+function return_convex_hull(ensemble, library)
+    # plot the first phase space
+    vertices = zeros(length(ensemble[1]), 4)
+    for i in 1:length(ensemble[1])
+        vertices[i, 1] = ensemble[1][i]
+        vertices[i, 2] = ensemble[2][i]
+        vertices[i, 3] = ensemble[3][i]
+        vertices[i, 4] = ensemble[4][i]
+    end
+    v_1 = vrep(vertices)
+    p_1 = polyhedron(v_1, library)
+    vol_1 = Polyhedra.volume(p_1)
+
+    return p_1, vol_1
+end
+
+# returns the new ensemble after t_span time
+function simulate_ensembles(ensemble, tspan, p)
+    q_1s = []
+    q_2s = []
+    p_1s = []
+    p_2s = []
+    for i in 1:length(ensemble[1])
+        u_0 = [ensemble[1][i], ensemble[2][i], ensemble[3][i], ensemble[4][i]]
+        println("init cond" ,u_0)
+        sol = run_simulation(u_0, tspan, p)
+        q_1s = sol[1, :]
+        q_2s = sol[2, :]
+        p_1s = sol[3, :]
+        p_2s = sol[4, :]
+
+        push!(q_1s, q_1s[end])
+        push!(q_2s, q_2s[end])
+        push!(p_1s, p_1s[end])
+        push!(p_2s, p_2s[end])
+    end
+
+    return [q_1s, q_2s, p_1s, p_2s]
+end
+
+function calculate_ensemble_volume(ensemble)
+    dΓ = 1;
+    for i in 1:length(ensemble[1])
+        dΓ *= (ensemble[1][i]) * (ensemble[2][i]) * (ensemble[3][i]) * ensemble[4][i]
+    end
+    return dΓ
+end
+
+# PROBLEM 1 -> show that two systems with the same intial energy can have different chaotic results
+# u_1_0 = [0, 0, 6.26, 0.0]
+# u_2_0 = [0.0, 0.0, 0.0, 6.26]
+
+# # print the initial energies
+# println("energy 1 = ", energy(u_1_0, p))
+# println("energy 2 = ", energy(u_2_0, p))
+
+# sol_1 = run_simulation(u_1_0, tspan, p)
+# sol_2 = run_simulation(u_2_0, tspan, p)
+
+# # plot the phase space
+# p1 = scatter(sol_1[1, :], sol_1[3, :], label="q_1", xlabel="q_1 (radians)", ylabel="p_1 (momentum)", title="Phase Space i=1 for Pendulum 1", legend=false, color=:blue, msw=0, ms=.75)
+# p2 = scatter(sol_2[1, :], sol_2[3, :], label="q_2", xlabel="q_1 (radians)", ylabel="p_1 (momentum)", title="Phase Space i=1 for Pendulum 2", legend=false, color=:red, msw=0, ms=.75)
+
+# # plot the other phase space
+# p3 = scatter(sol_1[2, :], sol_1[4, :], label="q_1", xlabel="q_2 (radians)", ylabel="p_2 (momentum)", title="Phase Space i=2 for Pendulum 1", legend=false, color=:blue, msw=0, ms=.75)
+# p4 = scatter(sol_2[2, :], sol_2[4, :], label="q_2", xlabel="q_2 (radians)", ylabel="p_2 (momentum)", title="Phase Space i=2 for Pendulum 2", legend=false, color=:red, msw=0, ms=.75)
+
+# # plot the trajectories
+# p5 = scatter(sol_1[1, :], sol_1[2, :], label="q_1 (radians)", xlabel="q_1 (radians)", ylabel="q_2 (radians)", title="Relative Trajectory for Pend1", legend=false, color=:blue, msw=0, ms=.75)
+# p6 = scatter(sol_2[1, :], sol_2[2, :], label="q_2 (radians)", xlabel="q_1 (radians)", ylabel="q_2 (radians)", title="Relative Trajectory for Pend2", legend=false, color=:red, msw=0, ms=.75)
+
+# # plot the trajectories over time
+# p7 = plot(sol_1.t, sol_1[1, :], label="pendulum 1", xlabel="time (s)", ylabel="q_1 (radians)", title="q_1 vs time", legend=:topright)
+# plot!(p7, sol_2.t, sol_2[1, :], label="pendulum 2", xlabel="time (s)", ylabel="q_1 (radians)", title="q_1 vs time")
+
+# p8 = plot(sol_1.t, sol_1[2, :], label="pendulum 1", xlabel="time (s)", ylabel="q_2 (radians)", title="q_2 vs time", legend=:topright)
+# plot!(p8, sol_2.t, sol_2[2, :], label="pendulum 2", xlabel="time (s)", ylabel="q_2 (radians)", title="q_2 vs time")
+
+# plot(p1, p2, p3, p4, p5, p6, p7, p8, layout=(4, 2), size=(1000, 1000))
+
+# savefig("hw3/double_pendulum.png")
+
+# PROBLEM 2 -> show that the volume of the phase space is conserved
+
+# plot the initial ensemble in each phase space
+ensemble = build_ensemble(0.0, pi, -1.0, 2.0, 10)
+println("4d volume before = ", return_convex_hull(ensemble, lib)[2])
+
+# simulate the ensemble
+new_ensemble = simulate_ensembles(ensemble, tspan, p)
+println("4d volume after = ", return_convex_hull(new_ensemble, lib)[2])
+
+# println("energy before = ", energy(ensemble, p))
+# println("volume before = ", calculate_ensemble_volume(ensemble))