From c0c5ad459b4fe320d65db276b88534fb7c75b61c Mon Sep 17 00:00:00 2001 From: sotech117 Date: Mon, 26 Feb 2024 20:34:11 -0500 Subject: finish hw3 code and start hw4. add stencil code via upload --- .DS_Store | Bin 6148 -> 6148 bytes hw3/3-12.jl | 113 ++++++++++++++++++++++++++++++ hw3/3-13.jl | 3 + hw3/DrivenPendulum.jl | 110 ----------------------------- hw3/Hamiltonian.jl | 181 +++++++++++++++++++++++++----------------------- hw3/Lorenz.jl | 13 ++-- hw3/bigdata.png | Bin 0 -> 69727 bytes hw3/double_pendulum.png | Bin 672371 -> 604461 bytes hw3/ensemble.png | Bin 0 -> 133646 bytes hw3/l2.jl | 3 + hw4/4-10.jl | 0 hw4/SolarSystem1.jl | 83 ++++++++++++++++++++++ hw4/SolarSystem2.jl | 179 +++++++++++++++++++++++++++++++++++++++++++++++ 13 files changed, 480 insertions(+), 205 deletions(-) create mode 100644 hw3/3-12.jl delete mode 100644 hw3/DrivenPendulum.jl create mode 100644 hw3/bigdata.png create mode 100644 hw3/ensemble.png create mode 100644 hw4/4-10.jl create mode 100644 hw4/SolarSystem1.jl create mode 100644 hw4/SolarSystem2.jl diff --git a/.DS_Store b/.DS_Store index f47adcd..ada1a6c 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/hw3/3-12.jl b/hw3/3-12.jl new file mode 100644 index 0000000..c83e4c3 --- /dev/null +++ b/hw3/3-12.jl @@ -0,0 +1,113 @@ +# FOR PROBLEM 3.12 +# author: sotech117 + +#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia + +# Simulate driven pendulum to find chaotic regime + +using Plots # for plotting trajectory +using DifferentialEquations # for solving ODEs + +ω0 = 1.0 # ω0^2 = g/l +β = 0.5 # β = friction +f = 1.2 # forcing amplitude +ω = .66667 # forcing frequency +param = (ω0, β, f, ω) # parameters of anharmonic oscillator + +function tendency!(dθp::Vector{Float64}, θp::Vector{Float64}, param, t::Float64) + + (θ, p) = θp # 2d phase space + (dθ, dp) = dθp # 2d phase space derviatives + + (ω0, β, f, ω) = param + + a = -ω0^2 * sin(θ) - β * dθ + f * forcing(t, ω) # acceleration with m = 1 + + dθp[1] = p + dθp[2] = a +end + +function forcing(t::Float64, ω::Float64) + + return sin(ω * t) + +end + +function energy(θp::Vector{Float64}, param) + + (θ, p) = θp + + (ω0, β, f, ω) = param + + pe = ω0^2 * (1.0 - cos(θ)) + ke = 0.5 * p^2 + + return pe + ke + +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 get_poincare_sections(sample_θ, sample_p, sample_t, Ω_d, ϵ::Float64, phase_shift=0.0::Float64) + n = 0 + + poincare_θ = [] + poincare_p = [] + + for i in 1:length(sample_θ) + if abs(sample_t[i] * Ω_d - (2 * π * n + phase_shift)) < ϵ / 2 + push!(poincare_θ, sample_θ[i]) + push!(poincare_p, sample_p[i]) + n += 1 + end + end + + return (poincare_θ, poincare_p) +end + +θ0 = 0.2 # initial position in meters +p0 = 0.0 # initial velocity in m/s +θp0 = [θ0, p0] # initial condition in phase space +t_final = 1000.0 # final time of simulation + +tspan = (0.0, t_final) # span of time to simulate + +prob = ODEProblem(tendency!, θp0, tspan, param) # specify ODE +sol = solve(prob, Tsit5(), reltol=1e-12, abstol=1e-12) # solve using Tsit5 algorithm to specified accuracy + +sample_times = sol.t +println("\n\t Results") +println("final time = ", sample_times[end]) +println("Initial energy = ", energy(sol[:,1], param)) +println("Final energy = ", energy(sol[:, end], param)) + +(ω0, β, f, ω) = param + +# Plot of position vs. time +# θt = plot(sample_times, [sol[1, :], f * forcing.(sample_times, ω)], xlabel = "t", ylabel = "θ(t)", legend = false, title = "θ vs. t") + +# Phase space plot +cleaned = clean_θ(sol[1, :]) +θp = scatter(cleaned, sol[2, :], xlabel = "θ (radians)", ylabel = "ω (radians/s)", legend = false, title = "Phase Space Plot", mc=:black, ms=.35, ma=1) + + +# plot the poincare sections +(poincare_θ, pointcare_p) = get_poincare_sections(cleaned, sol[2, :], sol.t, ω, 0.1) +s1 = scatter(poincare_θ, pointcare_p, xlabel = "θ (radians)", ylabel = "ω (radians/s)", label="2nπ", title = "Poincare Sections", mc=:red, ms=2, ma=0.75) +s2 = scatter(get_poincare_sections(cleaned, sol[2, :], sol.t, ω, 0.1, π / 2.0), mc=:blue, ms=2, ma=0.75, label="2nπ + π/2", title="Poincare Sections", xlabel = "θ (radians)", ylabel = "ω (radians/s)", legend=:bottomleft) +s3 = scatter(get_poincare_sections(cleaned, sol[2, :], sol.t, ω, 0.1, π / 4.0), mc=:green, ms=2, ma=0.75, label="2nπ + π/4", title="Poincare Sections", xlabel = "θ (radians)", ylabel = "ω (radians/s)") + +plot(θp, s1, s2, s3) \ No newline at end of file diff --git a/hw3/3-13.jl b/hw3/3-13.jl index aacfa44..dc8dec6 100644 --- a/hw3/3-13.jl +++ b/hw3/3-13.jl @@ -1,3 +1,6 @@ +# FOR PROBLEM 3.13 +# author: sotech117 + #!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia # Simulate driven pendulum to find chaotic regime diff --git a/hw3/DrivenPendulum.jl b/hw3/DrivenPendulum.jl deleted file mode 100644 index 53a1af9..0000000 --- a/hw3/DrivenPendulum.jl +++ /dev/null @@ -1,110 +0,0 @@ -#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia - -# Simulate driven pendulum to find chaotic regime - -using Plots # for plotting trajectory -using DifferentialEquations # for solving ODEs - -ω0 = 1.0 # ω0^2 = g/l -β = 0.5 # β = friction -f = 1.2 # forcing amplitude -ω = .66667 # forcing frequency -param = (ω0, β, f, ω) # parameters of anharmonic oscillator - -function tendency!(dθp::Vector{Float64}, θp::Vector{Float64}, param, t::Float64) - - (θ, p) = θp # 2d phase space - (dθ, dp) = dθp # 2d phase space derviatives - - (ω0, β, f, ω) = param - - a = -ω0^2 * sin(θ) - β * dθ + f * forcing(t, ω) # acceleration with m = 1 - - dθp[1] = p - dθp[2] = a -end - -function forcing(t::Float64, ω::Float64) - - return sin(ω * t) - -end - -function energy(θp::Vector{Float64}, param) - - (θ, p) = θp - - (ω0, β, f, ω) = param - - pe = ω0^2 * (1.0 - cos(θ)) - ke = 0.5 * p^2 - - return pe + ke - -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 get_poincare_sections(sample_θ, sample_p, sample_t, Ω_d, ϵ::Float64, phase_shift=0.0::Float64) - n = 0 - - poincare_θ = [] - poincare_p = [] - - for i in 1:length(sample_θ) - if abs(sample_t[i] * Ω_d - (2 * π * n + phase_shift)) < ϵ / 2 - push!(poincare_θ, sample_θ[i]) - push!(poincare_p, sample_p[i]) - n += 1 - end - end - - return (poincare_θ, poincare_p) -end - -θ0 = 0.2 # initial position in meters -p0 = 0.0 # initial velocity in m/s -θp0 = [θ0, p0] # initial condition in phase space -t_final = 1000.0 # final time of simulation - -tspan = (0.0, t_final) # span of time to simulate - -prob = ODEProblem(tendency!, θp0, tspan, param) # specify ODE -sol = solve(prob, Tsit5(), reltol=1e-12, abstol=1e-12) # solve using Tsit5 algorithm to specified accuracy - -sample_times = sol.t -println("\n\t Results") -println("final time = ", sample_times[end]) -println("Initial energy = ", energy(sol[:,1], param)) -println("Final energy = ", energy(sol[:, end], param)) - -(ω0, β, f, ω) = param - -# Plot of position vs. time -# θt = plot(sample_times, [sol[1, :], f * forcing.(sample_times, ω)], xlabel = "t", ylabel = "θ(t)", legend = false, title = "θ vs. t") - -# Phase space plot -cleaned = clean_θ(sol[1, :]) -θp = scatter(cleaned, sol[2, :], xlabel = "θ (radians)", ylabel = "ω (radians/s)", legend = false, title = "Phase Space Plot", mc=:black, ms=.35, ma=1) - - -# plot the poincare sections -(poincare_θ, pointcare_p) = get_poincare_sections(cleaned, sol[2, :], sol.t, ω, 0.1) -s1 = scatter(poincare_θ, pointcare_p, xlabel = "θ (radians)", ylabel = "ω (radians/s)", label="2nπ", title = "Poincare Sections", mc=:red, ms=2, ma=0.75) -s2 = scatter(get_poincare_sections(cleaned, sol[2, :], sol.t, ω, 0.1, π / 2.0), mc=:blue, ms=2, ma=0.75, label="2nπ + π/2", title="Poincare Sections", xlabel = "θ (radians)", ylabel = "ω (radians/s)", legend=:bottomleft) -s3 = scatter(get_poincare_sections(cleaned, sol[2, :], sol.t, ω, 0.1, π / 4.0), mc=:green, ms=2, ma=0.75, label="2nπ + π/4", title="Poincare Sections", xlabel = "θ (radians)", ylabel = "ω (radians/s)") - -plot(θp, s1, s2, s3) \ No newline at end of file diff --git a/hw3/Hamiltonian.jl b/hw3/Hamiltonian.jl index bc8f31a..7cefe55 100644 --- a/hw3/Hamiltonian.jl +++ b/hw3/Hamiltonian.jl @@ -1,11 +1,12 @@ +# FOR PHYS2600 HW3 +# author: sotech117 + 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 +t_span = (0.0, 1000.0) +p = g function tendency!(du, u, p, t) q_1, q_2, p_1, p_2 = u @@ -17,10 +18,10 @@ function tendency!(du, u, p, t) 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) + k_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 + dp_1 = -2 * sin(q_1) - k_1 + k_2 * sin(2*(q_1 - q_2)) + dp_2 = - sin(q_2) + k_1 - k_2 * sin(2*(q_1 - q_2)) du[1] = dq_1 du[2] = dq_2 @@ -38,21 +39,16 @@ function energy(u, p) end function build_ensemble(q_start, q_end, p_start, p_end, n) - q_1s = [] - q_2s = [] - p_1s = [] - p_2s = [] + q_1s = zeros(n) + q_2s = zeros(n) + p_1s = zeros(n) + p_2s = zeros(n) 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) + q_1s[i] = q_start + (q_end - q_start) * rand() + q_2s[i] = q_start + (q_end - q_start) * rand() + p_1s[i] = p_start + (p_end - p_start) * rand() + p_2s[i] = p_start + (p_end - p_start) * rand() end - return [q_1s, q_2s, p_1s, p_2s] end @@ -79,97 +75,110 @@ function run_simulation(u_0, tspan, p) 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 = [] + q_1s = zeros(length(ensemble[1])) + q_2s = zeros(length(ensemble[2])) + p_1s = zeros(length(ensemble[3])) + p_2s = zeros(length(ensemble[4])) + println("length ensemble = ", length(ensemble[1])) 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 + q_1s[i] = sol[1, end] + q_2s[i] = sol[2, end] + p_1s[i] = sol[3, end] + p_2s[i] = sol[4, end] + end + println("q1_s", length(q_1s)) return [q_1s, q_2s, p_1s, p_2s] end -function calculate_ensemble_volume(ensemble) - dΓ = 1; +function sum_all_positions_and_momenta(ensemble) + q_1 = 0 + q_2 = 0 + p_1 = 0 + p_2 = 0 for i in 1:length(ensemble[1]) - dΓ *= (ensemble[1][i]) * (ensemble[2][i]) * (ensemble[3][i]) * ensemble[4][i] + q_1 += ensemble[1][i] + q_2 += ensemble[2][i] + p_1 += ensemble[3][i] + p_2 += ensemble[4][i] end - return dΓ + return [q_1, q_2, p_1, p_2] end +function estimate_phase_space_volume(ensemble, n) + # make 100000 random points in the 4d phase space from [-pi, pi] x [-pi, pi] x [-3, 3] x [-3, 3] + points = build_ensemble(-pi, pi, -3, 3, 100000) + + # iterate over the points and find the number of points that hit the ensemble + count = 0 + for i in 1:length(points[1]) + for j in 1:length(ensemble[1]) + if abs(points[1][i] - ensemble[1][j]) < n && abs(points[2][i] - ensemble[2][j]) < n && abs(points[3][i] - ensemble[3][j]) < n && abs(points[4][i] - ensemble[4][j]) < n + count += 1 + break + end + end + end + + # return the volume + return (count / 100000) # note units don't matter, just comparing +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] +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)) +# 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) +sol_1 = run_simulation(u_1_0, t_span, p) +sol_2 = run_simulation(u_2_0, t_span, 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 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 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 +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") +# 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") +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)) +plt = plot(p1, p2, p3, p4, p5, p6, p7, p8, layout=(4, 2), size=(1000, 1000)) -# savefig("hw3/double_pendulum.png") +savefig(plt, "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]) - +ensemble = build_ensemble(-1, 1, -1, 1, 1000) +# plot these points +s1 = scatter(ensemble[1], ensemble[3], label="t=0", xlabel="q_1 (radians)", ylabel="p_1 (momentum)", title="Phase Space (i=1)", color=:blue, msw=.5) +s2 = scatter(ensemble[2], ensemble[4], label="t=0", xlabel="q_2 (radians)", ylabel="p_2 (momentum)", title="Phase Space (i=2)", color=:red, msw=.5) # simulate the ensemble -new_ensemble = simulate_ensembles(ensemble, tspan, p) -println("4d volume after = ", return_convex_hull(new_ensemble, lib)[2]) +new_ensemble = simulate_ensembles(ensemble, t_span, p) +# plot these points +scatter!(s1, new_ensemble[1], new_ensemble[3], label="t=1000", xlabel="q_1 (radians)", ylabel="p_1 (momentum)", title="Phase Space (i=1)", color=:green, msw=.5) +scatter!(s2, new_ensemble[2], new_ensemble[4], label="t=1000", xlabel="q_2 (radians)", ylabel="p_2 (momentum)", title="Phase Space (i=2)", color=:green, msw=.5) + +# print the volumes +println("\n\nOriginal Volume:\t", estimate_phase_space_volume(ensemble, .1)) +println("Final Volume:\t", estimate_phase_space_volume(new_ensemble, .1)) + +plt_2 = plot(s1, s2) +savefig(plt_2, "hw3/ensemble.png") -# println("energy before = ", energy(ensemble, p)) -# println("volume before = ", calculate_ensemble_volume(ensemble)) diff --git a/hw3/Lorenz.jl b/hw3/Lorenz.jl index 4a7fa0e..b1dd8f0 100644 --- a/hw3/Lorenz.jl +++ b/hw3/Lorenz.jl @@ -1,3 +1,6 @@ +# FOR PROBLEM 3.25 +# author: sotech117 + using DifferentialEquations using Plots @@ -25,8 +28,7 @@ r_maxes = [] z_maxes = [] r_mins = [] z_mins = [] -let - r_diverge_1 = -1 +let for i in 1:r_steps r = r_values[i] p = [10.0, r, 8/3] # parameters of the Lorentz 63 system @@ -49,12 +51,6 @@ let push!(r_mins, r_values[i]) push!(z_mins, z_values[j]) end - - # calcualte dz - # if r_diverge_1 < 0 && length(z_mins) > 1 && length(z_maxes) > 1 && abs(z_mins[end] - z_maxes[end]) > 20.0 - # r_diverge_1 = r - # println("Divergence starts at r = $r") - # end end if i % 500 == 0 @@ -70,6 +66,5 @@ let vline!([122], label="Critical r value ≈ 122", lc=:orange, lw=1.5, ls=:dash) vline!([146], label="Critical r value ≈ 146", lc=:purple, lw=1.5, ls=:dash) - savefig("hw3/test11.png") end \ No newline at end of file diff --git a/hw3/bigdata.png b/hw3/bigdata.png new file mode 100644 index 0000000..589f43a Binary files /dev/null and b/hw3/bigdata.png differ diff --git a/hw3/double_pendulum.png b/hw3/double_pendulum.png index e95959a..1024c62 100644 Binary files a/hw3/double_pendulum.png and b/hw3/double_pendulum.png differ diff --git a/hw3/ensemble.png b/hw3/ensemble.png new file mode 100644 index 0000000..dd9fb70 Binary files /dev/null and b/hw3/ensemble.png differ diff --git a/hw3/l2.jl b/hw3/l2.jl index 64a9cc0..2e1519c 100644 --- a/hw3/l2.jl +++ b/hw3/l2.jl @@ -1,3 +1,6 @@ +# FOR PROBLEM 3.14 +# author: sotech117 + using DifferentialEquations using Plots diff --git a/hw4/4-10.jl b/hw4/4-10.jl new file mode 100644 index 0000000..e69de29 diff --git a/hw4/SolarSystem1.jl b/hw4/SolarSystem1.jl new file mode 100644 index 0000000..e62ee4f --- /dev/null +++ b/hw4/SolarSystem1.jl @@ -0,0 +1,83 @@ +#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia + +# Simulate planet around stationary star + +using Plots # for plotting trajectory +using DifferentialEquations # for solving ODEs +using LinearAlgebra +#using Unitful +#using UnitfulAstro # astrophysical units + +G = 4.0*pi^2 # time scale = year and length scale = AU +δ = 0.0 # deviation from inverse-square law + +param = (G, δ) # parameters of force law + +function tendency!(drp, rp, param, t) + + # 6d phase space + r = rp[1:3] + p = rp[4:6] + + (G, δ) = param + + r2 = dot(r, r) + + a = -G * r * r2^(-1.5-δ) # acceleration with m = 1 + + drp[1:3] = p[1:3] + drp[4:6] = a[1:3] + +end + +function energy(rp, param) + + r = rp[1:3] + p = rp[4:6] + + (G, δ) = param + + r2 = dot(r, r) + + pe = -G * r2^(-0.5 - δ)/(1.0 + 2.0 * δ) + ke = 0.5 * dot(p, p) + + return pe + ke + +end + +function angularMomentum(rp) + + r = rp[1:3] + p = rp[4:6] + + return cross(r, p) + +end + +r0 = [1.0, 0.0, 0.0] # initial position in AU +p0 = [0.0, 1.0*pi, 0.0] # initial velocity in AU / year +rp0 = vcat(r0, p0) # initial condition in phase space +t_final = 10.0 # final time of simulation + +tspan = (0.0, t_final) # span of time to simulate + + +prob = ODEProblem(tendency!, rp0, tspan, param) # specify ODE +sol = solve(prob, Tsit5(), reltol=1e-12, abstol=1e-12) # solve using Tsit5 algorithm to specified accuracy + +sample_times = sol.t +println("\n\t Results") +println("final time = ", sample_times[end]) +println("Initial energy = ", energy(sol[:,1], param)) +println("Final energy = ", energy(sol[:, end], param)) +println("Initial angular momentum = ", angularMomentum(sol[:,1])) +println("Final angular momentum = ", angularMomentum(sol[:, end])) + +# Plot of position vs. time +xt = plot(sample_times, sol[1, :], xlabel = "t", ylabel = "x(t)", legend = false, title = "x vs. t") + +# Plot of orbit +xy = plot(sol[1, :], sol[2, :], xlabel = "x", ylabel = "y", legend = false, title = "Orbit", aspect_ratio=:equal) + +plot(xt, xy) diff --git a/hw4/SolarSystem2.jl b/hw4/SolarSystem2.jl new file mode 100644 index 0000000..b2e02b7 --- /dev/null +++ b/hw4/SolarSystem2.jl @@ -0,0 +1,179 @@ +#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia + +# Simulate a solar system + +using Plots # for plotting trajectory +using DifferentialEquations # for solving ODEs +using LinearAlgebra # for dot and cross products + +G = 4.0*pi^2 # time scale = year and length scale = AU + +mutable struct body + name::String # name of star or planet + m::Float64 # mass + r::Vector{Float64} # position vector + p::Vector{Float64} # momentum vector +end + +function angularMomentum(b::body) + r = b.r + p = b.p + return cross(r, p) +end + +function kineticEnergy(b::body) + p = b.p + m = b.m + return dot(p, p) / (2.0 * m) +end + +function rp(b::body) + return [b.r; b.p] +end + +function force(body1::body, body2::body) + r = body1.r - body2.r + rSquared = dot(r, r) + return -G * body1.m * body2.m * r * rSquared^(-1.5) +end + +function potentialEnergy(body1::body, body2::body) + r = body1.r - body2.r + rSquared = dot(r, r) + return -G * body1.m * body2.m * rSquared^(-0.5) +end + +mutable struct SolarSystem + bodies::Vector{body} + numberOfBodies::Int64 + phaseSpace::Matrix{Float64} # 6N dimensional phase space +end + +function SolarSystem() + + bodies = Vector{body}() + push!(bodies, body("Sun", 1.0, zeros(3), zeros(3))) + #push!(bodies, body("Earth", 1.0, [-1.0, 0.0, 0.0], [0.0, 1.0 * pi, 0.0])) + push!(bodies, body("Star", 1.0, [1.0, 0.0, 0.0], [0.0, sqrt(1.95) * 2.0 * pi, 0.0])) + #push!(bodies, body("Jupiter", 1.0, [3.0, 0.0, 0.0], [0.0, 0.25 * pi, 0.0])) + numberOfBodies = size(bodies)[1] + + phaseSpace = zeros(6, 0) + for b in bodies + phaseSpace = [phaseSpace rp(b)] + end + + return SolarSystem(bodies, numberOfBodies, phaseSpace) + +end + +function TotalAngularMomentum(s::SolarSystem) + L = zeros(3) + for b in s.bodies + L += angularMomentum(b) + end + + return L +end + +function TotalEnergy(s::SolarSystem) + ke = 0.0 + pe = 0.0 + + for body1 in s.bodies + ke += kineticEnergy(body1) + for body2 in s.bodies + if (body1 != body2) + pe += 0.5 * potentialEnergy(body1, body2) + end + end + end + + return pe + ke +end + +function ZeroOutLinearMomentum!(s::SolarSystem) + + totalLinearMomentum = zeros(3) + totalMass = 0.0 + for body in s.bodies + totalLinearMomentum += body.p + totalMass += body.m + end + + s.phaseSpace = zeros(6, 0) + for body in s.bodies + body.p -= body.m * totalLinearMomentum / totalMass + s.phaseSpace = [s.phaseSpace rp(body)] + end + + return nothing +end + +function tendency!(dps, ps, s::SolarSystem, t) + + i = 1 # update phase space of individual bodies + for b in s.bodies + b.r = ps[1:3, i] + b.p = ps[4:6, i] + i += 1 + end + + # find velocities of bodies and forces on them. O(N^2) computational cost + N = s.numberOfBodies + for i in 1:N + b1 = s.bodies[i] + dps[1:3, i] = b1.p / b1.m #dr/dt + dps[4:6, i] = zeros(3) + for j in 1:i-1 + b2 = s.bodies[j] + f = force(b1, b2) # call only once + dps[4:6, i] += f + dps[4:6, j] -= f # Newton's 3rd law + end + end + + return nothing +end + +s = SolarSystem() +ZeroOutLinearMomentum!(s) +println(typeof(s)) +println("Initial total energy = ", TotalEnergy(s)) +println("Initial total angular momentum = ", TotalAngularMomentum(s)) +println("Number of bodies = ", s.numberOfBodies) +for b in s.bodies + println("body name = ", b.name) +end + +t_final = 10.0 # final time of simulation +tspan = (0.0, t_final) # span of time to simulate +prob = ODEProblem(tendency!, s.phaseSpace, tspan, s) # specify ODE +sol = solve(prob, Tsit5(), reltol=1e-12, abstol=1e-12) # solve using Tsit5 algorithm to specified accuracy + +sample_times = sol.t +println("\n\t Results") +println("Final time = ", sample_times[end]) +println("Final total energy = ", TotalEnergy(s)) +println("Final total angular momentum = ", TotalAngularMomentum(s)) + +body1 = 1 +body2 = 2 +# Plot of position vs. time +xt = plot(sample_times, sol[1, body1, :], xlabel = "t", ylabel = "x(t)", legend = false, title = "x vs. t") + +# Plot of orbit +xy = plot([(sol[1, body1, :], sol[2, body1, :]), (sol[1, body2, :], sol[2, body2, :])], colors = ("yellow", "green"), xlabel = "x", ylabel = "y", legend = false, title = "Orbit") +#= +plot(xy, aspect_ratio=:equal) +=# + +samples = 1000 +interval = floor(Int,size(sol.t)[1] / samples) +animation = @animate for i=1:samples-1 + plot([(sol[1, body1, i*interval], sol[2, body1, i*interval]), (sol[1, body1, (i+1)*interval], sol[2, body1, (i+1)*interval])], + aspect_ratio=:equal, colors = ("yellow", "green"), xlabel = "x", ylabel = "y", legend = false, title = "Orbit", xlims=(-1, 1), ylims = (-1, 1)) +end + +gif(animation, fps=15) + -- cgit v1.2.3-70-g09d2