1 files changed, 81 insertions, 0 deletions
diff --git a/hw3/t.jl b/hw3/t.jl
new file mode 100644
index 0000000..307d532
--- /dev/null
+++ b/hw3/t.jl
@@ -0,0 +1,81 @@
+#!/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
+
+θ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 = 150.0 # final time of simulation
+
+tspan = (0.0, t_final) # span of time to simulate
+
+prob1 = ODEProblem(tendency!, θp0, tspan, param) # specify ODE
+sol1 = solve(prob1, Tsit5(), reltol=1e-12, abstol=1e-12) # solve using Tsit5 algorithm to specified accuracy
+
+# do second prlblem with different initial conditions
+θ0 = 0.2001 # initial position in meters
+θp0 = [θ0, p0] # initial condition in phase space
+prob2 = ODEProblem(tendency!, θp0, tspan, param) # specify ODE
+sol2 = solve(prob2, Tsit5(), reltol=1e-12, abstol=1e-12) # solve using Tsit5 algorithm to specified accuracy
+
+function interpolate_Δθ(sol1, sol2)
+   Δθ = []
+   # find the time steps that are common to both solutions
+    common_times = range(0.00, stop=t_final, length=1000)
+    for t in common_times
+       # find the two time indices that are closest to t
+        i1 = argmin(abs.(sol1.t .- t))
+        i2 = argmin(abs.(sol2.t .- t))
+        push!(Δθ, abs(sol1[1, i1] - sol2[1, i2]))
+    end
+
+    return common_times, log.(Δθ)
+end
+
+(ω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")
+
+plot(interpolate_Δθ(sol1, sol2), xlabel = "t", ylabel = "Δθ(t)", legend = false, title = "Δθ vs. t")
+\ No newline at end of file