Computation of the Bifurcation Diagram
for the Three-Variable Autocatalator​
Authors: Housam Binous and Ahmed Bellagi

Consider a chemical reaction such as the shift conversion:
CO+
H
2
O
CO
2
+
H
2
. To study its chemical kinetics, it is written in abstract form as
X+Y→W+Z
. The rate of this reaction is given by
r=
dz
dt
=kxy
, where
k
is the rate constant and
x
,
y
, and
z
are the concentrations of the species
X
,
Y
, and
Z
.​Hypothetically, the autocatalator's steps are the following:​​
reactionrate​​P→A
k
0
p
0
​​P+C→A+C
k
c
p
0
c​​A→B
k
u
a​​A+2B→3B
k
1
a
2
b
​​B→C
k
2
b​​C→D
k
3
c
​​Here
P
is a precursor with constant concentration,
D
is the final product,
A
,
B
, and
C
are intermediate chemical species,
k
0
,
k
c
,
k
u
,
k
1
,
k
2
, and
k
3
are rate constants for the individual reactions, and
p
0
,
d
,
a
,
b
, and
c
are the concentrations of the hypothetical chemical species
P
,
D
,
A
,
B
, and
C
. ​This hypothetical reaction scheme is a model for a three-variable autocatalator. The autocatalytic reaction occurs in the step
A+2B→3B
, with
B
catalyzing its own formation. This step introduces the nonlinear term
k
1
a
2
b
in the governing equations, necessary in order to obtain sufficient complexity for chaos to occur. ​The rate equations for the three intermediate species can be written in the form:​​
da
dt
=
k
0
p
0
+
k
c
p
0
c-
k
u
a-
k
1
a
2
b
,​
db
dt
=
k
u
a+
k
1
a
2
b
-
k
2
b
,​
dc
dt
=
k
2
b-
k
3
c
.​The corresponding dimensionless rate equations are:​​
dα
dτ
=κ+μγ-α
2
β
-α
,​
σ
dβ
dτ
=α
2
β
+α-β
,​
δ
dγ
dτ
=β-γ
,​where
α
,
β
, and
γ
are the dimensionless concentrations of
A
,
B
, and
C
, and the four parameters
κ
,
μ
,
δ
, and
σ
depend on the rates of the individual reactions and the concentration of the precursor.​We compute the bifurcation diagram for this nonlinear chemical system. This bifurcation diagram, a remerging Feigenbaum tree, was first given in Ref. 1.​References:​[1] B. Peng, S. K. Scott, and K. Showalter, "Period Doubling and Chaos in a Three-Variable Autocatalator," The Journal of Physical Chemistry, 94(13), 1990 pp. 5243–5246.​[2] S. Lynch, Dynamical Systems with Applications using Mathematica, Boston: Birkhäuser, 2007.​
Bifurcation Diagram - Computation takes over one hour.
In[]:=
sol[μ_]:=NDSolve[{​​D[α[t],{t,1}]κ+μγ[t]-α[t]β[t]^2-α[t],​​σD[β[t],{t,1}]α[t]β[t]^2+α[t]-β[t],​​δD[γ[t],{t,1}]β[t]-γ[t],​​α[0]1,β[0]1,γ[0]1},{α[t],β[t],γ[t]},{t,0,20}]
In[]:=
κ=10;δ=210^-2;σ=510^-3;
In[]:=
iterate=Compile[{μ},{f=sol[μ];Map[{μ,#}&,FindMinValue[Log[10,γ[t]/.f],{t,#,#+0.01}]&/@Table[i,{i,4,10.0,0.01}]]}]
Out[]=
CompiledFunction
Argument count: 1
Argument types:
{_Real}

In[]:=
LaunchKernels[]
Out[]=
KernelObject
Local kernel
KernelID: 1
,KernelObject
Local kernel
KernelID: 2
,KernelObject
Local kernel
KernelID: 3
,KernelObject
Local kernel
KernelID: 4

In[]:=
pts=Parallelize[Table[iterate[μ],{μ,0.10,0.18,0.0001}]];
In[]:=
ListPlot[Flatten[pts,1],PlotStyleTable[{PointSize[0.003],Darker@Blue},{i,1,Length[pts]}],FrameTrue,FrameLabel{Style["μ",Italic],Style["Log[10,γ]",Italic]},ImageSize{400,400},AspectRatio1,GridLines->Automatic]
Out[]=