# How to make a bifurcation diagram of the Lorenz system under a varying parameter value?

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}

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

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 = InternalBag[]; (* 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, InternalStuffBag[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[InternalBagPart[res, All],
PlotStyle -> {Gray, Opacity[0.1], PointSize[0.001]},
AxesLabel -> {"r", "z"}]


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.

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

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

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

]
`