by Jean-Marc Ginoux (ginoux@univ-tln.fr)

© Copyright Eric Javoy,computer scientist, 2005

http://ginoux.univ-tln.fr

This notebook contains an applet elaborated in order to provide the bifurcation diagram of a dynamical system. It just needs Mathematica (version 3, 4 or 5).

Bifurcation diagram with Mathematica

A. Dynamical system

Let's consider the following dynamical system :

((dx)/(dt))=(Subscript[x, 1])=f(x, y, z, λ)
((dy)/(dt))=(Subscript[y, 1])=g (x, y, z, μ)
((dz)/(dt))=(Subscript[z, 1])=h(x,y,z,ρ)

where λ, μ and ρ are parameters.

B. Bifurcation parameter

Choosing a bifurcation parameter between the set of paramters of the dynamical system. λ for example.

C. Details

This applet provides the bifurcation diagram of a three variables 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 dynamicla 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]) in order to remove the transient regime.

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

μ=;
ρ=;

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];
\ \ \ \ \ \ \ ]\)
]

\(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]\)

D. Example

This example provides the bifurcation diagram of a dynamical system elaborated by Ginoux et al. [2005]. It highlights the period doubling induced by the parameter (Subscript[δ, 1]). Other parametres are: ξ = 0.866, ε = 1.428, (Subscript[δ, 2]) = 0.376.

ξ ((d x)/(d t))=x (1 - x) - (x^((1/2))) y
   ((d y)/(d t))=-(Subscript[δ, 1]) y + (x^((1/2))) y - (y^((1/2))) 
            z     (12)
   ((d z)/(d t))=ε ((y^((1/2))) - (Subscript[δ, 2])) z

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

ξ = 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];
\ \ \ \ \ \ \ ]
]\)

Export["c:\TableauPlot.csv", TableauPlot]

\(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]\)