|
| 1 | +# # Steady state sweeps |
| 2 | + |
| 3 | +using HarmonicBalance, SteadyStateDiffEq, ModelingToolkit |
| 4 | +using BenchmarkTools, Plots, StaticArrays, OrdinaryDiffEq, LinearAlgebra |
| 5 | +using HarmonicBalance: OrderedDict |
| 6 | + |
| 7 | +@variables α ω ω0 F γ η t x(t); |
| 8 | + |
| 9 | +diff_eq = DifferentialEquation( |
| 10 | + d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) + η * d(x, t) * x^2 ~ F * cos(ω * t), x |
| 11 | +) |
| 12 | +add_harmonic!(diff_eq, x, ω) |
| 13 | +harmonic_eq = get_harmonic_equations(diff_eq) |
| 14 | + |
| 15 | +fixed = (ω0 => 1.0, γ => 1e-2, F => 0.02, α => 1.0, η => 0.3) |
| 16 | +ω_span = (0.8, 1.3); |
| 17 | +ω_range = range(ω_span..., 201); |
| 18 | +varied = ω => ω_range |
| 19 | + |
| 20 | +result_HB = get_steady_states(harmonic_eq, varied, fixed) |
| 21 | +plot(result_HB, "sqrt(u1^2+v1^2)") |
| 22 | + |
| 23 | +# ## Steady state sweep using `SteadyStateDiffEq.jl` |
| 24 | +param = OrderedDict(merge(Dict(fixed), Dict(ω => 1.1))) |
| 25 | +x0 = [0, 0.0]; |
| 26 | +prob_ss = SteadyStateProblem( |
| 27 | + harmonic_eq, x0, param; in_place=false, u0_constructor=x -> SVector(x...) |
| 28 | +) |
| 29 | + |
| 30 | +varied = 6 => ω_range |
| 31 | +result_ss = steady_state_sweep( |
| 32 | + prob_ss, DynamicSS(Rodas5()); varied, abstol=1e-5, reltol=1e-5 |
| 33 | +) |
| 34 | + |
| 35 | +plot(result_HB, "sqrt(u1^2+v1^2)") |
| 36 | +plot!(ω_range, norm.(result_ss); c=:gray, ls=:dash) |
| 37 | + |
| 38 | +# ## Adiabatic evolution |
| 39 | + |
| 40 | +timespan = (0.0, 50_000) |
| 41 | +sweep = AdiabaticSweep(ω => (0.8, 1.3), timespan) # linearly interpolate between two values at two times |
| 42 | +ode_problem = ODEProblem(harmonic_eq, fixed; u0=[0.01; 0.0], timespan, sweep) |
| 43 | +time_soln = solve(ode_problem, Tsit5(); saveat=250) |
| 44 | + |
| 45 | +plot(result_HB, "sqrt(u1^2+v1^2)") |
| 46 | +plot(time_soln.t, norm.(time_soln.u)) |
| 47 | + |
| 48 | +# ## using follow_branch |
| 49 | + |
| 50 | +followed_branch, Ys = follow_branch(1, result_HB; y="√(u1^2+v1^2)") |
| 51 | +Y_followed_gr = |
| 52 | + real.([Ys[param_idx][branch] for (param_idx, branch) in enumerate(followed_branch)]); |
| 53 | + |
| 54 | +plot(result_HB, "sqrt(u1^2+v1^2)") |
| 55 | +plot!(ω_range, Y_followed_gr; c=:gray, ls=:dash) |
| 56 | + |
| 57 | +# ## comparison |
| 58 | + |
| 59 | +@btime result_ss = steady_state_sweep( |
| 60 | + prob_ss, DynamicSS(Rodas5()); varied, abstol=1e-5, reltol=1e-5 |
| 61 | +) |
| 62 | + |
| 63 | +@btime time_soln = solve(ode_problem, Tsit5(); saveat=250) |
| 64 | + |
| 65 | +@btime begin |
| 66 | + followed_branch, Ys = follow_branch(1, result_HB; y="√(u1^2+v1^2)") |
| 67 | + Y_followed_gr = |
| 68 | + real.([Ys[param_idx][branch] for (param_idx, branch) in enumerate(followed_branch)]) |
| 69 | +end |
| 70 | + |
| 71 | +# Plotting them together gives: |
| 72 | +plot(ω_range, norm.(result_ss)) |
| 73 | +plot!(ω_range, norm.(time_soln.u)) |
| 74 | +plot!(ω_range, Y_followed_gr) |
0 commit comments