GitHub

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.NLDYN

System 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.1;
cpars = (parm=:arclength, nmax=400, 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.20 with step 0.0995 (0.9973) converged in 1 iterations. 0.174321 Δα.
3. 0.32 with step 0.1202 (0.8514) converged in 1 iterations. 0.138045 Δα.
4. 0.48 with step 0.1623 (0.7642) converged in 1 iterations. 0.138009 Δα.
5. 0.70 with step 0.2187 (0.6859) converged in 2 iterations. 1.110848 Δα.
6. 0.89 with step 0.1845 (0.4140) converged in 3 iterations. 1.859535 Δα.
7. 1.02 with step 0.1349 (0.2230) converged in 2 iterations. 0.329273 Δα.
8. 1.15 with step 0.1387 (0.1995) converged in 2 iterations. 0.402833 Δα.
9. 1.28 with step 0.1423 (0.1995) converged in 2 iterations. 0.431273 Δα.
10. 1.40 with step 0.1317 (0.1995) converged in 2 iterations. 0.401732 Δα.
11. 1.50 with step 0.1139 (0.1995) converged in 2 iterations. 0.402192 Δα.
12. 1.58 with step 0.0955 (0.1995) converged in 2 iterations. 1.569040 Δα.
13. 1.65 with step 0.0793 (0.1995) converged in 2 iterations. 2.754150 Δα.
14. 1.70 with step 0.0652 (0.1995) converged in 2 iterations. 1.381343 Δα.
15. 1.74 with step 0.0514 (0.1995) converged in 3 iterations. 1.009740 Δα.
16. 1.77 with step 0.0369 (0.1995) converged in 2 iterations. 0.743120 Δα.
17. 1.79 with step 0.0250 (0.1995) converged in 2 iterations. 0.542324 Δα.
18. 1.81 with step 0.0174 (0.1995) converged in 2 iterations. 0.408723 Δα.
19. 1.82 with step 0.0128 (0.1995) converged in 2 iterations. 0.321505 Δα.
20. 1.83 with step 0.0099 (0.1995) converged in 2 iterations. 0.262470 Δα.
21. 1.83 with step 0.0080 (0.1995) converged in 2 iterations. 0.220706 Δα.
22. 1.84 with step 0.0068 (0.1995) converged in 2 iterations. 0.192668 Δα.
23. 1.84 with step 0.0063 (0.1995) converged in 2 iterations. 0.182747 Δα.
24. 1.85 with step 0.0061 (0.1995) converged in 2 iterations. 0.179944 Δα.
25. 1.86 with step 0.0060 (0.1995) converged in 2 iterations. 0.180988 Δα.
26. 1.86 with step 0.0060 (0.1995) converged in 2 iterations. 0.185655 Δα.
27. 1.87 with step 0.0060 (0.1995) converged in 2 iterations. 0.195853 Δα.
28. 1.87 with step 0.0061 (0.1995) converged in 2 iterations. 0.216223 Δα.
29. 1.88 with step 0.0062 (0.1995) converged in 2 iterations. 0.256397 Δα.
30. 1.88 with step 0.0064 (0.1995) converged in 2 iterations. 0.337636 Δα.
31. 1.89 with step 0.0066 (0.1995) converged in 2 iterations. 0.519420 Δα.
32. 1.90 with step 0.0069 (0.1995) converged in 2 iterations. 1.079731 Δα.
33. 1.90 with step 0.0074 (0.1995) converged in 2 iterations. 0.354930 Δα.
34. 1.91 with step 0.0081 (0.1995) converged in 2 iterations. 1.291400 Δα.
35. 1.92 with step 0.0089 (0.1995) converged in 2 iterations. 0.699971 Δα.
36. 1.93 with step 0.0099 (0.1995) converged in 2 iterations. 0.501397 Δα.
37. 1.94 with step 0.0109 (0.1995) converged in 2 iterations. 0.406848 Δα.
38. 1.95 with step 0.0113 (0.1995) converged in 3 iterations. 0.362077 Δα.
39. 1.96 with step 0.0104 (0.1995) converged in 2 iterations. 0.344552 Δα.
40. 1.97 with step 0.0087 (0.1995) converged in 2 iterations. 0.341964 Δα.
41. 1.98 with step 0.0068 (0.1995) converged in 2 iterations. 0.346703 Δα.
42. 1.98 with step 0.0053 (0.1995) converged in 2 iterations. 0.355206 Δα.
43. 1.99 with step 0.0040 (0.1995) converged in 2 iterations. 0.366276 Δα.
44. 1.99 with step 0.0030 (0.1995) converged in 2 iterations. 0.380148 Δα.
45. 1.99 with step 0.0022 (0.1995) converged in 1 iterations. 0.398028 Δα.
46. 1.99 with step 0.0016 (0.1995) converged in 1 iterations. 0.423551 Δα.
47. 2.00 with step 0.0015 (0.1995) converged in 1 iterations. 0.597198 Δα.
48. 2.00 with step 0.0015 (0.1995) converged in 1 iterations. 0.007166 Δα.
49. 2.00 with step 0.0022 (0.2992) converged in 1 iterations. 0.002937 Δα.
50. 2.00 with step 0.0034 (0.4488) converged in 2 iterations. 0.593666 Δα.
51. 2.01 with step 0.0020 (0.2709) converged in 1 iterations. 0.303618 Δα.
52. 2.01 with step 0.0016 (0.2080) converged in 1 iterations. 0.230889 Δα.
53. 2.01 with step 0.0015 (0.1995) converged in 1 iterations. 0.196755 Δα.
54. 2.01 with step 0.0017 (0.1995) converged in 1 iterations. 0.186357 Δα.
55. 2.01 with step 0.0020 (0.1995) converged in 1 iterations. 0.183955 Δα.
56. 2.01 with step 0.0023 (0.1995) converged in 1 iterations. 0.184878 Δα.
57. 2.02 with step 0.0027 (0.1995) converged in 2 iterations. 0.187595 Δα.
58. 2.02 with step 0.0032 (0.1995) converged in 2 iterations. 0.191394 Δα.
59. 2.02 with step 0.0037 (0.1995) converged in 2 iterations. 0.196695 Δα.
60. 2.03 with step 0.0045 (0.1995) converged in 2 iterations. 0.204209 Δα.
61. 2.03 with step 0.0053 (0.1995) converged in 2 iterations. 0.215592 Δα.
62. 2.04 with step 0.0065 (0.1995) converged in 2 iterations. 0.234349 Δα.
63. 2.05 with step 0.0079 (0.1995) converged in 2 iterations. 0.268407 Δα.
64. 2.06 with step 0.0096 (0.1995) converged in 2 iterations. 0.336514 Δα.
65. 2.07 with step 0.0115 (0.1995) converged in 3 iterations. 0.480075 Δα.
66. 2.08 with step 0.0125 (0.1995) converged in 3 iterations. 0.793388 Δα.
67. 2.10 with step 0.0119 (0.1995) converged in 3 iterations. 1.938434 Δα.
68. 2.11 with step 0.0107 (0.1995) converged in 2 iterations. 2.508885 Δα.
69. 2.12 with step 0.0097 (0.1995) converged in 2 iterations. 0.727269 Δα.
70. 2.13 with step 0.0091 (0.1995) converged in 2 iterations. 0.438129 Δα.
71. 2.14 with step 0.0087 (0.1995) converged in 2 iterations. 0.354290 Δα.
72. 2.15 with step 0.0085 (0.1995) converged in 2 iterations. 0.327453 Δα.
73. 2.16 with step 0.0084 (0.1995) converged in 2 iterations. 0.317052 Δα.
74. 2.17 with step 0.0084 (0.1995) converged in 2 iterations. 0.310239 Δα.
75. 2.18 with step 0.0084 (0.1995) converged in 2 iterations. 0.302837 Δα.
76. 2.19 with step 0.0104 (0.1995) converged in 2 iterations. 0.377712 Δα.
77. 2.21 with step 0.0142 (0.1995) converged in 2 iterations. 0.517837 Δα.
78. 2.24 with step 0.0210 (0.1995) converged in 2 iterations. 0.779882 Δα.
79. 2.29 with step 0.0350 (0.1995) converged in 3 iterations. 1.263138 Δα.
80. 2.36 with step 0.0580 (0.1995) converged in 3 iterations. 2.114778 Δα.
81. 2.46 with step 0.0805 (0.1995) converged in 3 iterations. 4.157232 Δα.
82. 2.59 with step 0.1034 (0.1995) converged in 2 iterations. 0.532017 Δα.
83. 2.76 with step 0.1304 (0.1995) converged in 2 iterations. 0.556374 Δα.
84. 2.96 with step 0.1615 (0.1995) converged in 2 iterations. 0.606976 Δα.
85. 3.20 with step 0.1960 (0.1995) converged in 3 iterations. 0.616528 Δα.
86. 3.47 with step 0.2334 (0.1995) converged in 2 iterations. 0.605422 Δα.
87. 3.80 with step 0.2734 (0.1995) converged in 2 iterations. 0.585819 Δα.
88. 4.17 with step 0.3163 (0.1995) converged in 2 iterations. 0.564063 Δα.

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 = sols.p;

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

The 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 = 1.8;
    cparsb = (parm=:arclength, nmax=500, save_jacs=true, angopt=deg2rad(30));

    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       	3.44474909e-01      	NaN
1       	5.66308616e-01      	6.97325573e-01
2       	9.55852563e-02      	1.00037857e-01
3       	5.05119252e-03      	1.20639218e-02
4       	1.32882116e-06      	3.28883375e-04
5       	1.59963802e-12      	2.64082886e-07
6       	7.84118341e-16      	9.89129377e-14
Final   	7.84118341e-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       	3.44474909e-01      	0.00000000e+00
1       	5.78420741e-01      	2.47646059e+00
2       	1.56069464e-04      	3.25297702e-01
3       	4.21863789e-11      	7.06513294e-05
4       	5.20417043e-18      	6.93685014e-07
Final   	5.20417043e-18
----------------------

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       	7.84118341e-16      	NaN
1       	7.84118341e-16      	8.02333993e-17
Final   	7.84118341e-16
----------------------
2. 1.48 with step -0.0122 (0.0990) converged in 2 iterations. 0.340666 Δα.
3. 1.46 with step -0.0180 (0.1109) converged in 2 iterations. 0.388931 Δα.
4. 1.44 with step -0.0250 (0.1197) converged in 2 iterations. 0.425901 Δα.
5. 1.40 with step -0.0351 (0.1263) converged in 2 iterations. 0.480481 Δα.
6. 1.34 with step -0.0488 (0.1290) converged in 2 iterations. 0.548092 Δα.
7. 1.27 with step -0.0651 (0.1276) converged in 2 iterations. 0.613157 Δα.
8. 1.18 with step -0.0806 (0.1227) converged in 2 iterations. 0.651478 Δα.
9. 1.08 with step -0.0909 (0.1164) converged in 2 iterations. 0.650226 Δα.
10. 0.98 with step -0.0945 (0.1104) converged in 2 iterations. 0.617862 Δα.
11. 0.88 with step -0.0930 (0.1060) converged in 2 iterations. 0.570627 Δα.
12. 0.79 with step -0.0890 (0.1038) converged in 2 iterations. 0.520560 Δα.
13. 0.71 with step -0.0843 (0.1039) converged in 2 iterations. 0.473782 Δα.
14. 0.63 with step -0.0792 (0.1066) converged in 3 iterations. 0.427645 Δα.
15. 0.56 with step -0.0694 (0.1025) converged in 2 iterations. 0.399960 Δα.
16. 0.49 with step -0.0674 (0.1099) converged in 1 iterations. 0.371676 Δα.
17. 0.41 with step -0.0780 (0.1429) converged in 1 iterations. 0.349569 Δα.
18. 0.32 with step -0.0891 (0.1885) converged in 1 iterations. 0.335404 Δα.
19. 0.22 with step -0.0977 (0.2511) converged in 1 iterations. 0.332631 Δα.
20. 0.12 with step -0.0986 (0.3350) converged in 1 iterations. 0.350404 Δα.
21. 0.04 with step -0.0839 (0.4417) converged in 3 iterations. 1.084036 Δα.

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.76831535e-01      	1.00000000e+00
1       	1.91403875e-01      	4.84441428e-01
2       	2.95865087e-02      	6.23477517e-02
3       	8.35111610e-05      	5.26915228e-03
4       	2.29527117e-08      	7.92745727e-05
5       	1.13710841e-15      	6.86478020e-09
Final   	1.13710841e-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.39150554e-01      	8.58504819e-01
1       	1.54337530e+00      	2.54669662e+00
2       	3.45765657e-04      	5.37612939e-02
3       	1.18971278e-09      	1.15163992e-04
4       	1.11022302e-16      	1.95671946e-05
Final   	1.11022302e-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.13710841e-15      	6.65389791e+01
1       	1.13710841e-15      	3.41975271e-16
Final   	1.13710841e-15
----------------------
2. 2.48 with step 0.0145 (0.0655) converged in 2 iterations. 0.468692 Δα.
3. 2.50 with step 0.0186 (0.0674) converged in 2 iterations. 0.534294 Δα.
4. 2.52 with step 0.0213 (0.0670) converged in 2 iterations. 0.560637 Δα.
5. 2.55 with step 0.0239 (0.0659) converged in 2 iterations. 0.584805 Δα.
6. 2.58 with step 0.0261 (0.0641) converged in 2 iterations. 0.603962 Δα.
7. 2.61 with step 0.0279 (0.0619) converged in 2 iterations. 0.617070 Δα.
8. 2.64 with step 0.0293 (0.0595) converged in 2 iterations. 0.624492 Δα.
9. 2.67 with step 0.0302 (0.0570) converged in 2 iterations. 0.627346 Δα.
10. 2.70 with step 0.0308 (0.0545) converged in 1 iterations. 0.626823 Δα.
11. 2.74 with step 0.0379 (0.0634) converged in 1 iterations. 0.625691 Δα.
12. 2.79 with step 0.0471 (0.0739) converged in 2 iterations. 0.626215 Δα.
13. 2.84 with step 0.0485 (0.0707) converged in 2 iterations. 0.624407 Δα.
14. 2.90 with step 0.0495 (0.0677) converged in 2 iterations. 0.619712 Δα.
15. 2.95 with step 0.0502 (0.0650) converged in 2 iterations. 0.613491 Δα.
16. 3.00 with step 0.0536 (0.0625) converged in 2 iterations. 0.641087 Δα.
17. 3.07 with step 0.0575 (0.0595) converged in 2 iterations. 0.680443 Δα.
18. 3.13 with step 0.0608 (0.0559) converged in 2 iterations. 0.721718 Δα.
19. 3.20 with step 0.0633 (0.0517) converged in 2 iterations. 0.763912 Δα.
20. 3.26 with step 0.0646 (0.0473) converged in 2 iterations. 0.805865 Δα.
21. 3.33 with step 0.0648 (0.0427) converged in 2 iterations. 0.846408 Δα.
22. 3.40 with step 0.0637 (0.0381) converged in 2 iterations. 0.884515 Δα.
23. 3.46 with step 0.0614 (0.0337) converged in 2 iterations. 0.919433 Δα.
24. 3.52 with step 0.0582 (0.0296) converged in 2 iterations. 0.950733 Δα.
25. 3.58 with step 0.0542 (0.0258) converged in 1 iterations. 0.977417 Δα.
26. 3.64 with step 0.0616 (0.0277) converged in 1 iterations. 1.001468 Δα.
27. 3.71 with step 0.0700 (0.0296) converged in 1 iterations. 1.024554 Δα.
28. 3.79 with step 0.0794 (0.0315) converged in 2 iterations. 1.048345 Δα.
29. 3.87 with step 0.0723 (0.0269) converged in 1 iterations. 1.069056 Δα.
30. 3.95 with step 0.0808 (0.0284) converged in 1 iterations. 1.087201 Δα.
31. 4.04 with step 0.0900 (0.0299) converged in 1 iterations. 1.101983 Δα.

Plotting

set_theme!(theme_latexfonts())
fsz = 24;
fig = Figure(fontsize=fsz);

ax = Axis(fig[1, 1], xlabel=L"Excitation Frequency $\Omega$", ylabel="RMS Response",
    yscale=log10);
scatterlines!(ax, Oms./(stab.==0), [norm(u)/sqrt(2) for u in eachcol(uh)], label="Stable")
scatterlines!(ax, Oms./(stab.==2), [norm(u)/sqrt(2) for u in eachcol(uh)], label="Unstable")

for (i, (om,uhb)) in enumerate(zip(Omsb, uhbs))
    scatterlines!(ax, om, [norm(u)/sqrt(2) for u in eachcol(uhb)],
        label="B-$(i)")
    # scatterlines!(ax, om[2:end], abs.(uhb[12,2:end]),
    #     label="Branch $(i)")
end

Legend(fig[2, 1], ax, nbanks=4, position=:cb, tellheight=true, tellwidth=false)
xlims!(Om0, Om1)

    fig
Example block output

Branches B-1 and B-2 are the quasi-periodic branches that we have computed.

This page was generated using Literate.jl.