1 files changed, 71 insertions, 0 deletions
diff --git a/hw5/HarmonicOscillator.jl b/hw5/HarmonicOscillator.jl
new file mode 100644
index 0000000..adfaff8
--- /dev/null
+++ b/hw5/HarmonicOscillator.jl
@@ -0,0 +1,71 @@
+#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia
+
+# Simulate anharmonic oscillator that may be damped and driven
+
+using Plots # for plotting trajectory
+using DifferentialEquations # for solving ODEs
+
+ω0 = 1.0 # ω0^2 = k/m
+β = 0.0 # β = b/m = friction
+c = 10.0 # anharmonic term
+f = 0.3 # forcing amplitude
+ω = 2.0 # forcing frequency
+param = (ω0, β, c, f, ω) # parameters of anharmonic oscillator
+
+function tendency!(dxp::Vector{Float64}, xp::Vector{Float64}, param, t::Float64)
+   
+   (x, p) = xp # 2d phase space
+
+   (ω0, β, c, f, ω) = param
+
+   a = -ω0^2 * x - β * p - c * x^3 + f * forcing(t, ω) # acceleration with m = 1
+
+   dxp[1] = p
+   dxp[2] = a
+ 
+end
+
+function forcing(t::Float64, ω::Float64)
+
+   return cos(ω * t)
+
+end
+
+function energy(xp::Vector{Float64}, param)
+
+   (x, p) = xp
+
+   (ω0, β, c, f, ω) = param
+
+   pe = 0.5 * ω0^2 * x^2 + (c/4.0) * x^4
+   ke = 0.5 * p^2
+
+   return pe + ke
+
+end
+
+x0 = 0.0 # initial position in meters
+p0 = 0.0 # initial velocity in m/s
+xp0 = [x0, p0] # initial condition in phase space 
+t_final = 100.0 # final time of simulation
+
+tspan = (0.0, t_final) # span of time to simulate
+
+prob = ODEProblem(tendency!, xp0, tspan, param) # specify ODE
+sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # 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, β, c, f, ω) = param
+
+# Plot of position vs. time
+xt = plot(sample_times, [sol[1, :] f * forcing.(sample_times, ω)], xlabel = "t", ylabel = "x(t)", legend = false, title = "x vs. t")
+
+# Phase space plot
+xp = plot(sol[1, :], sol[2, :], xlabel = "x", ylabel = "p", legend = false, title = "phase space")
+
+plot(xt, xp)
+\ No newline at end of file