Study of the dynamic behavior of 
Three-Variable Autocatalator
Authors' Data: Housam Binous & Nasri Zakia
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Reference :
Bo Peng, Stephen K. Scott, and Kenneth Showalter, "Period Doubling and Chaos in a Three-Variable Autocatalator", The Journal of Physical Chemistry. Vol. 94, No. 13, 1990.

Off[General::"spell1"]

Period 1

κ=10;δ=2 10^-2;σ=5 10^-3;

μ=0.1;

sol1=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] == 0.01,β[0]== 0,γ[0]== 0},{α[t],β[t],γ[t]},{t,0,25},MaxSteps->1000000,AccuracyGoal-> 10]

Plot[Log[10, β[t]/.sol1],{t,20.1,20.5},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

Power Spectrum

Plot[Log[10, β[t]/.sol1],{t,20.1,20.2},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol1]== 0,{t,20.2}]

Plot[Log[10, β[t]/.sol1],{t,24,24.1},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol1]== 0,{t,24.08}]

tbl=Transpose[Table[Evaluate[β[t]/.sol1],{t,20.1896,24.0819,0.01}]];

fdata=Abs[Fourier[tbl[[1]]]];

ListPlot[fdata,Joined->True,PlotStyle-> RGBColor[0,0,1],PlotRange->{{0,200},{0,300}}]

Prominent Peaks

For[i=2,i<200,{i++,If[fdata[[i]]≥50,Print[i," ",fdata[[i]]]]}]

Period 2

μ=0.14;

sol2=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] == 0.01,β[0]== 0,γ[0]== 0},{α[t],β[t],γ[t]},{t,0,25},MaxSteps->1000000,AccuracyGoal-> 10]

Plot[Log[10, β[t]/.sol2],{t,20.15,21},PlotPoints-> 1000,PlotStyle-> RGBColor[1,0,0]]

Power Spectrum

Plot[Log[10, β[t]/.sol2],{t,20.1,20.2},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol2]== 0,{t,20.16}]

Plot[Log[10, β[t]/.sol2],{t,24,24.1},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol2]== 0,{t,24.08}]

tbl=Transpose[Table[Evaluate[β[t]/.sol2],{t,20.1601,24.0715,0.01}]];

fdata=Abs[Fourier[tbl[[1]]]];

ListPlot[fdata,Joined->True,PlotStyle-> RGBColor[0,0,1],PlotRange->{{0,200},{0,300}}]

Prominent Peaks

For[i=2,i<200,{i++,If[fdata[[i]]≥33,Print[i," ",fdata[[i]]]]}]

Period 4

μ=0.15;

sol3=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] == 0.01,β[0]== 0,γ[0]== 0},{α[t],β[t],γ[t]},{t,0,25},MaxSteps->1000000,AccuracyGoal-> 10]

Plot[Log[10, β[t]/.sol3],{t,20,20.8},PlotPoints-> 1000,PlotStyle-> RGBColor[1,0,0]]

Power Spectrum

Plot[Log[10, β[t]/.sol3],{t,20.1,20.2},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol3]== 0,{t,20.18}]

Plot[Log[10, β[t]/.sol3],{t,24.6,24.8},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol3]== 0,{t,24.71}]

tbl=Transpose[Table[Evaluate[β[t]/.sol3],{t,20.1837,24.7109,0.001}]];

fdata=Abs[Fourier[tbl[[1]]]];

ListPlot[fdata,Joined->True,PlotStyle-> RGBColor[0,0,1],PlotRange-> {{0,300},{0,400}}]

Prominent Peaks

For[i=2,i<200,{i++,If[fdata[[i]]≥50,Print[i," ",fdata[[i]]]]}]

Period 8

μ=0.151;

sol4=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] == 0.01,β[0]== 0,γ[0]== 0},{α[t],β[t],γ[t]},{t,0,25},MaxSteps->1000000,AccuracyGoal-> 10]

Plot[Log[10, β[t]/.sol4],{t,20,20.8},PlotPoints-> 1000,PlotStyle-> RGBColor[1,0,0]]

Power Spectrum

Plot[Log[10, β[t]/.sol4],{t,20.1,20.3},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol4]== 0,{t,20.2}]

Plot[Log[10, β[t]/.sol4],{t,24.1,24.3},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol4]== 0,{t,24.2}]

tbl=Transpose[Table[Evaluate[β[t]/.sol4],{t,20.1982,24.2033,0.001}]];

fdata=Abs[Fourier[tbl[[1]]]];

ListPlot[fdata,Joined->True,PlotStyle-> RGBColor[0,0,1],PlotRange->{{0,300},{0,300}}]

Chaos

μ=0.153;

sol5=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] == 0.01,β[0]== 0,γ[0]== 0},{α[t],β[t],γ[t]},{t,0,28},MaxSteps->1000000,AccuracyGoal-> 10]

Plot[Log[10, β[t]/.sol5],{t,20,21.5},PlotPoints-> 1000,PlotStyle-> RGBColor[1,0,0]]

Power Spectrum

Plot[Log[10, β[t]/.sol5],{t,20.2,20.4},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol5]== 0,{t,20.25,20.30}]

Plot[Log[10, β[t]/.sol5],{t,20.6,20.8},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol5]== 0,{t,20.6,20.62}]

tbl=Transpose[Table[Evaluate[β[t]/.sol5],{t,20.2618,20.6256,0.001}]];

fdata=Abs[Fourier[tbl[[1]]]];

ListPlot[fdata,Joined->True,PlotStyle-> RGBColor[0,0,1],PlotRange->{{0,200},{0,350}}]

Period 5

μ=0.155;

sol6=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] == 0.01,β[0]== 0,γ[0]== 0},{α[t],β[t],γ[t]},{t,0,25},MaxSteps->1000000,AccuracyGoal-> 10]

Plot[Log[10, β[t]/.sol6],{t,20.2,21.1},PlotPoints-> 1000,PlotStyle-> RGBColor[1,0,0]]

Power Spectrum

Plot[Log[10, β[t]/.sol6],{t,20.15,20.25},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol6]== 0,{t,20.22}]

Plot[Log[10, β[t]/.sol6],{t,24.6,24.8},PlotStyle-> RGBColor[1,0,0],PlotPoints-> 1000]

FindRoot[Log[10, β[t]/.sol6]== 0,{t,24.73}]

tbl=Transpose[Table[Evaluate[β[t]/.sol6],{t,20.2109,24.7262,0.001}]];

fdata=Abs[Fourier[tbl[[1]]]];

ListPlot[fdata,Joined->True,PlotStyle-> RGBColor[0,0,1],PlotRange->{{0,300},{0,250}}]

Prominent Peaks

For[i=2,i<200,{i++,If[fdata[[i]]≥50,Print[i," ",fdata[[i]]]]}]

Phase portrait of strange attractor in α-β-γ phase space 
for μ = 0.153.

μ=0.153;

sol7=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] == 0.01,β[0]== 0,γ[0]== 0},{α[t],β[t],γ[t]},{t,0,25},MaxSteps->1000000,AccuracyGoal-> 10]

ParametricPlot3D[Flatten[{Log[10, α[t]/.sol7],Log[10, β[t]/.sol7],Log[10, γ[t]/.sol7]}],{t,20,22},PlotPoints-> 10000,ViewPoint-> {-1,2,1}]

Phase portrait of period 5 limit cycle in α-β-γ phase space 
for μ = 0.155.

μ=0.155;

sol8=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] == 0.01,β[0]== 0,γ[0]== 0},{α[t],β[t],γ[t]},{t,0,25},MaxSteps->1000000,AccuracyGoal-> 10]

ParametricPlot3D[Flatten[{Log[10, α[t]/.sol8],Log[10, β[t]/.sol8],Log[10, γ[t]/.sol8]}],{t,20,22},PlotPoints-> 10000,ViewPoint-> {-1,2,1}]