I have been trying to generate a bifurcation map but I am having trouble understanding my errors. I think my problem relies in "iterate" and "mps" since I copied the formatting from a wolfram demonstration code. I am just not familiar with making maps in mathematica.

Module[{deltad = 0.35, deltav = 2.3, deltaz = 0.35, alpha = 1,

uv = 5*10^10, \[Alpha] = 1, a = 1, gammay = 1,
ud = 1*10^6,
n = 3500,
s = 0.4,
r = 0.385,
\[Kappa] = 0.000001,
rho = 1,
sigmay = 0.9, sol, pts, iterate},

sol[\[Beta]_] := 
NDSolve[{x'[t] == 
  r*x[t] - \[Kappa]*((x[t]*z[t])/(
     x[t] + y[t] + z[t] + d[t])) - (\[Beta]*v[t]*x[t])/(
   x[t] + y[t] + z[t] + d[t]), 
 y'[t] == 
  s*y[t] + (\[Beta]*v[t]*x[t])/(x[t] + y[t] + z[t] + d[t]) - 
   alpha*y[t] - \[Kappa]*(y[t]*z[t])/(x[t] + y[t] + z[t] + d[t]), 
 v'[t] == 
  n*alpha*y[t] - deltav*v[t] + Piecewise[{{1*10^6, 3 <= t <= 5}}],
  d'[t] == 
  sigmay*y[t] - deltad*d[t] + Piecewise[{{5*10^10, 0 <= t <= 2}}],
  z'[t] == rho*d[t - \[Tau]] - deltaz*z[t], x[0] == 1.5*10^8, 
 y[0] == 0, v[0] == 0, d[t /; t <= 0] == th, z[0] == 0}, {x[t], 
 y[t], v[t], d[t], z[t]}, {t, 0, 310}];
iterate = 
Compile[{\[Beta]}, {fsol = sol[\[Beta]]; 
 Map[{\[Beta], #} &, 
  FindMaxValue[{x[t]} /. fsol, {t, #, # + 10}] & /@ 
   Table[i, {i, 200, 300.0, 10}]]}];
pts = Quiet[
  iterate[\[Beta]], {\[Beta], 0.00001, 0.005, 0.0005}], 1]]; 
PlotStyle -> 
Table[{PointSize[0.01], RGBColor[.49, 0, 0]}, {i, 1, 
  Length[pts]}], Frame -> True, ImageSize -> {400, 350}, 
PlotRange -> All, 
FrameLabel -> {Style["c", Italic], Style["z", Italic]}, 
ImageSize -> {500, 500}, AspectRatio -> 1, 
ImagePadding -> {{35, 10}, {35, 10}}]],
{{\[Tau], 0}, 0, 3}, {{th, 0}, 0, 2}, ControlPlacement -> Top, 
SynchronousUpdating -> False]

Here is the link :http://demonstrations.wolfram.com/BifurcationDiagramForTheRoesslerAttractor/. The code is on the corner. Thank you!

  • $\begingroup$ Just added the link. $\endgroup$
    – peace96
    Jul 27, 2017 at 21:44
  • $\begingroup$ Do you need the Manipulate? What do you want on the x- and y-axes of your bifurcation diagram? $\endgroup$
    – Chris K
    Jul 27, 2017 at 21:58
  • $\begingroup$ I've had a quick play with your ODE (ignoring iterate and pts for the moment) and, for the ranges you give for \[Tau] and th, can't find anything looking like an equilibrium or limit cycle. What kinds of bifurcations are you expecting to see in your system? I may well be missing something, but iterate (and the appearance of FindMaxValue in it) doesn't look like it would produce what I would think of as a bifurcation diagram. $\endgroup$ Jul 28, 2017 at 0:06
  • $\begingroup$ I am trying to vary beta (X axis) vs x[t]+y[t]. The manipulate kinda helps me to see if the system changes as we vary the delay. If beta=0.00341, I get a a limit cycle / steady states for x[t]+y[t] (time=200). I was kinda thinking that iterate isn't what I should be using. Any way of approaching this better? $\endgroup$
    – peace96
    Jul 28, 2017 at 0:40
  • $\begingroup$ Okay, that makes a lot more sense. iterate seems to be finding a list of maximal amplitudes, which seems a bit odd. And it looks like the oscillations are transient. But there's definitely bifurcations in there. $\endgroup$ Jul 28, 2017 at 1:05

1 Answer 1


Note: This is not really an answer to the original question, and as such is a little off-topic. I've tried to give a bit of context (and code) to show how to make a start approaching bifurcations with Mathematica.

First off: If you're new to bifurcation theory I would recommend not starting here. It's got forcing, delay, branches of equilibria going off to infinity, Hopf bifurcations, etc. It's pretty messy. What I will do is show you (a) how I'd first approach a bifurcation problem, and (b) where this one starts getting ugly. It will, unfortunately, be a very brief sketch because I don't have all that much time.

For simplicity I'm going to ignore the forcing terms on v and d, and the delay on d (the initial forcing won't make any difference to the asymptotic behaviour anyway). Using your parameter values define F by

F[x_, y_, v_, d_, z_, \[Beta]_] := {
  r*x - \[Kappa]*((x*z)/(x + y + z + d)) - (\[Beta]*v*x)/(x + y + z + d),
  s*y (\[Beta]*v*x)/(x + y + z + d) - alpha*y - \[Kappa]*(y*z)/(x + y + z + d),
  n*alpha*y - deltav*v,
  sigmay*y - deltad*d,
  rho*d - deltaz*z}

Then we can find the (unique) branch of equilibria of $du/dt = F(u,\beta)$ (where $u = (x, y, v, d, z)$) as a function of $\beta$ by

eqsol = FullSimplify[
  Solve[F[x, y, v, d, z, \[Beta]] == 0, {x, y, v, d, z}][[1]]]

$$ {x \rightarrow 6.49351 + \frac{4.34205*10^{-8}}{\beta},\\ y \rightarrow -\frac{0.00164286}{0.00276234 - \beta} - \frac{3.97684*10^{-9}}{\beta},\\ v \rightarrow -\frac{2.50001}{0.00276234 - \beta} - \frac{6.05172*10^{-6}}{\beta}, \\ d \rightarrow -\frac{0.0042245}{0.00276234 - \beta} - \frac{1.02262*10^{-8}}{\beta}, \\ z \rightarrow -\frac{0.01207}{0.00276234 - \beta} - \frac{2.92176*10^{-8}}{\beta}} $$

  Plot[# /. eqsol, {\[Beta], 10^-5, 0.005}, AxesLabel -> {\[Beta], #}] & /@ {x, y, v, d, z}, 

enter image description here

Clearly this branch of equilibria is not very well behaved, and is undefined at $\beta = 0.00276234$.

Pressing on, let's consider the stability of these equilibria, and in particular, look for values of $\beta$ at which one or more eigenvalues change sign (indicating that a bifurcation occurs).

The general form of the linearized system (the Jacobian) is D[F[x, y, v, d, z, \[Beta]], {{x, y, v, d, z}}]. Evaluating on the branch of equilibria:

Leq = Simplify[D[F[x, y, v, d, z, \[Beta]], {{x, y, v, d, z}}] /. eqsol];

Mathematica has trouble finding the eigenvectors of Leq, but I'm not going to dwell on it here. It can find the eigenvalues, though, which are more important for getting a rough sketch of what's going on.

evals = Simplify[Eigenvalues[Leq]];

Plot the real part of the eigenvalues against $\beta$ (the last two eigenvalues are a complex conjugate pair):

GraphicsGrid[Partition[Plot[Re[#], {\[Beta], 10^-5, 0.005}] & /@ evals, UpTo[3]]]

enter image description here

Again, nasty (but not unexpected) stuff going on at $\beta = 0.00276234$. The first thing to notice is that there's always at least one eigenvalue with positive real part, so this branch of equilibria is always unstable. (This may or may not be of interest to you -- sometimes existence of solutions is the only thing that matters). However, understanding how other solutions branch off this one can still shed light on what you might observe in numerical solutions. The main thing to note about the system is that there's a Hopf bifurcation as the real part of the complex eigenvalues 4 and 5 change signs at around $\beta=0.001$. You can pin this down with

FindRoot[Re[evals[[4]]], {\[Beta], 0.001}]
FindRoot[Re[evals[[5]]], {\[Beta], 0.001}]

{[Beta] -> 0.00105973}

{[Beta] -> 0.00105973}

That is, there is (generically) a branch of periodic solutions with amplitude growing as $\sqrt{\beta-\beta_0}$. This may or may not be stable -- the eigenvalues I found earlier are for the branch of equilibria and so don't apply here -- but it's probably going to be unstable close to the bifurcation point.

Summary: Without doing an in-depth analysis, I really can't say a whole lot more. In general, a bifurcation analysis would start by identifying a reference state -- a branch of equilibria to study bifurcations from. Typically, after a change of coordinates, this branch would be at the origin for all parameter values, ie: $F(0, \lambda) \equiv 0$, and referred to as the "trivial branch". Then you look at the linearized system on the trivial branch and compute eigenvalues and eigenvectors. A bifurcation occurs when the real part of an eigenvalue changes sign. Generically (in the absence of symmetry or other degeneracies) only one eigenvalue (or a complex conjugate pair in a Hopf bifurcation) will change sign at a time, so the nullspace of the linearization at the bifurcation point will be one-dimensional (or one complex dimension). Then you might project the bifurcation onto the nullspace or the centre manifold to reduce the dimensions in a large system. That's usually what gets plotted in a bifurcation diagram -- the projection of the bifurcating branch onto one-dimensional subspace (it's complicated for Hopf, but you can still get it down to one dimension). Then it's an orgy of analysis and you follow wherever it takes you.

I don't know what the physical interpretation of your system is, but I'm unsure of how well it can be describing a physical system. There are, as far as I can tell, no stable equilibria anywhere. Does the unstable branch of equilibria have a convincing explanation? It should, before you start any in-depth analysis. For a lot of parameters in the ranges you give, numerical solutions tend to blow up to infinity. Does that makes sense? My point is that, if I were you, I would makes very sure that the system is doing what it's supposed to be doing before embarking on an in-depth bifurcation analysis. The features shown up by this preliminary study must make sense (or, ideally, be expected before you even do the calculations). Interesting bifurcation analyses never proceed along the "typical" lines I mentioned above -- that's what makes them interesting. But I've gone down far too many pointless rabbit holes that turn out to be dead ends because I didn't get things solid before I started.

If you're looking for references to get going with, you can't go far wrong with Steve Strogatz's "Nonlinear Dynamics and Chaos". I'm partial to the singularity approach in Golubitsky and Schaeffer "Singularities and groups..., vol 1", but that's just a personal preference.

Good luck. Bifurcation theory is awesome.

Edit: Hopf bifurcations. Some authors do indeed define a Hopf bifurcation as having all other eigenvalues negative, but it is by no means necessary and a number of things get conflated to confuse the issue.

  1. Examples of Hopf bifurcation usually have other eigenvalues negative so they can demonstrate the numerical behaviour of the system.
  2. Some authors put negative eigenvalues in the definition of the Hopf bifurcation (like in the Wikipedia page). This ensures that the branch of periodic orbits will be observable in numerical simulations because of the change of stability at the bifurcation point. So it's kind of just an extension of point 1.
  3. Some texts are talking about the normal form of the Hopf bifurcation (which any Hopf bifurcation is locally equivalent to on the centre manifold), in which case all other eigenvalues are irrelevant.

If the branch of equilibria has any positive eigenvalues then neither the equilibria nor the periodic orbits close to the bifurcation point will be observable in numerical simulations. That may make them uninteresting in many contexts. But there is still a branch of (unstable) periodic orbits. And a Hopf bifurcation, in its most general form, is concerned primarily with the existence of this branch of solutions, not its stability.

Two things that will mess up a Hopf bifurcation:

  1. Other eigenvalues on the imaginary axis.
  2. The complex pair of eigenvalues not crossing the imaginary axis with nonzero speed; ie, the real parts just touch the axis, but don't actually cross.
  • $\begingroup$ Is there really a Hopf bifurcation around β=0.001? The first eigenvalue has positive real part there, so the equilibrium is already unstable. I believe Hopf bifurcation requires all other eigenvalues to be negative. $\endgroup$
    – Chris K
    Jul 30, 2017 at 17:47
  • $\begingroup$ @ChrisK See my edit. It's appears to be a matter of taste how to define a Hopf bifurcation, with authors choosing a form that suits their purpose. $\endgroup$ Jul 30, 2017 at 23:04
  • $\begingroup$ Thanks, I wasn't aware of this possibility! $\endgroup$
    – Chris K
    Jul 31, 2017 at 1:59

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.