Hopf Bifurcation for the Brusselator

The differential equations for the “Brusselator”.  It’s an ODE in the plane with two parameters.

xdot = a-(b+1)*x+x^2*y;
ydot = b*x-x^2*y;

There is an equilibrium point at (a,b/a).  We change coordinates so this becomes the origin.

{xdot1,ydot1} = {xdot,ydot}/.{x->a+x,y->b/a+y}//Expand

Find the eigenvalues of the linearized system.

df = {{D[xdot1,x],D[xdot1,y]},{D[ydot1,x],D[ydot1,y]}}/.{x->0,y->0}//Simplify;
df//MatrixForm
Collect[CharacteristicPolynomial[df,z],z]
evals=Solve[%==0,z]//Simplify

There are some interesting bifurcations along the curve where b=1+a^2.  For simplicity we will set a=1 and consider a one-parameter family with parameter b.  The bifurcation will be at b=2.

{xdot2,ydot2} ={xdot1,ydot1}/.a->1
df2=df/.a->1;
df2//MatrixForm
Collect[CharacteristicPolynomial[df2,z],z]
evals=Solve[%==0,z]//Simplify;
evals/.b->2

The eigenvalues at b=2 are I, -I.  We can check that they eigenvalues pass through the imaginary axis with nonzero velocity as b varies.  In this case, the equilibrium goes from attractor to repeller as b increases.

realpart = (b-2)/2;alpha = D[realpart,b]/.b->2

Next we put the linearized matrix into real canonical form when a=1, b=2 using a linear change of variables.

Eigenvectors[df2/.b->2]
P={{-1,1},{2,0}};
{P//MatrixForm,Inverse[P].(df2/.b->2).P//MatrixForm}

Find the new nonlinear ODE after this coordinate change.

{xdot3,ydot3} =Inverse[P].({xdot2,ydot2}/.{b->2,x->-x+y,y->2*x})//Expand

Now plug the derivatives of this vectorfield into the magic formulas to get the coefficients a1, b1 of the normal form.

a1=((fxx fxy+fxy fyy-fxx gxx-gxx gxy+fyy gyy-gxy gyy)/(16 b))+(1/16) (fxxx+fxyy+gxxy+gyyy);
b1=(1/16) (-fxxy-fyyy+gxxx+gxyy)-(1/(48 b))(2 (fxx^2)+2 (fxy^2)+5 (fyy^2)+5 (gxx^2)+5 fxx (fyy-gxy)-fyy gxy+2 (gxy^2)+5 gxx gyy+2 (gyy^2)-fxy (gxx+5 gyy));
{f,g}={xdot3,ydot3};
sub = {fxx->D[f,{x,2}],fxy->D[f,{x,1},{y,1}],fyy->D[f,{y,2}],gxx->D[g,{x,2}],gxy->D[g,{x,1},{y,1}],gyy->D[g,{y,2}],
fxxx->D[f,{x,3}],fxxy->D[f,{x,2},{y,1}],fxyy->D[f,{x,1},{y,2}],fyyy->D[f,{y,3}],gxxx->D[g,{x,3}],gxxy->D[g,{x,2},{y,1}],gxyy->D[g,{x,1},{y,2}],gyyy->D[g,{y,3}]}/.{x->0,y->0}

Since we have alpha > 0 and a1 < 0 we have a supercritical Hopf bifurcation as b increases through b=2.

a10=a1/.sub/.b->2
b10 = b1/.sub/.b->2

Now we will have a look at the bifurcation using numerical integration.

BrusselatorOrbit[{x0_,y0_},{a_,b_},t0_]:={x[t],y[t]}/.NDSolve[{x'[t]==a-(b+1)*x[t]+x[t]^2*y[t],y'[t]==b*x[t]-x[t]^2*y[t],x[0]==x0,y[0]==y0},{x[t],y[t]},{t,0,t0}][[1]]

Here is the plot with a = 1 and b=1.8, b= 2.2 showing the birth of a periodic attractor.

MakeBrusselatorPlot[{x0_,y0_},{a_,b_},t0_]:= With[{orb=BrusselatorOrbit[{x0,y0},{a,b},t0]},
ParametricPlot[orb,{t,0,t0},PlotRange->{{0,4},{0,4}}]]

{MakeBrusselatorPlot[{1,0},{1.0,1.8},50],MakeBrusselatorPlot[{1,0},{1.0,2.2},50]}

Another way to visualize the bifurcation is to plot graphs of x(t), y(t) as functions of t.

MakeBrusselatorGraphs[{x0_,y0_},{a_,b_},t0_]:= With[{orb=BrusselatorOrbit[{x0,y0},{a,b},t0]},
Plot[orb,{t,0,t0},PlotRange->{{0,t0},{0,4}}]]

{MakeBrusselatorGraphs[{1,0},{1.0,1.8},200],MakeBrusselatorGraphs[{1,0},{1.0,2.2},200]}

We can use the Manipulate command to produce an interactive orbit calculator.

MakeBrusselatorPlot[{x0_,y0_},{a_,b_},t0_]:= With[{orb=BrusselatorOrbit[{x0,y0},{a,b},t0]},
ParametricPlot[orb,{t,0,t0},PlotRange->{{0,4},{0,4}}]]

Manipulate[MakeBrusselatorPlot[p,{1.0,b},t0],{{p,{1,0}},Locator},{{t0,0.1},0,50},{{b,1},1,3}]

MakeBrusselatorGraphs[{x0_,y0_},{a_,b_},t0_]:= With[{orb=BrusselatorOrbit[{x0,y0},{a,b},t0]},
Plot[orb,{t,0,t0},PlotRange->{{0,t0},{0,4}}]]

Manipulate[MakeBrusselatorGraphs[{1,0},{1.0,b},200],{{b,1},1,3}]