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 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)
|
