3: Branch Switching Using Normal Forms

Example D3: Branch Switching Using Normal Forms

This example is a direct continuation of Example D2. In that example we did a branch switching analysis by arbitrarily choosing the amplitude of perturbation. Here, we use the NORMALFORMFIT function to fit the system to a cubic nonlinear "normal form" (through a similarity transformation eliminating certain terms).

The normal form can easily be interpreted to provide an estimate for the bifurcated branch amplitudes. It can be seen that the normal form approach provides a near exact estimate of the bifurcated branch amplitude in the end of this example.

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

Om0 = 0.1;
Om1 = 4.0;
dOm = 0.1;
cpars = (parm=:arclength, nmax=200, 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      	1.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

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:3,2:3], Om*E0[2:3,2:3]);
    stab[i] = sum(real(evs).>=0);
end

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];

bi = 2;
bifi = bifis[bi];
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
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)]]];

Fit the HB Residue to a Local Normal form

Fl1 = Fl[[rinds[1], iinds[1]]];  # Single harmonic excitation
uh1 = sols.up[bifi][[rinds[1], iinds[1]]];  # Bifurcated Solution
vh1 = [normalize(real(eVecs[:,ei])) normalize(imag(eVecs[:,ei]))];  # Projector

uh0 = sols.up[bifi-dxis[bi]][[rinds[1], iinds[1]]];  # Previous Solution
d_amp = norm(vh1'*(uh1-uh0));

E0f, _ = HARMONICSTIFFNESS(0, -2mdl.M, 0, Oms[bifi], 1)
xyfun = (u) -> vh1\ (E0f\HBRESFUN!([uh1+vh1*u; Oms[bifi]], mdl, Fl1, 1, N));

A, B, C, H = NORMALFORMFIT(xyfun, d_amp/10, 100);
λ = sum(diag(A))/2;
α = sum(C[[1,4,5,8]])/4;

pVal, pVec = eigen(A);
vi = argmax(abs.(pVec[1,:]+1im*pVec[2,:]));
Pvec = vh1*pVec[:,vi];
Pvec = normalize([real(pVec[1]+1im*pVec[2]), imag(pVec[1]+1im*pVec[2])]);

Perturbation Vector

Phq0 = zeros(Nhcq);
Phq0[[rindsq[hq1s[1]]; iindsq[hq1s[1]]]] = Pvec;

Perturbation Amplitude

qamp = √(-λ/α);

Apply phase constraint and converge with deflation

cL = I(Nhcq+2)[:, setdiff(1:Nhcq+2, iindsq[hq1s[1]])];
uC = cL'*[Uhq0; Wself; Oms[bifi]];

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

u0 = cL'*[Uhq0+qamp*Phq0; Wself; Oms[bifi]];
probq = NonlinearProblem(funq, u0[1:end-1], u0[end]);
solq = solve(probq, 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.05849397e-01      	NaN
1       	3.63871578e-03      	4.45247145e-02
2       	1.14523870e-05      	1.00234607e-03
3       	3.29219135e-11      	2.59894716e-06
4       	1.14411740e-15      	1.17463572e-11
Final   	1.14411740e-15
----------------------

We can now verify that the converged solution does indeed have an amplitude that's almost exactly equal to the normal form prediction.

[cL[2,:]'*[solq.u; Oms[bifi]], qamp]

2-element Vector{Float64}:
 0.5782945663043191
 0.5781654860030309

This page was generated using Literate.jl.