This program calculates the amplitude eqation for a Hopf bifurcation in a 2 component reaction diffusion equation.

> restart;

> mkoeff:=proc(exp,var1,max1,var2,max2)
global koeff2:
local mk_exp, ausdruck, ivar1, mk_exp1, ivar2:
mk_exp:=collect(exp,var1):
for ivar1 from 0 to max1 do
mk_exp1[ivar1]:=coeff(mk_exp,var1,ivar1):
mk_exp1[ivar1]:=collect(mk_exp1[ivar1],var2):
for ivar2 from 0 to max2 do
koeff2[ivar1,ivar2]:=coeff(mk_exp1[ivar1],var2,ivar2):
ausdruck:=var1^ivar1*var2^ivar2:
print(ausdruck,koeff2[ivar1,ivar2]):
od:
od:
end:

> abl:=proc(func,var)
if (var=t) then
diff(func,t)+epsilon*diff(func,t1)+epsilon**2*diff(func,t2)+
epsilon**3*diff(func,t3)+epsilon**4*diff(func,t4)
elif (var=x) then
diff(func,x)+epsilon*diff(func,x1)
elif (var=y) then
diff(func,y)+epsilon*diff(func,y1)
else
diff(func,var)
fi
end:

> alias(x=xf(t,t2),y=yf(t,t2),A1=Af(t2),B21=B21f(t2),B22=B22f(t2),C21=C21f(t2),C22=C22f(t2),A21=A21f(t2),A22=A22f(t2),A31=A31f(t2),A32=A32f(t2),B31=B31f(t2),B32=B32f(t2),C31=C31f(t2),C32=C32f(t2),D31=D31f(t2),D32=D32f(t2),A1s=Asf(t2),B21s=B21sf(t2),B22s=B22sf(t2),C21s=C2xsf(t2),C22s=C22sf(t2),A21s=A21sf(t2),A22s=A22sf(t2),A31s=A31sf(t2),A32s=A32sf(t2),B31s=B31sf(t2),B32s=B32sf(t2),C31s=C31sf(t2),C32s=C32sf(t2),D31s=D31sf(t2),D32s=D32sf(t2));

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> #A:=1; with A=1 there are no DC-contributions and the results become much simpler

> gl1:=-abl(x,t)+A-(B+1)*x+x^2*y;

> gl2:=-abl(y,t)+B*x-x^2*y;

determine the critical eigenvector

defxyl:=[x=A+x0*exp(I*omega*t),y=B/A+y0*exp(I*omega*t)];

> expand(subs(defxyl,gl1)):

> coeff(%,exp(I*omega*t));

> defy0:=y0=solve(%,y0);

> expand(subs(defxyl,gl2));

> subs(defy0,%):

> lin2:=factor(coeff(%,exp(I*omega*t)));

> defom:=omega=solve(evalc(Re(lin2)),omega)[1];

> defBc:=Bc=solve(subs(defom,lin2),B);

> subs([B=Bc],defy0);

> defyp:=yp=expand(subs([defBc,defom,x0=1],rhs(%)));

> defyps:=yps=subs(I=-I,rhs(defyp)); yp is the y-component of the right-eigenvector

> defB:=B=Bc+B2*epsilon^2;

> defypp:=ypp=factor(subs(defom,A^2/(I*omega+A^2))); ypp is the y-component of the left-eigenvector

> defypps:=ypps=factor(subs(defom,A^2/(-I*omega+A^2)));

expansion of x and y

> defx:=x=A+epsilon*(A1*exp(I*omega*t)+A1s*exp(-I*omega*t))+epsilon^2*((A21+A22)*exp(I*omega*t)+(A21s+A22s)*exp(-I*omega*t)+B21+B21s+(C21+C22)*exp(2*I*omega*t)+(C21s+C22s)*exp(-2*I*omega*t))+epsilon^3*((A31+A32)*exp(I*omega*t)+(A31s+A32s)*exp(-I*omega*t)+B31+(C31+C32)*exp(2*I*omega*t)+(C31s+C32s)*exp(-2*I*omega*t)+(D31+D32)*exp(3*I*omega*t)+(D31s+D32s)*exp(-3*I*omega*t));

> defy:=y=B/A+epsilon*(yp*A1*exp(I*omega*t)+yps*A1s*exp(-I*omega*t))+epsilon^2*((yp*A21+yps*A22)*exp(I*omega*t)+(yps*A21s+yp*A22s)*exp(-I*omega*t)+yp*B21+yps*B21s+(C21*yp+C22*yps)*exp(2*I*omega*t)+(yps*C21s+yp*C22s)*exp(-2*I*omega*t))+epsilon^3*((yp*A31+yps*A32)*exp(I*omega*t)+(yps*A31s+yp*A32s)*exp(-I*omega*t)+yp*B31+yps*B31s+(yp*C31+yps*C32)*exp(2*I*omega*t)+(yps*C31s+yp*C32s)*exp(-2*I*omega*t)+(yp*D31+yps*D32)*exp(3*I*omega*t)+(yps*D31s+yp*D32s)*exp(-3*I*omega*t));

>

> subs([defx,defy,defB],gl1):

> subs([defyp,defyps,defB],%):subs([defBc],%):

> gl1e:=collect(expand(convert(taylor(%,epsilon=0,4),polynom)),[epsilon,exp(I*omega*t)]):

>

>

> subs([defx,defy,defB],gl2):

> subs([defyp,defyps,defB],%): subs([defBc],%):

> gl2e:=collect(expand(convert(taylor(%,epsilon=0,4),polynom)),[epsilon,exp(I*omega*t)]):

check whether eigenvector correct

> collect(expand(coeff(gl1e,epsilon,1)),exp(I*omega*t)):

> expand(subs([defBc,defom],%));

> collect(expand(coeff(gl2e,epsilon,1)),exp(I*omega*t)):

> expand(subs([defBc,defom],%));

> linop1:=collect(coeff(expand(subs(defxyl,gl1)),exp(I*omega*t)),[x0,y0]);

> linop2:=collect(coeff(expand(subs(defxyl,gl2)),exp(I*omega*t)),[x0,y0]);

> collect(subs([defom,defB,defypp,epsilon=0],linop1+ypp*linop2),[x0,y0],factor);

> subs(defBc,%);

> gl1e2:=collect(expand(coeff(gl1e,epsilon,2)),exp(I*omega*t));

>

> gl1e20:=coeff(gl1e2,exp(I*omega*t),0);

> gl1e21:=coeff(gl1e2,exp(I*omega*t),1);

> gl1e22:=coeff(gl1e2,exp(I*omega*t),2);

> gl1e2m1:=coeff(gl1e2,exp(I*omega*t),-1);

> gl1e2m2:=coeff(gl1e2,exp(I*omega*t),-2);

> gl2e2:=collect(expand(coeff(gl2e,epsilon,2)),exp(I*omega*t));

> gl2e20:=coeff(gl2e2,exp(I*omega*t),0);

> gl2e21:=coeff(gl2e2,exp(I*omega*t),1);

> gl2e22:=coeff(gl2e2,exp(I*omega*t),2);gl2e2m1:=coeff(gl2e2,exp(I*omega*t),-1);

> gl2e2m2:=coeff(gl2e2,exp(I*omega*t),-2);

> defA22:=A22=subs(defom,solve(gl2e21,A22));

> defA22s:=A22s=subs(defom,solve(gl2e2m1,A22s));

> defB2:=solve({gl1e20,gl2e20},{B21,B21s});

> defC2:=solve({gl1e22,gl2e22},{C21,C22});

> defC2s:=solve({gl1e2m2,gl2e2m2},{C21s,C22s});

test whether equations at quadratic order solved correctly

> subs(defom,normal(subs([defA22,defA22s,defB2[1],defB2[2],defC2[1],defC2[2],defC2s[1],defC2s[2]],gl2e2)));

> subs(defom,normal(subs([defA22,defA22s,defB2[1],defB2[2],defC2[1],defC2[2],defC2s[1],defC2s[2]],gl1e2)));

>

> gl1e3:=collect(expand(coeff(gl1e,epsilon,3)),exp(I*omega*t));

> gl2e3:=collect(expand(coeff(gl2e,epsilon,3)),exp(I*omega*t));

> subs([defB2[1],defB2[2],defC2[1],defC2[2],defC2s[1],defC2s[2]],gl1e3):

> collect(expand(%),exp(I*omega*t)):

> gl1e31:=collect(coeff(%,exp(I*omega*t),1),[A1,A1s],distributed,factor);

>

> subs([defB2[1],defB2[2],defC2[1],defC2[2],defC2s[1],defC2s[2]],gl2e3):

> collect(expand(%),exp(I*omega*t)):

> gl2e31:=collect(coeff(%,exp(I*omega*t),1),[A1,A1s],distributed,factor);

>

project on left-zero-eigenvector to yield the solvability condition, i.e. the amplitude equation

> normal(gl1e31+ypp*gl2e31):

> subs(defypp,%):

> subs([defom,defBc],%):

> normal(evalc(collect(%,[A1,A1s],distributed,factor))):

> cgl:=diff(A1,t2)=solve(%,diff(A1,t2)):

> collect(%,[A1,A1s],factor);

Test whether equations are actually satisfied to this order by the solutions determined

> gl1e31e:=subs(cgl,gl1e31);

> gl2e31e:=subs(cgl,gl2e31);

> defA32:=A32=collect(subs(defom,solve(gl2e31e,A32)),[A1,A1s],distributed,factor);

> normal(subs(defom,subs(defA32,gl1e31e)));

> collect(subs([defA22,defA22s,defB2[1],defB2[2],defC2[1],defC2[2],defC2s[1],defC2s[2],cgl,defA32,defBc],gl1e),[epsilon,exp(I*omega*t)],factor):subs(defom,%);test:=coeff(%,exp(I*A*t));

> collect(subs([defA22,defA22s,defB2[1],defB2[2],defC2[1],defC2[2],defC2s[1],defC2s[2],cgl,defA32,defBc],gl2e),[epsilon,exp(I*omega*t)],factor):subs(defom,%);test:=coeff(%,exp(I*A*t));

>