Example D2: Bifurcation Analysis of the Forced Response of a Van der Pol Oscillator
This example is designed to demonstrate stability analysis and branch switching across a (secondary) Hopf bifurcation that occurs in a forced Van der Pol oscillator.
The dynamical system that will be studied here is:
\[\ddot{x} - c \dot{x} + k x + \eta x^2 \dot{x} = F\cos\Omega t.\]
This is an SDoF nonlinear oscillator with negative viscous damping. The nonlinearity provides the dissipation that saturates the response. Classically known as the Van der Pol oscillator, this is a useful system for learning about limit cycles. In the considered version, we also have an external periodic excitation.
The system will have a stable periodic response close to resonance, which will lose its stability as we move away from it.
Preamble: Load Packages
using GLMakie
using LinearAlgebra
using SparseArrays
using NonlinearSolve
using DSP
using Revise
using juliajim.HARMONIC
using juliajim.CONTINUATION
using juliajim.MDOFUTILS
using juliajim.NLDYNSystem Setup
Here we setup the system with parameters.
pars = (c=0.02, k=4., F=1., eta=0.1);
mdl = MDOFGEN(1.0, -pars.c, pars.k);
fnl = (t,u,ud) -> return pars.eta.*u.^2 .*ud, 2pars.eta.*u.*ud, pars.eta.*u.^2;
mdl = ADDNL(mdl, :Inst, fnl, 1.0);Setup HB
h = 0:5;
N = 128;
t = range(0, 2π, N+1)[1:N];
Nhc = NHC(h);
inds0, hinds, zinds, rinds, iinds = HINDS(mdl.Ndofs, h);
Fl = zeros(Nhc*mdl.Ndofs, 1);
Fl[rinds[1]] = pars.F;
Om0 = 0.1;
Om1 = 2.0;
E, dEdw = HARMONICSTIFFNESS(mdl.M, mdl.C, mdl.K, Om0, h);
Uw0 = [E\Fl; Om0];Setup HB Residue, Get First Point
fun = NonlinearFunction((r,u,p)->HBRESFUN!([u;p], mdl, Fl, h, N; R=r),
jac=(J,u,p)->HBRESFUN!([u;p], mdl, Fl, h, N; dRdU=J),
paramjac=(Jp,u,p)->HBRESFUN!([u;p], mdl, Fl, h, N; dRdw=Jp));
prob = NonlinearProblem(fun, Uw0[1:end-1], Uw0[end]);
sol = solve(prob, show_trace=Val(true));
Algorithm: NewtonRaphson(
descent = NewtonDescent(),
autodiff = AutoForwardDiff(),
vjp_autodiff = AutoFiniteDiff(
fdtype = Val{:forward}(),
fdjtype = Val{:forward}(),
fdhtype = Val{:hcentral}(),
dir = true
),
jvp_autodiff = AutoForwardDiff(),
concrete_jac = Val{false}()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 3.93569199e-05 0.00000000e+00
1 9.39503268e-13 1.40930725e-05
2 1.11022302e-16 2.66616371e-13
Final 1.11022302e-16
----------------------Continuation
Just like in other examples we call the CONTINUATE utility to obtain the periodic forced response.
Om0 = 0.1;
Om1 = 4.0;
dOm = 0.4;
cpars = (parm=:arclength, nmax=100, save_jacs=true);
sols, _, _, _, _ = CONTINUATE(Uw0[1:end-1], fun, [Om0, Om1], dOm; cpars...);
Algorithm: NewtonRaphson(
descent = NewtonDescent(),
autodiff = AutoForwardDiff(),
vjp_autodiff = AutoFiniteDiff(
fdtype = Val{:forward}(),
fdjtype = Val{:forward}(),
fdhtype = Val{:hcentral}(),
dir = true
),
jvp_autodiff = AutoForwardDiff(),
concrete_jac = Val{false}()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 3.93569199e-05 0.00000000e+00
1 9.39503268e-13 1.40930725e-05
Final 9.39503268e-13
----------------------
2. 0.50 with step 0.3999 (0.4000) converged in 2 iterations.
3. 0.90 with step 0.3990 (0.4000) converged in 2 iterations.
4. 1.28 with step 0.3939 (0.4000) converged in 3 iterations.
5. 1.54 with step 0.2966 (0.3266) converged in 3 iterations.
6. 1.68 with step 0.1724 (0.2666) converged in 3 iterations.
7. 1.75 with step 0.0832 (0.2177) converged in 2 iterations.
8. 1.79 with step 0.0540 (0.2177) converged in 2 iterations.
9. 1.82 with step 0.0368 (0.2177) converged in 2 iterations.
10. 1.85 with step 0.0263 (0.2177) converged in 2 iterations.
11. 1.86 with step 0.0196 (0.2177) converged in 3 iterations.
12. 1.87 with step 0.0124 (0.1778) converged in 2 iterations.
13. 1.88 with step 0.0104 (0.1778) converged in 2 iterations.
14. 1.89 with step 0.0091 (0.1778) converged in 2 iterations.
15. 1.90 with step 0.0081 (0.1778) converged in 2 iterations.
16. 1.91 with step 0.0074 (0.1778) converged in 2 iterations.
17. 1.91 with step 0.0070 (0.1778) converged in 2 iterations.
18. 1.92 with step 0.0066 (0.1778) converged in 2 iterations.
19. 1.93 with step 0.0063 (0.1778) converged in 2 iterations.
20. 1.93 with step 0.0062 (0.1778) converged in 2 iterations.
21. 1.94 with step 0.0060 (0.1778) converged in 2 iterations.
22. 1.94 with step 0.0059 (0.1778) converged in 2 iterations.
23. 1.95 with step 0.0058 (0.1778) converged in 2 iterations.
24. 1.96 with step 0.0058 (0.1778) converged in 2 iterations.
25. 1.96 with step 0.0057 (0.1778) converged in 2 iterations.
26. 1.97 with step 0.0057 (0.1778) converged in 2 iterations.
27. 1.97 with step 0.0057 (0.1778) converged in 2 iterations.
28. 1.98 with step 0.0056 (0.1778) converged in 2 iterations.
29. 1.98 with step 0.0056 (0.1778) converged in 2 iterations.
30. 1.99 with step 0.0056 (0.1778) converged in 2 iterations.
31. 2.00 with step 0.0056 (0.1778) converged in 2 iterations.
32. 2.00 with step 0.0056 (0.1778) converged in 2 iterations.
33. 2.01 with step 0.0056 (0.1778) converged in 2 iterations.
34. 2.01 with step 0.0056 (0.1778) converged in 2 iterations.
35. 2.02 with step 0.0056 (0.1778) converged in 2 iterations.
36. 2.02 with step 0.0056 (0.1778) converged in 2 iterations.
37. 2.03 with step 0.0057 (0.1778) converged in 2 iterations.
38. 2.04 with step 0.0057 (0.1778) converged in 2 iterations.
39. 2.04 with step 0.0057 (0.1778) converged in 2 iterations.
40. 2.05 with step 0.0058 (0.1778) converged in 2 iterations.
41. 2.05 with step 0.0058 (0.1778) converged in 2 iterations.
42. 2.06 with step 0.0059 (0.1778) converged in 2 iterations.
43. 2.06 with step 0.0059 (0.1778) converged in 2 iterations.
44. 2.07 with step 0.0061 (0.1778) converged in 2 iterations.
45. 2.08 with step 0.0062 (0.1778) converged in 2 iterations.
46. 2.08 with step 0.0064 (0.1778) converged in 2 iterations.
47. 2.09 with step 0.0066 (0.1778) converged in 2 iterations.
48. 2.10 with step 0.0070 (0.1778) converged in 2 iterations.
49. 2.11 with step 0.0074 (0.1778) converged in 2 iterations.
50. 2.11 with step 0.0081 (0.1778) converged in 2 iterations.
51. 2.12 with step 0.0090 (0.1778) converged in 2 iterations.
52. 2.13 with step 0.0103 (0.1778) converged in 2 iterations.
53. 2.15 with step 0.0122 (0.1778) converged in 2 iterations.
54. 2.16 with step 0.0149 (0.1778) converged in 2 iterations.
55. 2.19 with step 0.0186 (0.1778) converged in 2 iterations.
56. 2.21 with step 0.0239 (0.1778) converged in 2 iterations.
57. 2.25 with step 0.0316 (0.1778) converged in 2 iterations.
58. 2.30 with step 0.0432 (0.1778) converged in 2 iterations.
59. 2.38 with step 0.0611 (0.1778) converged in 2 iterations.
60. 2.48 with step 0.0879 (0.1778) converged in 3 iterations.
61. 2.59 with step 0.0991 (0.1451) converged in 3 iterations.
62. 2.69 with step 0.0968 (0.1185) converged in 2 iterations.
63. 2.80 with step 0.1055 (0.1185) converged in 2 iterations.
64. 2.91 with step 0.1108 (0.1185) converged in 2 iterations.
65. 3.03 with step 0.1138 (0.1185) converged in 2 iterations.
66. 3.14 with step 0.1156 (0.1185) converged in 2 iterations.
67. 3.26 with step 0.1166 (0.1185) converged in 2 iterations.
68. 3.38 with step 0.1172 (0.1185) converged in 2 iterations.
69. 3.50 with step 0.1176 (0.1185) converged in 2 iterations.
70. 3.61 with step 0.1179 (0.1185) converged in 1 iterations.
71. 3.78 with step 0.1670 (0.1676) converged in 2 iterations.
72. 3.95 with step 0.1672 (0.1676) converged in 1 iterations.
73. 4.18 with step 0.2366 (0.2370) converged in 2 iterations.Obtain Harmonics
uh = zeros(Complex, maximum(h)+1, length(sols));
uh[[inds0; hinds], :] = hcat([[up[zinds];up[rinds,:]+1im*up[iinds,:]] for up in sols.up]...);
Oms = [up[end] for up in sols.up];Stability Analysis
We now use an averaging formulation to obtain the stability coefficients based on just the first harmonics. This will work if the response is dominantly single harmonic.
E0, _ = HARMONICSTIFFNESS(0., -2mdl.M, 0, 1, h);
E0 = collect(E0);
stab = zeros(length(Oms));
for (i,(J,Om)) in enumerate(zip(sols.J, Oms))
# evs = eigvals(J[2:end,2:end], Om*collect(E0[2:end,2:end])); # Multiharmonic
evs = eigvals(J[2:3,2:3], Om*E0[2:3,2:3]);
stab[i] = sum(real(evs).>=0);
endThe variable stab will store the number of unstable eigenvalues that have been detected. One may interpret these as the Floquet exponents of the system.
Bifurcation Analysis
Setup QPHB for Branch Switching
Nhmax = 3;
hq = HSEL(Nhmax, [1.,1.]);
inds0q,hindsq, zindsq,rindsq,iindsq = HINDS(1, hq);
Nhcq = NHC(hq);
Nq = 32;
h01 = findall(vec(all(hq.==0, dims=2)));
hq1s = setdiff(findall(hq[:,2].==0), h01).-length(h01); # w1 alone
hq2s = setdiff(findall(hq[:,1].==0), h01).-length(h01); # w2 alone
Flq = zeros(Nhcq);
Flq[rindsq[hq2s[1]]] = pars.F;Detection
bifis = findall(stab[1:end-1].!=stab[2:end]);
bifis[2] += 1;
dxis = [-1, 1];
uhbs = [];
Omsb = [];
Wselfs = [];We loop over the two detected bifurcation points
for (bi, bifi) in enumerate(bifis)
# Eigenanalysis to obtain unstable manifold (eigenvectors)
eVals, eVecs = eigen(sols.J[bifi][2:3, 2:3], Oms[bifi]*collect(E0[2:3,2:3]));
eVecsC = eVecs[1,:]-1im.*eVecs[2,:]; # Complexify
ei = argmax(abs.(eVecsC)); # Only one should "survive" for the Hopf bifurcation
sig = imag(eVals[ei]);
Wself = Oms[bifi] + sig; # Self Excited Frequency
eVecsC = eVecsC*exp(-1im*angle(eVecsC[ei])); # Normalize phase
Pvec = normalize([real(eVecsC[ei]), -imag(eVecsC[ei])]); # Perturbation Vector
# Setup Initial Guess
Uhq0 = zeros(Nhcq);
Uhq0[[rindsq[hq2s]; iindsq[hq2s]]] = sols.up[bifi][[rinds[1:length(hq2s)];
iinds[1:length(hq2s)]]];
# Perturbation Vector
Phq0 = zeros(Nhcq);
Phq0[[rindsq[hq1s[1]]; iindsq[hq1s[1]]]] = Pvec;
# Perturbation amplitude
qamp = 1.0;
# We choose this arbitrarily for now. It is possible to use the
# method of normal forms to fix this exactly, see [next example](@ref ex_d3).
### Apply phase constraint and converge with deflation
cL = I(Nhcq+2)[:, setdiff(1:Nhcq+2, iindsq[hq1s[1]])];
uC = cL'*[Uhq0; Wself; Oms[bifi]];
# The matrix `cL` is defined so that the original solution vector `Uw`
# (which includes harmonics and both the frequency components) can be
# recovered by $Uw = cL \hat{Uw}$.
funq = NonlinearFunction((r,u,p)-> QPHBRESFUN!([u;p], mdl, Flq, hq, Nq;
R=r,cL=cL,U0=uC),
jac=(J,u,p)->QPHBRESFUN!([u;p], mdl, Flq, hq, Nq;
dRdU=J,cL=cL,U0=uC),
paramjac=(Jp,u,p)->QPHBRESFUN!([u;p], mdl, Flq, hq, Nq;
dRdw=Jp,cL=cL,U0=uC));
# The [`QPHBRESFUN!`](@ref juliajim.MDOFUTILS.QPHBRESFUN!) function
# supports deflation specification through the keyword argument U0.
u0 = cL[1:end-1,:]'*[Uhq0+qamp*Phq0; Wself];
probq = NonlinearProblem(funq, u0[1:end-1], Oms[bifi]);
solq = solve(probq, show_trace=Val(true));
# Get one point before (without deflation)
funr = NonlinearFunction((r,u,p)-> QPHBRESFUN!([u;p], mdl, Flq, hq, Nq; R=r,cL=cL),
jac=(J,u,p)->QPHBRESFUN!([u;p], mdl, Flq, hq, Nq; dRdU=J,cL=cL));
probq = NonlinearProblem(funr, u0[1:end-1], Oms[bifi-dxis[bi]]);
solprev = solve(probq, show_trace=Val(true));
### Continue away from the bifurcation point
Om0b = Oms[bifi];
Om1b = (dxis[bi]<0) ? Om0 : Om1;
dOmb = 0.2;
cparsb = (parm=:arclength, nmax=100, save_jacs=true);
solsb, _, _, _, _ = CONTINUATE(solq.u, funq, [Om0b, Om1b], dOmb; cparsb...);
# Prepend previous point & expand constraint
sup = [cL*up for up in [[solprev.u;Oms[bifi-dxis[bi]]], solsb.up...]];
# Obtain Harmonics
uhq = zeros(Complex, size(hq,1), length(sup));
uhq[[inds0q; hindsq], :] = hcat([[up[zindsq];up[rindsq,:]+1im*up[iindsq,:]]
for up in sup]...);
push!(uhbs, uhq);
push!(Omsb, [up[end] for up in sup]);
push!(Wselfs, [up[end-1] for up in sup]);
end
Algorithm: NewtonRaphson(
descent = NewtonDescent(),
autodiff = AutoForwardDiff(),
vjp_autodiff = AutoFiniteDiff(
fdtype = Val{:forward}(),
fdjtype = Val{:forward}(),
fdhtype = Val{:hcentral}(),
dir = true
),
jvp_autodiff = AutoForwardDiff(),
concrete_jac = Val{false}()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 2.74324712e-01 NaN
1 2.37721922e+00 8.79693034e-01
2 6.66794766e-01 1.89389940e-01
3 4.95056588e-01 1.26110914e-01
4 1.22594686e-02 2.07557118e-02
5 2.58972407e-05 9.93024336e-04
6 4.26687126e-10 1.77127223e-06
7 1.30732344e-15 1.32874671e-11
Final 1.30732344e-15
----------------------
Algorithm: NewtonRaphson(
descent = NewtonDescent(),
autodiff = AutoForwardDiff(),
vjp_autodiff = AutoFiniteDiff(
fdtype = Val{:forward}(),
fdjtype = Val{:forward}(),
fdhtype = Val{:hcentral}(),
dir = true
),
jvp_autodiff = AutoForwardDiff(),
concrete_jac = Val{false}()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 2.74324712e-01 0.00000000e+00
1 4.14480092e-01 1.97444792e+00
2 1.59223021e-04 1.43328242e-01
3 1.95991844e-10 3.20911518e-04
4 2.22044605e-16 1.70504345e-06
Final 2.22044605e-16
----------------------
Algorithm: NewtonRaphson(
descent = NewtonDescent(),
autodiff = AutoForwardDiff(),
vjp_autodiff = AutoFiniteDiff(
fdtype = Val{:forward}(),
fdjtype = Val{:forward}(),
fdhtype = Val{:hcentral}(),
dir = true
),
jvp_autodiff = AutoForwardDiff(),
concrete_jac = Val{false}()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 1.30732344e-15 NaN
1 1.30732344e-15 7.39353061e-16
Final 1.30732344e-15
----------------------
2. 1.54 with step -0.0032 (0.0253) converged in 2 iterations.
3. 1.54 with step -0.0037 (0.0253) converged in 2 iterations.
4. 1.53 with step -0.0042 (0.0253) converged in 2 iterations.
5. 1.53 with step -0.0047 (0.0253) converged in 2 iterations.
6. 1.52 with step -0.0052 (0.0253) converged in 2 iterations.
7. 1.52 with step -0.0058 (0.0253) converged in 2 iterations.
8. 1.51 with step -0.0064 (0.0253) converged in 2 iterations.
9. 1.50 with step -0.0070 (0.0253) converged in 2 iterations.
10. 1.49 with step -0.0077 (0.0253) converged in 2 iterations.
11. 1.48 with step -0.0084 (0.0253) converged in 1 iterations.
12. 1.47 with step -0.0129 (0.0357) converged in 2 iterations.
13. 1.46 with step -0.0145 (0.0357) converged in 2 iterations.
14. 1.44 with step -0.0161 (0.0357) converged in 2 iterations.
15. 1.42 with step -0.0178 (0.0357) converged in 2 iterations.
16. 1.40 with step -0.0196 (0.0357) converged in 2 iterations.
17. 1.38 with step -0.0214 (0.0357) converged in 2 iterations.
18. 1.35 with step -0.0231 (0.0357) converged in 2 iterations.
19. 1.33 with step -0.0248 (0.0357) converged in 2 iterations.
20. 1.30 with step -0.0263 (0.0357) converged in 2 iterations.
21. 1.27 with step -0.0278 (0.0357) converged in 2 iterations.
22. 1.24 with step -0.0290 (0.0357) converged in 2 iterations.
23. 1.21 with step -0.0301 (0.0357) converged in 2 iterations.
24. 1.18 with step -0.0310 (0.0357) converged in 2 iterations.
25. 1.15 with step -0.0318 (0.0357) converged in 2 iterations.
26. 1.12 with step -0.0325 (0.0357) converged in 2 iterations.
27. 1.08 with step -0.0330 (0.0357) converged in 2 iterations.
28. 1.05 with step -0.0335 (0.0357) converged in 1 iterations.
29. 1.00 with step -0.0479 (0.0505) converged in 2 iterations.
30. 0.95 with step -0.0485 (0.0505) converged in 2 iterations.
31. 0.90 with step -0.0489 (0.0505) converged in 2 iterations.
32. 0.85 with step -0.0493 (0.0505) converged in 1 iterations.
33. 0.78 with step -0.0701 (0.0715) converged in 2 iterations.
34. 0.71 with step -0.0705 (0.0715) converged in 2 iterations.
35. 0.64 with step -0.0707 (0.0715) converged in 2 iterations.
36. 0.57 with step -0.0709 (0.0715) converged in 2 iterations.
37. 0.50 with step -0.0711 (0.0715) converged in 1 iterations.
38. 0.40 with step -0.1007 (0.1011) converged in 2 iterations.
39. 0.30 with step -0.1009 (0.1011) converged in 1 iterations.
40. 0.17 with step -0.1262 (0.1264) converged in 2 iterations.
41. 0.04 with step -0.1263 (0.1264) converged in 2 iterations.
Algorithm: NewtonRaphson(
descent = NewtonDescent(),
autodiff = AutoForwardDiff(),
vjp_autodiff = AutoFiniteDiff(
fdtype = Val{:forward}(),
fdjtype = Val{:forward}(),
fdhtype = Val{:hcentral}(),
dir = true
),
jvp_autodiff = AutoForwardDiff(),
concrete_jac = Val{false}()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 1.19949781e-01 1.48604067e+01
1 3.57329495e+00 8.50479327e-01
2 1.50841143e+00 9.92192320e-02
3 7.38530106e-02 4.86705387e-02
4 1.44373704e-02 2.27196843e-02
5 1.02437126e-04 1.04384384e-03
6 1.89557744e-10 3.85094377e-06
7 1.79483716e-15 1.44870752e-10
Final 1.79483716e-15
----------------------
Algorithm: NewtonRaphson(
descent = NewtonDescent(),
autodiff = AutoForwardDiff(),
vjp_autodiff = AutoFiniteDiff(
fdtype = Val{:forward}(),
fdjtype = Val{:forward}(),
fdhtype = Val{:hcentral}(),
dir = true
),
jvp_autodiff = AutoForwardDiff(),
concrete_jac = Val{false}()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 2.11049557e-01 1.47121050e+01
1 9.86366454e-01 2.06015711e+00
2 5.17066673e-04 6.41802281e-02
3 4.87344383e-09 2.04408878e-04
4 1.38777878e-17 1.98810021e-07
Final 1.38777878e-17
----------------------
Algorithm: NewtonRaphson(
descent = NewtonDescent(),
autodiff = AutoForwardDiff(),
vjp_autodiff = AutoFiniteDiff(
fdtype = Val{:forward}(),
fdjtype = Val{:forward}(),
fdhtype = Val{:hcentral}(),
dir = true
),
jvp_autodiff = AutoForwardDiff(),
concrete_jac = Val{false}()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 1.79483716e-15 2.69364144e+01
1 1.79483716e-15 6.20426010e-16
Final 1.79483716e-15
----------------------
2. 2.38 with step 0.0025 (0.0225) converged in 2 iterations.
3. 2.38 with step 0.0028 (0.0225) converged in 2 iterations.
4. 2.38 with step 0.0031 (0.0225) converged in 2 iterations.
5. 2.39 with step 0.0034 (0.0225) converged in 1 iterations.
6. 2.39 with step 0.0052 (0.0318) converged in 2 iterations.
7. 2.40 with step 0.0059 (0.0318) converged in 2 iterations.
8. 2.41 with step 0.0066 (0.0318) converged in 2 iterations.
9. 2.41 with step 0.0073 (0.0318) converged in 2 iterations.
10. 2.42 with step 0.0082 (0.0318) converged in 2 iterations.
11. 2.43 with step 0.0090 (0.0318) converged in 2 iterations.
12. 2.44 with step 0.0100 (0.0318) converged in 1 iterations.
13. 2.46 with step 0.0156 (0.0450) converged in 2 iterations.
14. 2.48 with step 0.0179 (0.0450) converged in 2 iterations.
15. 2.50 with step 0.0203 (0.0450) converged in 2 iterations.
16. 2.53 with step 0.0230 (0.0450) converged in 2 iterations.
17. 2.55 with step 0.0257 (0.0450) converged in 2 iterations.
18. 2.58 with step 0.0284 (0.0450) converged in 2 iterations.
19. 2.61 with step 0.0311 (0.0450) converged in 2 iterations.
20. 2.65 with step 0.0335 (0.0450) converged in 2 iterations.
21. 2.69 with step 0.0357 (0.0450) converged in 2 iterations.
22. 2.72 with step 0.0375 (0.0450) converged in 2 iterations.
23. 2.76 with step 0.0390 (0.0450) converged in 2 iterations.
24. 2.80 with step 0.0402 (0.0450) converged in 2 iterations.
25. 2.85 with step 0.0411 (0.0450) converged in 2 iterations.
26. 2.89 with step 0.0419 (0.0450) converged in 2 iterations.
27. 2.93 with step 0.0425 (0.0450) converged in 2 iterations.
28. 2.97 with step 0.0430 (0.0450) converged in 2 iterations.
29. 3.02 with step 0.0433 (0.0450) converged in 2 iterations.
30. 3.06 with step 0.0436 (0.0450) converged in 2 iterations.
31. 3.10 with step 0.0439 (0.0450) converged in 1 iterations.
32. 3.17 with step 0.0623 (0.0637) converged in 2 iterations.
33. 3.23 with step 0.0626 (0.0637) converged in 2 iterations.
34. 3.29 with step 0.0628 (0.0637) converged in 1 iterations.
35. 3.38 with step 0.0891 (0.0901) converged in 2 iterations.
36. 3.47 with step 0.0893 (0.0901) converged in 2 iterations.
37. 3.56 with step 0.0895 (0.0901) converged in 1 iterations.
38. 3.67 with step 0.1120 (0.1126) converged in 2 iterations.
39. 3.79 with step 0.1122 (0.1126) converged in 1 iterations.
40. 3.90 with step 0.1123 (0.1126) converged in 1 iterations.
41. 4.01 with step 0.1124 (0.1126) converged in 1 iterations.Plotting
set_theme!(theme_latexfonts())
fsz = 24;
fig = Figure(fontsize=fsz);
ax = Axis(fig[1, 1], xlabel=L"Excitation Frequency $\Omega$", ylabel="Response");
scatterlines!(ax, Oms./(stab.==0), [norm(u) for u in eachcol(uh)], label="Stable")
scatterlines!(ax, Oms./(stab.==2), [norm(u) for u in eachcol(uh)], label="Unstable")
for (i, (om,uhb)) in enumerate(zip(Omsb, uhbs))
scatterlines!(ax, om, [norm(u) for u in eachcol(uhb)],
label="Branch $(i)")
end
axislegend(ax, nbanks=3, position=:ct)
xlims!(Om0, Om1)
ylims!(0, 3.75)
fig
This page was generated using Literate.jl.