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
72
73
74
75
76
77
78
79
80
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")
|
