Chaotic Dynamical Systems Tools

Bifurcation diagrams with Mathematica

Jean-Marc GINOUX

ginoux@univ-tln.fr

Bruno ROSSETTO

rossetto@univ-tln.fr

The aim of this work is to present applets elaborated with Mathematica in order to study two or three-dimensional autonomous dynamical systems depending on parameters. As the parameters vary, the phase portrait of such systems also varies. So, the appearance of a topologically nonequivalent phase portrait under the variation of parameters is called a bifurcation. More precisely, the Hopf-bifurcation plays a very important role in the study of dynamical systems. Hence, the first applet allows to compute the Hopf bifurcation parameter value of two or three-dimensional autonomous dynamical systems with a high accuracy. As pointed out by Glass and Mackey, the construction of a bifurcation diagram is a good means of locating the signature of chaos in a dynamical system. The second applet makes it possible to construct the bifurcation diagram of such systems. Thus, chaotic behaviours such as period doubling cascade may be highlighted. Various examples of applications are proposed to illustrate the usefulness of these applets.

Characterizing the evolution of a given physical phenomenon such as population dynamics, for example, may be done thanks to a mathematical model which may consist of system of differential equations involving state variables and parameters. If the variables of such system do not depend explicitly on time, the system is said to be autonomous. Under certain assumptions defined below such system may be considered as autonomous dynamical system. As time does not occur explicitly in equations, solution of a system of differential equations may be projected in a space called phase-space in which the behaviour of the state variables is described. The plot of the solution in such space is called phase portrait. The bifurcation of a system of differential equation, i.e., of autonomous dynamical system is concerned with changes in the qualitative behaviour of its phase portrait as parameters vary and more precisely, when such a bifurcation parameter reaches a certain value, called critical value. Thus, bifurcation theory is of great importance in dynamical systems study because it indicates stability changes, structural changes in a system etc... So, plotting the solution of autonomous dynamical system according to the bifurcation parameter leads to the construction of a bifurcation diagram. Such diagram provides knowledge on the behaviour of the solution: constant, periodic, nonperiodic or even chaotic as pointed out by Glass and Mackey. As there is many kinds of behaviour of solutions there is many kinds of bifurcations. In this work, we focus on Hopf bifurcation corresponding to periodic solutions and period doubling bifurcation or period doubling cascade which is one of the route to chaos for dynamical systems.

In the following we consider a system of differential equations defined in a compact U included in (ℝ^n):

\(\(d  X\&⟶\)\/dt = \(\(ℱ\&⟶\) \((X\&⟶, μ)\)\(\ \)\)\)

with

(Overscript[X, ⟶])=(([(Subscript[x, 1]), (Subscript[x, 2]),  ... , 
                    (Subscript[x, n])])^t) ϵ U⊂(ℝ^n), μ ϵ V⊂(ℝ^p)

and

\(\(ℱ\&⟶\) \((X\&⟶, μ)\)\)=\(\(\([\(f\_1\) \((X\&⟶, μ)\), \(f\_2\) \((X\&⟶, μ)\),  ... , \(f\_n\) \((X\&⟶, μ)\)]\)\^t\) ϵ\ U\)⊂(ℝ^n)

The vector (Overscript[ℱ, ⟶]) defines a velocity vector field in U whose components (Subscript[f, i]) which are supposed to be continuous and infinitely derivable with respect to all (Subscript[x, i]) and t, i.e., are (C^∞) functions in U and with values included in ℝ, check the assumptions of the Cauchy-Lipschitz theorem. For more details, see for example Coddington & Levinson, [1955]. A solution of this system is an integral curve (Overscript[X, ⟶])(t) tangent to (Overscript[ℱ, ⟶]) whose values define the states of the dynamical system described by the Eq. (1). Since none of the components (Subscript[f, i]) of the velocity vector field depends here explicitly on time, the system is said to be autonomous.

Let's consider a dynamical system defined under the same conditions as above and having a fixed point \(\(X\&⟶\) = \(\(X\&⟶\)\^\(\(\ \)\(*\)\)\) \((μ)\)\)and let's make the following assumptions.

The jacobian functional matrix associated with this system has a pair of complex-conjugate eigenvalues \(σ\_\(1, 2\) = α \((μ)\) + \(\(i\)\(\ \)\(ω\) \((μ)\)\(\ \)\)\)such that:

- for a critical value of the bifurcation parameter μ = (Subscript[μ, C]), α ((Subscript[μ, C])) = 0 and (Subscript[((d(α/d)μ)), μ  =  (Subscript[μ, C])]) ≠ 0
- the n - 2 others eigenvalues evaluated in μ have their real parts strictly negative. 

If these assumptions are checked, then \(\((\(X\&⟶\^\(\(\ \)\(*\)\)\) \((μ\_C)\), μ\_C)\)\)is a Hopf bifurcation point which leads to the appearance, from the equilibrium state ((Overscript[X, ⟶])^( *)) ((Subscript[μ, C])), of a limit cycle.

Computation of the Hopf bifurcation parameter value

Using the following procedure it is possible to compute the Hopf bifurcation parameter value of two or three-dimensional dynamical systems. Since the two-dimensional procedure may be obtained by a simple reduction, the three-dimensional procedure is only presented.

Dynamical system

In dimension three the autonomous dynamical system (1) may be defined under the same conditions as:

\(\(d  X\&⟶\)\/dt = \(\(ℱ\&⟶\) \((X\&⟶)\)\(\ \)\)\)

with

(Overscript[X, ⟶])=(([x, y, 
                    z])^t) ϵ U⊂(ℝ^3), (λ, μ, ρ) ϵ V⊂(ℝ^3)

and

\(\(ℱ\&⟶\) \((X\&⟶)\)\)=(([f ((Overscript[X, ⟶]), λ), 
                g ((Overscript[X, ⟶]), μ), 
                h ((Overscript[X, ⟶]), ρ)])^t) ϵ U⊂(ℝ^3)

where λ, μ and ρ are real parameters.

((d  (Overscript[X, ⟶]))/dt)(({{((d x)/dt)}, {((d y)/dt)

                ((d z)/dt)}}))=(Overscript[ℱ, ⟶])(({{f (x, y, z, λ)}, {g (x, y, z, μ)
h(x,y,z,ρ)}}))

Bifurcation parameter

Choose a bifurcation parameter between the set of parameters of the dynamical system (2). λ for example.

Fixed points

Choosing λ as the Hopf bifurcation parameter and while "freezing" all the others, i.e., μ and ρ, each fixed point may be expressed according to the parameter λ.

Solve[{f (x, y, z, λ) == 0, g (x, y, z, μ) == 0, 
        h (x, y, z, ρ) == 0}, {x, y, z}]

Choose the fixed point in the vicinity of which the Hopf bifurcation may occurs. Let's call I this fixed point expressed according to the parameter λ.

Functional jacobian matrix

Since the functions (Subscript[f, i]) are supposed to be (C^∞) functions in a compact U included in (ℝ^n), it is possible to compute the functional jacobian matrix 𝒥 of system (1) or (2), defined by: 𝒥 = ((d  (Overscript[ℱ, ⟶]))/(d  
                  (Overscript[X, ⟶]))). The following procedure provides each elements of the functional jacobian matrix.

(Subscript[a, 11]) = Simplify[D[f, x]]

(Subscript[a, 12]) = Simplify[D[f, y]]

(Subscript[a, 13]) = Simplify[D[f, z]]

(Subscript[a, 21]) = Simplify[D[g, x]]

(Subscript[a, 22]) = Simplify[D[g, y]]

(Subscript[a, 23]) = Simplify[D[g, z]]

(Subscript[a, 31]) = Simplify[D[h, x]]

(Subscript[a, 32]) = Simplify[D[h, y]]

(Subscript[a, 33]) = Simplify[D[h, z]]

J = {{(Subscript[a, 11]), (Subscript[a, 12]), (Subscript[a, 13])}, {(Subscript[a, 21]), (Subscript[a, 22]), (Subscript[a, 23])}, {(Subscript[a, 31]), (Subscript[a, 32]), 
          (Subscript[a, 33])}}

MatrixForm[J]

Plugging expression of the fixed point I in the functional jacobian matrix and while using the Mathematica function Eigenvalues, it is possible to compute the eigenvalues according to the chosen parameter λ. Let's denote Sp (for spectrum) these eigenvalues.

Sp:=Eigenvalues[J];

Graphical estimation of the Hopf bifurcation parameter value

In order to find the range of this parameter, the real part of one of the complex-conjugate eigenvalue is plotted (with Plot function) according to the parameter chosen (λ). Then, the location of the points where this function vanishes provides an estimation of the Hopf bifurcation parameter value.

Plot[Re[Sp[[2]]],{λ, (Subscript[x, min]), (Subscript[x, max])},PlotRange -> {(Subscript[y, min]), (Subscript[y, max])}]

It provides the approximate value: λ\(TraditionalForm\`\_estimate\)

To confirm this, the other real part of the other complex-conjugate eigenvalue is also plotted according to the parameter chosen (λ).

Plot[Re[Sp[[3]]],{λ, (Subscript[x, min]), (Subscript[x, max])},PlotRange -> {(Subscript[y, min]), (Subscript[y, max])}]

Numerical computation of the Hopf bifurcation parameter value

Then, the Mathematica function: FindRoot, which uses a damped Newton's method, the secant method and Brent's method, provides a numerical computation of the Hopf bifurcation parameter value with a high accuracy depending on the choice of the operator and on the WorkingPrecision.

FindRoot[Re[Sp[[2]]]==0,{λ,λ\(\_estimate\)}]

This provides the Hopf bifurcation parameter value starting from the vicinity of its graphical estimation:  λ\(\_estimate\)

Then, it can be confirmed also for the other real part of the other complex conjugate eigenvalue.

FindRoot[Re[Sp[[3]]]==0,{λ,λ\(\_estimate\)}]

Example: Hopf bifurcation in the Lorenz attractor

The purpose of the model established by Edward Lorenz [1963] was in the beginning to analyze the impredictible behaviour of weather. After having developed non-linear partial derivative equations starting from the thermal equation and Navier-Stokes equations, Lorenz truncated them to retain only three modes. The most widespread form of the Lorenz model is as follows:

(Overscript[V, →])(({{((d x)/dt)}, {((d y)/dt)

                ((d z)/dt)}}))=(Overscript[ℱ, ⟶])(({{f (x, y, z, σ)}, { g (x, y, z, r)

                  h (x, y, z, β)}}))=(({{σ (y - x)}, {-x z + r x - y

                x y - β z}}))

with σ, r, and β are real parameters:  σ = 10, β=(8/3), r = 28

where the functions f, g and h are infinitely derivable with respect to all (Subscript[x, i]), and t, i.e., are (C^∞)functions in a compact U included in (ℝ^3) and with values in ℝ. Lorenz model has been extensively studied and the Hopf bifurcation parameter value r in the vicinity of the fixed points labelled \(\(C\^±\) \((\(±\@\(β \((r - 1)\)\)\), \(±\@\(β \((r - 1)\)\)\), r - 1)\)\)is well-known:

r = (Subscript[r, C]) = σ (( σ + β + 3)/(σ - β - 1)) ≃ 24.7368

So, let's apply the above procedure.

Solve[{σ (y - x) == 0, -x z + r x - y == 0, 
        x y - β z == 0}, {x, y, z}]

\({{y -> 0, z -> 0, 
        x -> 0}, {y -> \(-\@\(\(-β\) + r\ β\)\), 
        z -> \(-1\) + r, 
        x -> \(-\@\(\(-β\) + r\ β\)\)}, {y -> \@\(\(-β\) + r\ β\), z -> \(-1\) + r, 
        x -> \@\(\(-β\) + r\ β\)}}\)

Fixing parameters σ and β, one seeks the Hopf bifurcation parameter value r in the vicinity of any fixed point (C^±). Let's choose arbitrary \(\(C\^-\)\) and compute the functional jacobian matrix.

J=(({{-10, 10, 0}, {r - z, -1, -x}, {y, x, -8/3}}))

Evaluated at \(\(C\^-\)\) the functional jacobian matrix takes the form:

J=(({{-10, 10, 0}, {1, -1, 2 (Sqrt[(2/3)]) (Sqrt[-1 + r])}, {-2 (Sqrt[(2/3)]) (Sqrt[-1 + r]), -2 (Sqrt[(2/3)]) (Sqrt[-1 + r]), -(8/3)}}))

Then, it is possible to compute its eigenvalues according to the chosen parameter r.

Sp:=Eigenvalues[J];

Now, let's plot the real part of one of the complex-conjugate eigenvalues according to r.

Plot[Re[Sp[[2]]],{r, 0, 30},PlotRange -> {-3, 3},\(AxesLabel -> {"r", \*""Re[(Subscript[λ, 2])]""}\),PlotStyle -> {RGBColor[1, 0, 0]}]

The other provides the following plot.

Plot[Re[Sp[[3]]],{r, 0, 30},PlotRange -> {-3, 3},\(AxesLabel -> {"r", \*""Re[(Subscript[λ, 3])]""}\),PlotStyle -> {RGBColor[1, 0, 0]}]

Then, the location of the points where these functions (Re[(Subscript[λ, 1, 2])]=ϕ(r)) vanish provides an estimation of the Hopf bifurcation parameter value, i.e., r ≃ 24. Using FindRoot function provides a numerical computation of the Hopf bifurcation value r.

FindRoot[Re[Sp[[2]]]==0,{r, 24}]

{r -> 24.736842105263143`}

Of course the other provides the same result.

FindRoot[Re[Sp[[3]]]==0,{r, 24}]

{r -> 24.736842105263143`}

Both values of r are corresponding to the one theorically computed.

Since it has been established [Ruelle et al., 1971] that chaos may occur in three-dimensional dynamical systems, analytical methods such as the existence of Shil'nikov orbit [Shil'nikov, 1965] and numerical tools such as the construction of a bifurcation diagram [Glass et al., 1988] have been developed. Let's consider as previously, a three-dimensional autonomous dynamical system defined under the same conditions as above.

\(\(d  X\&⟶\)\/dt = \(\(ℱ\&⟶\) \((X\&⟶)\)\(\ \)\)\)

with

(Overscript[X, ⟶])=(([x, y, 
                    z])^t) ϵ U⊂(ℝ^3), (λ, μ, ρ) ϵ V⊂(ℝ^3)

and

\(\(ℱ\&⟶\) \((X\&⟶)\)\)=(([f ((Overscript[X, ⟶]), λ), 
                g ((Overscript[X, ⟶]), μ), 
                h ((Overscript[X, ⟶]), ρ)])^t) ϵ U⊂(ℝ^3)

where λ, μ and ρ are real parameters.

((d  (Overscript[X, ⟶]))/dt)(({{((d x)/dt)}, {((d y)/dt)

                ((d z)/dt)}}))=(Overscript[ℱ, ⟶])(({{f (x, y, z, λ)}, {g (x, y, z, μ)
h(x,y,z,ρ)}}))

Choose the bifurcation parameter (λ for example) between the set of parameters of the dynamical system and "freeze" all the others (μ and ρ). In order to construct a bifurcation diagram, the system is numerically integrated (with NDSolve function) according to the bifurcation parameter λ. Then, each variable x[t], y[t] and z[t], solution of this dynamical system (3) may be numerically computed according to the bifurcation parameter λ. A bifurcation diagram is the plot of the successive maxima of one of these variables (z[t], for example) as a function of the bifurcation parameter λ. The following applet developed by Eric Javoy, computer scientist, is based on dichotomy principle and uses a follower which is "tracking" the curve until it reaches its maximum. This applet provides the bifurcation diagram of a three-dimensional dynamical system. It plots the maximum (Subscript[Z, max]) of the solution of the third variable Z of this dynamical system in function of the bifurcation parameter chosen λ. The bifurcation parameter λ is in the range: \(λ\_min\  ≤ \ λ\  ≤ \ λ\_max\). The step (= Pasλ) of evaluation may be adjusted in function of the region of the diagram: \(λ < λ\_int1, 
      0.0001, \((λ ≥ λ\_int1)\) && \((λ < λ\_int2)\), 0.0001, λ ≥ λ\_int2, 
      0.0001\).
The dynamical system is first integrated with a classical method of NDSolve between the range of time: {t, (Subscript[t, 0]), (Subscript[t, max])}.
Then, the parameter is varying with the chosen step (= Pasλ) and Z is evaluated in order to find its maximum between (Subscript[t, min]) and (Subscript[t, max]) after elimination of the transient behaviour.

Applet by Eric Javoy, computer scientist

Clear["Global`*"]
Off[General::spell1]
$RecursionLimit = Infinity;

μ=;
ρ=;

Pasλ = 
        0.001;
TableauPlot = {};

For[λ = (Subscript[λ, min]),λ ≤ (Subscript[λ, max]), λ = λ + Pasλ,
solution=NDSolve[

{x'[t]==f(x, y, z, λ),x[0] == (Subscript[x, 0]),
y'[t]==g (x, y, z, μ),y[0] == (Subscript[y, 0]),
z'[t]==h(x,y,z,ρ),z[0] == (Subscript[z, 0])},

{x, y, z},{t, (Subscript[t, 0]), (Subscript[t, max])},Method -> ExplicitRungeKutta];

Pasλ = 
            Which[λ < (Subscript[λ, int1]), 
              0.0001, (λ ≥ (Subscript[λ, int1])) && (λ < (Subscript[λ, int2])), 
              0.0001, λ ≥ (Subscript[λ, int2]), 0.0001];

Print["λ=", λ];
Pas =  .05;
Tendance =  0;
BorneBasse = (Subscript[t, min]);
BorneHaute = (Subscript[t, max]);
Precedent = Evaluate[{z[BorneBasse]} /. 
                    solution][[1]][[1]];
Tendance = 
            Which[Precedent >= Evaluate[{z[BorneBasse + Pas]} /. 
                        solution][[1]][[1]], -1, 
              Precedent < Evaluate[{z[BorneBasse + Pas]} /. 
                        solution][[1]][[1]], 1];
\(For\ [Valeur = BorneBasse + Pas, \ 
            Valeur ≤ BorneHaute, 
            Valeur = 
              Valeur + 
                Pas, 
\ \ \ \ \ \ \ \ Actuel = \(\(Evaluate[{z[Valeur]} /. 
                      solution]\)[\([1]\)]\)[\([1]\)]; 
\ \ \ \ \ \ \ \ If[\((Actuel >= Precedent)\), 
              Goto\ [NonMaximum]]; 
\ \ \ \ \ \ \ \ If[
              Actuel < Precedent\  && \ Tendance == \(-1\), 
              Goto[ActualisePrecedent]]; 
\ \ \ \ \ \ \ \ TableauPlot = 
              Append[TableauPlot, {N[FromDigits[RealDigits[λ]], 6], 
                  N[FromDigits[
                      RealDigits[\(\(Evaluate[{z[Valeur - Pas]} /. 
                                solution]\)[\([1]\)]\)[\([1]\)]]], 
                    6]}]; 
Tendance = \(-1\); 
\ \ \ \ \ \ \ \ Goto\ [
              ActualisePrecedent]; \ \ \ \ \ \ \ 
Label[
              NonMaximum]; 
\ \ \ \ \ \ \ \ Tendance = 
              1; 
\ \ \ \ \ \ \ \ Goto[
              ActualisePrecedent]; \ \ \ \ 
Label[
              RechercheFine]; 
\ \ \ \ \ \ \ \ Print["t="; Valeur]; \ \ \ 
Label[
              PremieresValeurs]; 
\ \ \ \ \ \ \ \ If[
              Valeur != BorneBasse, 
              Tendance = 
                Which[Actuel < Precedent, \(-1\), 
                  Actuel\  ≥ Precedent, 
                  1]]; 
Label[
              ActualisePrecedent]; \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ Precedent = Actuel; 
Label[
              FinTraitement];
\ \ \ \ \ \ \ ]\)
]

The following function plots the bifurcation diagram.

\(ListPlot[TableauPlot, PlotStyle -> {PointSize[0.005], Hue[ .6]}, 
      ImageSize -> {400, 400}, 
      PlotRange -> {{λ\_min, λ\_max}, {z\_a, z\_b}}, 
      Frame -> True, 
      FrameLabel -> {"λ ", \*""\!\(Z\_max\)""}, 
      RotateLabel -> False]\)

Example: bifurcation diagram for the Volterra-Gause model

This three-dimensional predator-prey model, elaborated by Ginoux et al. [2005], consisted of a prey, a predator and top-predator has been named Volterra-Gause because it combines the original model of V.Volterra [1926] incorporating a logisitic limitation of P.F.Verhulst [1838] type on the growth of the prey and a limitation of G.F.Gause [1935] type on the intensity of the predation of the predator on the prey and of top-predator on the predator.

(Overscript[V, →])(({{((d x)/dt)}, {((d y)/dt)

                ((d z)/dt)}}))=(Overscript[ℱ, ⟶])(({{f (x, y, z, ξ)}, { g (x, y, z, (Subscript[δ, 1]))

                  h (x, y, z, ε, (Subscript[δ, 2]))}}))=(({{(1/ξ) (x (1 - x) - (x^((1/2))) y)}, {-(Subscript[δ, 1]) y + (x^((1/2))) y - (y^((1/2))) 
                    z

                ε z ((y^((1/2))) - (Subscript[δ, 2]))}}))

with ξ, ε, (Subscript[δ, 1]) and (Subscript[δ, 2]) are real parameters: ξ = 0.866, ε = 1.428, (Subscript[δ, 1]) = 0.577, (Subscript[δ, 2]) = 0.376 and where the functions f, g and h are infinitely derivable with respect to all (Subscript[x, i]), and t, i.e., are (C^∞)functions in a compact U included in (ℝ^3) and with values in ℝ.  Variations of the parameter (Subscript[δ, 1]) induce a period doubling cascade highlighted in the bifurcation diagram and which transcribes the chaotic feature of this model. Thus, choosing (Subscript[δ, 1]) as the bifurcation parameter and "freezing" all the others with the corresponding values: ξ = 0.866, ε = 1.428, (Subscript[δ, 2]) = 0.376, the bifurcation diagram of this system may be built with the above applet.

ξ = 0.8660254037844386`;
ε = 1.4288690166235205`;
(Subscript[δ, 2]) = 0.37602114655785496`;
k = 0.5;
p = 0.5;
Pasλ = 0.001;
TableauPlot = {};

\(For[λ = 0.6, λ ≤ 
        0.67, \ λ = λ + 
          Pasλ, 
solution = 
        NDSolve[{ξ\ \(x'\)[t] == \ x[t] - x[t]^2 - 0.25\ x[t]^k\ y[t], 
            x[0] == 
              0.15, \(y'\)[t] == \(-λ\)\ y[t] + x[t]^k\ y[t] - 
                0.25^0.5\ y[t]^p\ z[t], 
            y[0] == 
              0.5, \(z'\)[
                t] == ε\ \((0.25^0.5\ y[t]^
                        p - δ\_2)\)\ z[t], 
            z[0] == 0.5}, 
{x, y, z}, {t, 0, 700}, 
          Method -> 
            ExplicitRungeKutta]; 
Pasλ = 
        Which[λ < 0.628, 
          0.0001, \((λ ≥ 0.628)\) && \((λ < 
                0.634)\), 0.0001, λ ≥ 0.634, 
          0.0001]; 
Print["λ=", λ]; 
Pas =  .05; 
Tendance = \ 
        0; 
BorneBasse = 
        600; 
BorneHaute = 
        700; 
Precedent = \(\(Evaluate[{z[BorneBasse]} /. 
                solution]\)[\([1]\)]\)[\([1]\)]; 
Tendance = Which[Precedent >= \(\(Evaluate[{z[BorneBasse + Pas]} /. 
                    solution]\)[\([1]\)]\)[\([1]\)], \(-1\), 
          Precedent < \(\(Evaluate[{z[BorneBasse + Pas]} /. 
                    solution]\)[\([1]\)]\)[\([1]\)], 
          1]; 
For\ [Valeur = BorneBasse + Pas, \ 
        Valeur ≤ BorneHaute, 
        Valeur = 
          Valeur + 
            Pas, 
\ \ \ \ \ \ \ \ Actuel = \(\(Evaluate[{z[
                      Valeur]} /. 
                  solution]\)[\([1]\)]\)[\([1]\)]; 
\ \ \ \ \ \ \ \ If[\((Actuel >= Precedent)\), 
          Goto\ [NonMaximum]]; 
\ \ \ \ \ \ \ \ If[
          Actuel < Precedent\  && \ Tendance == \(-1\), 
          Goto[ActualisePrecedent]]; 
\ \ \ \ \ \ \ \ TableauPlot = 
          Append[TableauPlot, {N[FromDigits[RealDigits[λ]], 6], 
              N[FromDigits[
                  RealDigits[\(\(Evaluate[{z[Valeur - Pas]} /. 
                            solution]\)[\([1]\)]\)[\([1]\)]]], 
                6]}]; 
Tendance = \(-1\); 
\ \ \ \ \ \ \ \ Goto\ [
          ActualisePrecedent]; \ \ \ \ \ \ \ 
Label[
          NonMaximum]; 
\ \ \ \ \ \ \ \ Tendance = 
          1; 
\ \ \ \ \ \ \ \ Goto[
          ActualisePrecedent]; \ \ \ \ 
Label[
          RechercheFine]; 
\ \ \ \ \ \ \ \ Print["t="; 
          Valeur]; \ \ \ 
Label[
          PremieresValeurs]; 
\ \ \ \ \ \ \ \ If[
          Valeur != BorneBasse, 
          Tendance = 
            Which[Actuel < Precedent, \(-1\), 
              Actuel\  ≥ Precedent, 
              1]]; 
Label[
          ActualisePrecedent]; \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ Precedent = Actuel; 
Label[
          FinTraitement];
\ \ \ \ \ \ \ ]
]\)

\(ListPlot[TableauPlot, PlotStyle -> {PointSize[0.005], Hue[ .6]}, 
      ImageSize -> {400, 400}, 
      PlotRange -> {{0.6, 0.67}, {0.5, 1.1}}, Frame -> True, 
      FrameLabel -> {\*""(Subscript[δ, 1])"", \*""(Subscript[Z, max])""}, RotateLabel -> False]\)

The above bifurcation diagram highlights the period-doubling cascade of this model and its chaotic behaviour for parameter values (Subscript[δ, 1]) ≲ 0.63. The phase portrait of the Volterra-Gause model with the following parameters: ξ = 0.866, ε = 1.428, (Subscript[δ, 1]) = 0.577, (Subscript[δ, 2]) = 0.376 may be obtained by classical numerical integration and provides a chaotic attractor in the snail shell shape presented below.

Gause, G.F. [1935] The struggle for existence, Williams and Wilkins, Baltimore.

Ginoux, J.M., Rossetto, B. & Jamet, J.L. [2005] “Chaos in a three-dimensional Volterra-Gause model of predator-prey type,” Int. J. Bifurcation and Chaos  5, Vol. 15, 1689-1708.

Glass, L. & Mackey, M.M. [1988] From clocks to chaos (Princeton University Press,Princeton,New Jersey,USA).

Hopf, E., [1942] “Abzweigung einer periodischen Lösung von einer stationären Lösung eines differential-systems,” Bereich. Math. Phys. Kl. Sächs. Acad. Wiss. Leipzig, 94, 1-22, et Bereich. Math. Phys. Kl. Sächs. Acad. Wiss. Leipzig, 95(1), 3-22.

Lorenz, E.N. [1963] “Deterministic non-periodic flows,” J. Atmos. Sc 20, 130-141.

Shil'nikov, L. [1963] “A case of existence of a countable number of periodic motions,” Soviet. Math. Dokl. 6, 163-166.

Ruelle, D, & Takens, F. [1971] “On the nature of turbulence,” Communications of Mathematical Physics, 20, 167-192.

Verhulst, P.F. [1838] “Notice sur la loi que suit la population dans son accroissement,” Corresp. Math. Phys., X, 113-121.

Volterra, V. [1926] “Variazioni e fluttuazioni del numero d’individui in specie animali conviventi,” Mem. Acad. Lincei  III, 6, 31-113.