1
$\begingroup$

I’ve been trying to plot bifurcation diagram of 2D map following the example by Chris K in here, where I managed to plot few of periodic cycles

EDIT

(*2-cycle*)
A = 1.6; b = 2.47; c = 0.9; d = 1.19;
EM[{x_, y_}] := {x + y + A x y - b x, y + x/A - x y + c y^2 - d y}
EM2[{x_, y_}] = Simplify[EM[EM[{x, y}]]];
j2 := D[EM2[{x, y}], {{x, y}, 1}]
eqq2 = FindRoot[EM2[{x, y}] == {x, y}, {{x, 1.3}, {y, 0.644}}]
Eigenvalues[j2 /. eqq2] (*{-2.6888, 0.0128774} eigenvalues *)
ListLinePlot[Transpose[NestList[EM, {x + 10^-4, y} /. eqq2, 100]]]

(* 4-cycle *)
EM4[{x_, y_}] = Simplify[EM[EM[EM[EM[{x, y}]]]]];
j4 := D[EM4[{x, y}], {{x, y}, 1}]
eqq4 = FindRoot[EM4[{x, y}] == {x, y}, {{x, 1.3}, {y, 0.644}}]
Eigenvalues[j4 /. eqq4] (*{11.6416, 0.00122683} eigenvalues *)
ListLinePlot[Transpose[NestList[EM, {x + 10^-4, y} /. eqq4, 100]]]

Here is the 2-cycle

enter image description here

and 4-cycle

enter image description here

From these results, I try to code few attempts to analyse its bifurcation diagram with $b,c,d$ is fixed.

EDIT from Chris’s advice, I made some adjustment on Take : which now take last 50 elements from the list.

bifurcationPoints[AStart_, AStop_, m_, n_] := 
Flatten[Table [{A, #} & /@ Take[Quiet@NestList[{#1[[1]] 
+ #1[[2]] + A #1[[1]] #1[[2]] - 2.47 #1[[1]], 
#1[[2]] + #1[[1]]/A - #1[[1]] #1[[2]] + 0.9 #1[[2]]^2 - 1.2 #1[[2]]} &, 
{1.3, 0.6}, n][[All, 1]], -50], 
{A, AStart, AStop, (AStop - AStart)/m}], 1]
ListPlot[bifurcationPoints[1, 2.85, 270, 300]]

the code is working, but need to make some adjustment on the parameters.

However, my expectation is the plot ($x_t, y_t$ againts $A$) would be something like this or here.

Thus, from those 1D examples (with some modification), I tried

EM[A_, b_, c_, d_][{x_,y_}] := {x + y + A x y - b x, y + x/A - x y + c y^2 - d y}
pts = 2000;  
ListPlot[Flatten[Table[Transpose[{Table[A, pts],NestList[EM[A, 2.47, 0.9, 1.2], 
Nest[EM[A, 2.47, 0.9, 1.2], {1.3, 0.6}, 2000], pts - 1]}], {A,1.4, 2.0, 0.01}], 1], 
PlotStyle -> {Black, Opacity[0.2], PointSize[0.001]}, 
AxesLabel -> {A, xy}]

and from this example, I tried

eq1 = x[t] == x[t - 1] + y[t - 1] + A x[t - 1] y[t - 1] - 2.47 x[t - 1];
eq2 = y[t] == y[t - 1] + x[t - 1]/A - x[t - 1] y[t - 1] + 0.9 y[t - 1]^2 - 1.2 y[t - 1];
ic = {x[0] == 1.3, y[0] == 0.6};
eqtry = Flatten[Table[sol = RecurrenceTable[{eq1, eq2, ic}, {X[t], Y[t]}, {t, 0, 10}];
Replace[DeleteDuplicates[Take[sol, -2]], {{X_ -> {a, X}}, {Y_ -> {a, Y}}},1], {a, 1.4, 2.0, 0.1}], 1]      
ListPlot[eqtry]

The range of A could be from $(0, 1.5]$, but the code unable to produced any result. I believe this due to the selected parameters

Any suggestion anyone? Thanks in advance!

Improve this question
$\endgroup$
11
  • $\begingroup$ Could you include all the code you used in the first steps where you found cycles and eigenvalues? $\endgroup$
    – Chris K
    Commented Oct 20, 2020 at 16:31
  • $\begingroup$ Yes, definitely. It’s now edited to the original post. $\endgroup$
    – nightcape
    Commented Oct 20, 2020 at 22:55
  • 1
    $\begingroup$ Is there a reason you expect a period doubling cascade to occur? Do you expect chaos to emerge? From what I gather from a few minutes of playing with the bifurcation diagram, there's a whole lot of critical curves, but I don't see any chaotic regions. This of course doesn't mean there is no chaos; but you should first attempt an analytical investigation of the system. Moreover, is this system already known? Are there any papers in which it was studied? $\endgroup$
    – corey979
    Commented Oct 20, 2020 at 23:11
  • $\begingroup$ @corey979 Yeah, I did analytically analysed the system. Especially on the sensitivity of parameter A, which rise my suspicion on the chaotic behaviour. But no, this system is not known studied, and was motivated from predator-prey $\endgroup$
    – nightcape
    Commented Oct 20, 2020 at 23:20
  • $\begingroup$ Does this model have any physical basis? You mention that it was motivated by predator-prey dynamics, but those shouldn't become negative. $\endgroup$
    – Chris K
    Commented Oct 22, 2020 at 3:49

1 Answer 1

Reset to default
1
$\begingroup$

Some thoughts:

  1. The 4-cycle finding code does not seem to actually reach a 4-cycle because FindRoot gives errors.
  2. The Take in your first bifurcation code takes the first 15 time steps. I think you want to take the last ones.
  3. Here's a chaotic-looking result, taking the last 50 steps:
ListPlot[bifurcationPoints[1.53, 2.85, 270, 250]]

Mathematica graphics

  1. Your system seems to diverge for some range of $A<1.53$. For example,
A = 1.4;
ListLinePlot[Transpose[NestList[EM, {x + 10^-4, y} /. eqq2, 43]], 
 PlotRange -> All]

Mathematica graphics

This is why the bifurcation code doesn't work in this range. If I had to guess the source of the problem, it would be the term + c y^2 in the y equation, since $dy/dt=y^2$ diverges to infinity faster than exponentially.

Improve this answer
$\endgroup$
2
  • $\begingroup$ This is great! I guess I need to have A less than 3, even less than 1 if I want to keep the same parameters for b, c, d. $\endgroup$
    – nightcape
    Commented Oct 22, 2020 at 5:29
  • $\begingroup$ I honestly didn’t realised [Take] in the first bifurcation code. I’ll do the correction, thanks for pointing that out! $\endgroup$
    – nightcape
    Commented Oct 22, 2020 at 5:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.