1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
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)
|
