diff options
Diffstat (limited to 'hw3/Hamiltonian.jl')
| -rw-r--r-- | hw3/Hamiltonian.jl | 175 |
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)) |