# Bifurcation with a system of equations

I'm trying to plot a flip bifurcation diagram for a dynamical system of equations as follows:

x'[t] == v[t]
v'[t] == x[t] - A x[t]^3 + R*Cos[ω*t] - B v[t]


but am having trouble writing the code to get a bifurcation diagram to plot. I have numerical values for A,B, & Omega, and am trying to plot the diagram for R.

Other similar questions and references online don't seem to be helpful (either not nonlinear or not a system of equations).

• some code here mathematica.stackexchange.com/questions/13277/… – Nasser Nov 21 '13 at 4:44
• This seems to be the Rayleigh equation. You might be interested to look up its bifurcations. They should be studied pretty well. By the way, you might want to remove one parameter by rescaling. – Alexei Boulbitch Nov 21 '13 at 10:02

You will have to solve this equation using NDSolve and the several plot functions in Mathematica can take care of the visualization. Here is a sample code:

Let me put some number for the parameters in the equation:

A = 10;
B = 5;
R = 2;
\[Omega] = 2;


Using NDSolve to find the numeric solution, which is a interpolation object in Mathematica:

sol =
NDSolve[{
x'[t] == v[t], v'[t] == x[t] - A x[t]^3 + R*Cos[\[Omega]*t] - B v[t],
x[0] == 0, v[0] == 0
},
{x[t], v[t]}, {t, 0, 10}]


To extract the solution, use ReplaceAll (/.) and Rule (->) functions. Here you have

{x[t], v[t]} /. sol[[1]]


Now insert the above to ParametricPlot function:

ParametricPlot[Evaluate[{x[t], v[t]} /. sol[[1]]], {t, 0, 10}]


For this particular set of parameters, I have:

More interesting plots can be obtained for the values of {A -> 10, B -> 5, ω -> 2} used in the earlier answer but with larger R, for instance with R == 10000,

eq = {x'[t] == v[t], v'[t] == x[t] - A x[t]^3 + R*Cos[ω*t] - B v[t], x[0] == 0, v[0] == 0};
s = ParametricNDSolve[eq /. {A -> 10, B -> 5, ω -> 2}, {x, v}, {t, 0, 20}, {R}];
ParametricPlot[{x[10000][t], v[10000][t]} /. s, {t, 4, 20},
AspectRatio -> 1, ImageSize -> Large]


The curve is symmetric about the origin. The bifurcation diagram requested in the question can be obtained from (code patterned after that in 96004)

tab = Table[{sol, points} = Reap@NDSolveValue[{v'[t] ==
x[t] - A x[t]^3 + R*Cos[ω*t] - B v[t] /. {A -> 10, B -> 5, ω -> 2},
x'[t] == v[t], x[0] == 0, v[0] == 0,
WhenEvent[v'[t] > 0, If[t > 4, Sow[x[t]]]]}, {x, v}, {t, 0, 20}]; {R, #} & /@
Union[Flatten[points], SameTest -> (Abs[#1 - #2] < .02 &)], {R, 25, 10000, 25}];
ListPlot[Flatten[tab, 1], ImageSize -> Large, AxesLabel -> {R, x}]