Bifurcation with a system of equations

Question

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

Answer 1

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:

Answer 2

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