1
$\begingroup$

I'm trying to represent the behaviour of the system $\color{red}{r}$ varies from $\color{red}{1}$ to $\color{red}{100}$ statically (for instance something similar to a bifurcation graph)?

\begin{align} \frac{dx}{dt} &= \color{green}{10} \left(y-x\right) \\ \frac{dy}{dt}&= \color{red}{r}x -y - xz\\\ \frac{dz}{dt}&= xy -\color{green}{6}z \end{align}

How was this done?enter image description here

$\endgroup$
  • 1
    $\begingroup$ you have $\frac{dx}{dt}$ twice in there? Also what is $s$ and $b$? Do you have values for these? $\endgroup$ – Nasser Nov 3 '17 at 22:11
  • $\begingroup$ my bad. fixed! @Nasser $\endgroup$ – Conor Cosnett Nov 3 '17 at 22:28
  • 1
    $\begingroup$ If I'm understanding you right, you want a picture like this picture from here? $\endgroup$ – J. M. will be back soon Nov 4 '17 at 2:33
  • $\begingroup$ Dear @J.M. yes that is exactly what I want. What are they plotting on the vertical axis in that bifurcation graph? $\endgroup$ – Conor Cosnett Nov 4 '17 at 14:16
7
$\begingroup$

Looks like you want a plot of the extrema of one of the variables (maybe z[t] vs r). We can step through r values in a loop and save the extrema like this:

res = Internal`Bag[]; (* a place to store results *)

tmax = 100; (* how long to run for each r value *)
{x0, y0, z0} = {1, 1, 1}; (* initial ICs *)

Do[
  sol = NDSolve[{
    x'[t] == 10 (y[t] - x[t]),
    y'[t] == r x[t] - y[t] - x[t] z[t],
    z'[t] == x[t] y[t] - 6 z[t],
    x[0] == x0, y[0] == y0, z[0] == z0,
    (* save extrema of z[t] *)
    WhenEvent[z'[t] == 0, Internal`StuffBag[res, {r, z[t]}]]
   }, {x, y, z}, {t, 0, tmax}][[1]];

  (* save end value for next ICs *)
  {x0, y0, z0} = {x[tmax], y[tmax], z[tmax]} /. sol;
, {r, 200, 40, -0.1}];

ListPlot[Internal`BagPart[res, All],
  PlotStyle -> {Gray, Opacity[0.1], PointSize[0.001]},
  AxesLabel -> {"r", "z"}]

Mathematica graphics

A few notes:

  1. This takes a few minutes to run to get a nice looking result.

  2. Using final values from one run as initial conditions for the next is an easy way to stay near the attractor.

  3. Your value of b=6 is different than the b=8/3 used in the link, which is why the diagram is a little different.

  4. There may be alternative attractors for ranges of the parameter that this method will not find. It will also not find unstable limit sets.

$\endgroup$
  • $\begingroup$ So the axes are going to be labeled r and z? $\endgroup$ – Leif Oct 30 '18 at 14:17
  • $\begingroup$ @Leif Yes, I have labeled them now. $\endgroup$ – Chris K Mar 12 at 10:16
4
$\begingroup$

How about this? I just did direct solve and plot. Nothing fancy. Not sure if this is what you want or not. This can be made more use-friendly if needed

enter image description here

Manipulate[

 eq1 = x'[t] == s (y[t] - x[t]);
 eq2 = y'[t] == r x[t] - y[t] - x[t] z[t];
 eq3 = z'[t] == x[t] y[t] - b z[t];

 sol = NDSolveValue[{eq1, eq2, eq3, x[0] == x0, y[0] == y0, z[0] == z0}, 
       {x, y, z}, {t, 0, maxT}];

 ParametricPlot3D[{sol[[1]][t], sol[[2]][t], sol[[3]][t]}, {t, 0, maxT}, 
       AxesLabel -> {x, y, z}, 
       PlotRange -> All, 
       ImageSize -> {300, 300}],

 (*controls*)
 {{r, 15, "r"}, 1, 100, 1, Appearance -> "Labeled", ImageSize -> Tiny},
 {{s, 6, "s"}, 1, 100, 1, Appearance -> "Labeled", ImageSize -> Tiny},
 {{b, 10, "b"}, 1, 100, 1, Appearance -> "Labeled", ImageSize -> Tiny},
 {{maxT, 10, "time?"}, .1, 100, .1, Appearance -> "Labeled", ImageSize -> Tiny},

 Delimiter, (*initial conditions*)
 {{x0, 1, "x[0]"}, 0, 10, .1, Appearance -> "Labeled", ImageSize -> Tiny},
 {{y0, 0, "y[0]"}, 0, 10, .1, Appearance -> "Labeled", ImageSize -> Tiny},
 {{z0, 0, "z[0]"}, 0, 10, .1, Appearance -> "Labeled", ImageSize -> Tiny},
 ControlPlacement -> Left,
 ContinuousAction -> False

 ]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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