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.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      	3.31662479e+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 = 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       	4.85859077e-01      	1.13565321e+01
1       	4.64404850e-03      	2.96902753e-02
2       	9.24488481e-06      	3.59088833e-04
3       	1.74348449e-11      	5.27975403e-07
4       	3.58967432e-15      	1.40062807e-12
Final   	3.58967432e-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.24690326691878337
 0.24688194582741535

This page was generated using Literate.jl.