# How to get the frequency response with the unstable branches for a nonlinear driven system?

I'm working on a driven system and want to get the amplitude-frequency response curve with an unstable branch just like the following one

where the dashed lines correspond to the unstable branches and the solid lines correspond to the stable branches.
My codes are

Clear["Global*"]
α = 110.; β = 55.; δ = 1.; μ1 = 18.; μ2 = 42.; μ = μ2/μ1;
ηb = 10.;
ωlist = Range[2.5, 6., 0.1];
deltap = .22;
p0 = 0.2;
inipoint = 1.19;
tlength = 1000.;
w[λ_, ξ_] := (-((μ1*α)/2) Log[
1 - (λ^(-4) + 2*λ^2 -
3)/α] - (μ2*β)/2 Log[
1 - (λ^-4*ξ^4 + 2 λ^2*ξ^-2 -
3)/β])/μ1
dw[λ_, ξ_] = D[w[λ, ξ], λ];

ξin[λ_, ξ_,
x_] = (1 + (λ^3 - 1) (x^3 - 1)^-1 (ξ^3 - 1))^(1/3);
f[λ_, ξ_, x_] =
dw[λ, ξin[λ, ξ, x]]/(1 - λ^3);
sup[x_] := ((δ + x^3)/(1 + δ))^(1/3)

Get["NumericalDifferentialEquationAnalysis"];
np = 11; points = weights = Table[Null, {np}];
intf[x0_, ξ0_] :=
Block[{y = x0, ξ1 = ξ0},
Do[points[[i]] =
GaussianQuadratureWeights[np, y, sup[y]][[i, 1]], {i, 1, np}];
Do[weights[[i]] =
GaussianQuadratureWeights[np, y, sup[y]][[i, 2]], {i, 1, np}];
int = Sum[(f[λ, ξ1, y] /. λ -> points[[i]])*
weights[[i]], {i, 1, np}]; int]

eqns = {x'[t] == y[t],
y'[t] == -(1/
2 x'[t]^2 (3 - δ/
x[t]^3 (1 + δ/x[t]^3)^(-4/3) -
3 (1 + δ/x[t]^3)^(-1/3)) + intf[x[t], z[t]] -
deltap - p0*Sin[ω*t])/
x[t]/(1 - (1 + δ/x[t]^3)^(-1/3)),
z'[t] ==
z[t]*(μ (x[t]^2*z[t]^-2 -
x[t]^-4*z[t]^4))/(3 ηb*(1 - (x[t]^-4*z[t]^4 +
2 x[t]^2*z[t]^-2 - 3)/β))};

AbsoluteTiming@ListPlot[Level[#, {2}] &@(bifdata = ParallelTable[
reapx =
Reap[NDSolve[{eqns, {x[0] == inipoint, y[0] == 0,
z[0] == inipoint},
WhenEvent[y[t] == 0 && t > 0.95 tlength,
Sow[x[t]]]}, {}, {t, 0, tlength},
Method -> {"EquationSimplification" -> "Residual"}]][[2, 1]];

Transpose[{Table[ω, Length[reapx]],
reapx}], {ω, ωlist}]), PlotRange -> All]
data = Level[bifdata, {2}];
ave = {#[[1, 1]], Max@#[[All, 2]] - Min@#[[All, 2]]} & /@
GatherBy[data,
First];(*The amplitude is calculated by the difference of the \
maximum and minimum*)
ListPlot[ave, PlotRange -> All, PlotStyle -> Black,
PlotMarkers -> "\[Star]"]


where x is the displacement variable.
The amplitude-response curve I obtained from codes above is

where axis x is the external excitation frequency ω.
As depicted in my result, only the stable branch is given, how can I get the result with the unstable branch part?
Any suggestions would be much appreciated!

• The spectral properties of the model depend on the np. For example, when np=41, excitation occurs already at t> 60. – Alex Trounev Jul 18 at 7:09
• @AlexTrounev Thanks for your comment! Did you mean the integration method adopted here is responsible for that? Is there any way to get the results I want? – keanhy14 Jul 18 at 7:27
• To calculate the ambiguous dependence of amplitude on frequency, it is necessary to build a very accurate numerical model. – Alex Trounev Jul 18 at 7:33
• @AlexTrounev Yeah, my equations are quite complex. The strong nonlinearities make it hard to deal with it analytically, just like the Duffing equation. I learn from papers the unstable limit cycle is used to calculate the unstable branch. – keanhy14 Jul 18 at 7:48
• I am studying this in detail and would like to know how you got your differential equation from the standard Duffing equation. Can you give a reference or if it is straightforward put it in your question? Thanks – Hugh Jul 29 at 8:18

Here's a partial answer, because I suspect the OP's code doesn't correspond perfectly with their figure, resulting in a more complicated situation.

The main idea is to make a stroboscopic map that advances the system by one period. This can be treated as a discrete-time dynamical system. Equilibria (corresponding to limit cycles of the original system) can be found by FindRoot, including unstable cycles. Finally we can implement a crude continuation technique to track an equilibrium across the parameter, by using the previous answer as the initial guess for FindRoot.

First, a helper to find extrema:

FindExtrema[if_InterpolatingFunction, {tmin_?NumericQ, tmax_?NumericQ}] :=
Reap[NDSolve[{y'[t] == Evaluate[D[if[t], t]],
WhenEvent[y'[t] == 0, Sow[{t, y[t]}]],
y[tmin] == if[tmin]}, y[t], {t, tmin, tmax}]][[2]] /. {x_List} :> x;

FindExtrema[if_InterpolatingFunction] := FindExtrema[if, if["Domain"][[1]]];


The model is set up as in the original post (not copied here), except I changed the second equation to use y[t] instead of x'[t] to avoid the need for Method -> {"EquationSimplification" -> "Residual"}]:

eqns = {
x'[t] == y[t],
y'[t] == -(1/2 y[t]^2 (3 - δ/x[t]^3 (1 + δ/x[t]^3)^(-4/3) -
3 (1 + δ/x[t]^3)^(-1/3)) + intf[x[t], z[t]] - deltap - p0*Sin[ω*t])/x[t]/(1 - (1 + δ/x[t]^3)^(-1/3)),
z'[t] == z[t]*(μ (x[t]^2*z[t]^-2 - x[t]^-4*z[t]^4))/(3 ηb*(1 - (x[t]^-4*z[t]^4 + 2 x[t]^2*z[t]^-2 - 3)/β))};


Now, define the stroboscopic map F:

F[{x0_?NumericQ, y0_?NumericQ, z0_?NumericQ}] := (
sol = NDSolve[{eqns, {x[0] == x0, y[0] == y0, z[0] == z0}},
{x, y, z}, {t, 0, 2 π/ω}][[1]];
{x[2 π/ω], y[2 π/ω], z[2 π/ω]} /. sol
);


Next a function to find fixed points of the map:

findEq[{x0i_?NumericQ, y0i_?NumericQ, z0i_?NumericQ}, opts___] :=
{x0, y0, z0} /.
FindRoot[F[{x0, y0, z0}] == {x0, y0, z0}, {{x0, x0i}, {y0, y0i}, {z0, z0i}}, opts]


In action:

ω = 2.5;
eq = findEq[{1.2, 0.1, 1.2}]
(* {1.23676, 0.169015, 1.20327} *)
Plot[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, 2 π/ω}]


Since FindRoot may throw some FindRoot::lstol warnings, we can check whether the equilibrium we found is good:

F[eq] - eq
(* {-8.50313*10^-9, -4.05657*10^-8, 4.94302*10^-10} *)


Close enough!

We can also check the stability of a cycle with a finite-difference approximation to the Jacobian matrix of F:

j := {(F[eq + {ϵ, 0, 0}] - F[eq - {ϵ, 0, 0}])/(2 ϵ),
(F[eq + {0, ϵ, 0}] - F[eq - {0, ϵ, 0}])/(2 ϵ),
(F[eq + {0, 0, ϵ}] - F[eq - {0, 0, ϵ}])/(2 ϵ)};

ϵ = 10^-5;
Chop[Eigenvalues[j]]
(* {0.852387, 0.0131683 + 0.598595 I, 0.0131683 - 0.598595 I} *)


All real parts less than one in magnitude indicates stability. ϵ needs to be small enough to make this a good linearization of F' but not too small, to avoid numerical roundoff errors.

Now we'll make four tracks of the equilibria with different starting values of ω:

ω = 2.5;
ics = findEq[{1.2367572037134027, 0.169014899645155, 1.2032719485392172}];
Clear[ω];

res1 = Table[
ics = eq1[ω] = findEq[ics];
{ω, Abs[Differences[FindExtrema[x /. sol][[All, 2]]]][[1]]}
, {ω, 2.5, 3.27, 0.01}];

ω = 3.5;
ics = findEq[{1.7, -3, 1.9}];
Clear[ω];

res2 = Table[
ics = eq2[ω] = findEq[ics];
{ω, Abs[Differences[FindExtrema[x /. sol][[All, 2]]]][[1]]}
, {ω, 3.5, 2.7, -0.01}];

ω = 3.0;
ics = findEq[{0.81, -1.31, 1.51}, DampingFactor -> 0.1];
Clear[ω];

res3 = Table[
ics = eq3[ω] = findEq[ics];
{ω, Abs[Differences[FindExtrema[x /. sol][[All, 2]]]][[1]]}
, {ω, 3.0, 3.27, 0.01}];

ω = 3.0;
ics = findEq[{0.81, -1.31, 1.51}, DampingFactor -> 0.1];
Clear[ω];

res4 = Table[
ics = eq4[ω] = findEq[ics];
{ω, Abs[Differences[FindExtrema[x /. sol][[All, 2]]]][[1]]}
, {ω, 3.0, 2.5, -0.01}];


res3 and res4 are an unstable cycle. Getting the initial point for these tracks is the hardest part and involves some trial and error.

Plotting the results:

ListPlot[{res1, res2, res3, res4}]


Note the discrepancy between these results and OP's figure. The unstable branch doesn't connect up with the upper stable branch, so I suspect there is a lot more going on in this system. But without a reference to the source of the figure, it's hard to know how to proceed.

By the way, here's an alternative equilibrium tracking code that uses linear extrapolation of the last two answers to get an even better initial guess for FindRoot:

ω = 2.5;
ics′ = ics = findEq[{1.2367572037134027, 0.169014899645155, 1.2032719485392172}];
Clear[ω];

res1′ = Table[
eq1[ω] = findEq[2 ics - ics′];
ics′ = ics; ics = eq1[ω];
{ω, Abs[Differences[FindExtrema[x /. sol][[All, 2]]]][[1]]}
, {ω, 2.5, 3.27, 0.01}];

• Thank you very much for your suggestion, and it really helps a lot! My problem is similar to that of this paper. In section 4.3, the stability of a periodic solution is evaluated using Floquet theory. In the sixth paragraph of section 5, it is pointed out that "While the sequential continuation method fails at the bifurcation points of the diagram, it is enough to start the numerical procedure with a stable (unstable) limit cycle to calculate the whole stable (unstable) branch." – keanhy14 Jul 22 at 11:30
• The warning for codes you presented is that FindRoot: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. My MMA is V11.3. Does it affect the result? – keanhy14 Jul 22 at 11:34
• The stroboscopic map approach should be equivalent to the Floquet theory approach. I'll edit my answer to show you how to check stability of an equilibrium and also how to implement linear extrapolation of the last two answers for the initial guess, which should be a little better than just using the last answer. Finally, re: the FindRoot warnings, you might doublecheck the quality of your equilibria. – Chris K Jul 22 at 13:16
• BTW, don't forget that your model might also have other bifurcations, such as period-doubling or Hopf bifurcations, that won't be obvious like the turning points. Better check that with the Jacobian matrix. Could get complicated! – Chris K Jul 22 at 13:29
• I find a Demo here, and It gives the unstable solution of the limit cycle. – keanhy14 Jul 23 at 11:57

If I understand your question rigth, you are looking for steadystate solutions with frequency \[Omega]?
If so this approach might help you to find a numerical model:

First solve the system numerically (sufficient to use last period ) depending on \[Omega]

sol = ParametricNDSolveValue[{eqns, {x[0] == inipoint, y[0] == 0,z[0] == inipoint}}, {x, y, z}, {t, tlength - (2 Pi)/\[Omega], tlength}, {\[Omega]},Method -> {"EquationSimplification" -> "Residual"}]


Now plot solution for different frequencies

Show[
Table[Function[{\[Omega]},
ParametricPlot[{sol[\[Omega]][[2]][t], sol[\[Omega]][[1]][t]}, {t,tlength - (2 Pi)/\[Omega], tlength},PlotStyle ->Hue[\[Omega]]]][\[Omega]]
, {\[Omega], \[Omega]list}],
PlotRange -> {0,3}, Evaluated -> True, AxesOrigin -> {0, 0}, AxesLabel -> {"y[t]","x[t]"}]


The single curves contain the amplitudes you're looking for...

• Yes, this is a stable period solution for different excitation frequency, as shown in my result (fig. 2). But the unstable solution is what I exactly want (the dashed line in fig. 1). One can refer to the Duffing equation in Wikipedia. – keanhy14 Jul 18 at 10:32
• But you expect a stationary solution with frequency \[Omega] for the unstable case? – Ulrich Neumann Jul 18 at 10:48
• Yes, I need the unstable case. The more commonly seen case is the analytic form solution. With CounterPlot, the stable and unstable cases are both plotted, but my problem can only be solved by numerical method. – keanhy14 Jul 18 at 11:23