# FOR PROBLEM 3.13 # author: sotech117 #!/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)