From 80ef2a117a0bca5bb6dbf025660d6c924f815d54 Mon Sep 17 00:00:00 2001
From: sotech117 <michael_foiani@brown.edu>
Date: Wed, 21 Feb 2024 00:48:31 -0500
Subject: pull files from canvas to build off on

---
 hw3/DrivenPendulum.jl | 70 +++++++++++++++++++++++++++++++++++++++++++++++++++
 hw3/Lorenz.jl         | 32 +++++++++++++++++++++++
 2 files changed, 102 insertions(+)
 create mode 100644 hw3/DrivenPendulum.jl
 create mode 100644 hw3/Lorenz.jl

diff --git a/hw3/DrivenPendulum.jl b/hw3/DrivenPendulum.jl
new file mode 100644
index 0000000..8fe6a14
--- /dev/null
+++ b/hw3/DrivenPendulum.jl
@@ -0,0 +1,70 @@
+#!/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.0001 # β = friction
+f = 0.5 # forcing amplitude
+ω = 1.01 # 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
+
+   (ω0, β, f, ω) = param
+
+   a = -ω0^2 * sin(θ) - β * p + f * forcing(t, ω) # acceleration with m = 1
+
+   dθp[1] = p
+   dθp[2] = a
+ 
+end
+
+function forcing(t::Float64, ω::Float64)
+
+   return cos(ω * 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.0 # initial position in meters
+p0 = 0.0 # initial velocity in m/s
+θp0 = [θ0, p0] # initial condition in phase space 
+t_final = 10000.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
+θp = plot(sin.(sol[1, :]), sol[2, :], xlabel = "θ", ylabel = "p", legend = false, title = "phase space")
+
+plot(θt, θp)
\ No newline at end of file
diff --git a/hw3/Lorenz.jl b/hw3/Lorenz.jl
new file mode 100644
index 0000000..0144c92
--- /dev/null
+++ b/hw3/Lorenz.jl
@@ -0,0 +1,32 @@
+using DifferentialEquations
+using Plots
+
+# Simulate Lorenz 63 system and investigate sensitivity to initial conditions
+
+function tendency!(du, u, p, t)
+    x,y,z = u
+    σ,ρ,β = p
+    
+    du[1] = dx = σ*(y-x)
+    du[2] = dy = x*(ρ-z) - y
+    du[3] = dz = x*y - β*z
+end
+
+p = [10.0, 28.0, 8/3] # parameters of the Lorentz 63 system
+tspan = (0.0, 1000.0)
+
+u0 = [1.0, 0.0, 0.0]
+δu0 = [0.0, 0.001, 0.0] # small deviation in initial condition
+
+prob = ODEProblem(tendency!, u0, tspan, p)
+sol1 = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # control simulation
+
+
+prob = ODEProblem(tendency!, u0 + δu0, tspan, p)
+sol2 = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # perturbed initial condition
+
+plot(sol1, idxs = (1, 2, 3)) # 3d plot
+
+# plot(sol1, idxs = (1, 2))
+
+plot([sol1[1, :], sol2[1, :]], [sol1[2, :], sol2[2, :]])
\ No newline at end of file
-- 
cgit v1.2.3-70-g09d2