finish first 3 problems - codewise

1 files changed, 71 insertions, 0 deletions

diff --git a/hw3/3-13.jl b/hw3/3-13.jl
new file mode 100644
index 0000000..aacfa44
--- /dev/null
+++ b/hw3/3-13.jl
@@ -0,0 +1,71 @@
+#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia
+
+# Simulate driven pendulum to find chaotic regime
+
+using Plots # for plotting trajectory
+
+ω0 = 1.0 # ω0^2 = g/l
+β = 0.5 # β = friction
+f = 1.2 # forcing amplitude
+Ω_D = .66667 # forcing frequency
+param = (ω0, β, f, Ω_D) # parameters of anharmonic oscillator
+dt = 0.001 # time step
+t_final = 150.0 # final time of simulation
+
+function dynamics!(θ::Vector{Float64}, ω::Vector{Float64}, t::Vector{Float64}, Δt, param, steps)
+    (ω0, β, f, Ω_D) = param
+
+    for i in 1:steps
+        dω = -ω0^2 * sin(θ[end]) - β * ω[end] + f * forcing(t[end], Ω_D) # acceleration with m = 1
+        push!(ω, dω*dt + ω[end])
+
+        # euler-cromer method, using next ω
+        dθ = ω[end]
+        push!(θ, dθ*dt + θ[end])
+
+        push!(t, t[end] + Δt)
+    end
+end
+
+function forcing(t::Float64, ω::Float64)
+   return sin(ω * t)
+end
+
+function do_simulation(θ0, p0, t_final, Δt, param)
+    θs = [θ0]
+    ωs = [p0]
+    ts = [0.0]
+    steps = Int64(t_final/dt) # number of time steps
+    dynamics!(θs, ωs, ts, Δt, param, steps)
+    return (θs, ωs, ts)
+end
+
+(θ_1, _, t_1) = do_simulation(0.2, 0.0, t_final, dt, param)
+(θ_2, _, t_2) = do_simulation(0.2001, 0.0, t_final, dt, param)
+
+# absolute difference the two simulations thetas
+function abs_diff(θ1, θ2)
+    diff = zeros(length(θ1))
+    for i in 1:length(θ1)
+        diff[i] = abs(θ1[i] - θ2[i])
+    end
+    return diff
+end
+
+diff = abs_diff(θ_1, θ_2)
+plot(t_1, diff, yscale=:ln, xlabel="time (s)", ylabel="Δθ (radians)", title="Δθ vs time", legend=:bottomright, lw=2, label="|θ_1 - θ_2|")
+
+function linear_regression(x, y)
+    n = length(x)
+    x̄ = sum(x) / n
+    ȳ = sum(y) / n
+    a = sum((x .- x̄) .* (y .- ȳ)) / sum((x .- x̄).^2)
+    b = ȳ - a * x̄
+    return (a, b)
+end
+
+# linear regression of the log of the difference
+(a, b) = linear_regression(t_1, log.(diff))
+a = round(a, digits=4)
+b = round(b, digits=4)
+plot!(t_1, exp.(a*t_1 .+ b), label="exp($a*t + $b)", lw=1, color=:red, linestyle=:dash)