1/* Filename cm.mac 2 3 *************************************************************** 4 * * 5 * <package name> * 6 * <functionality description> * 7 * * 8 * from: Perturbation Methods, Bifurcation * 9 * Theory and Computer Algebra. * 10 * by Rand & Armbruster (Springer 1987) * 11 * Programmed by Richard Rand * 12 * These files are released to the public domain * 13 * * 14 *************************************************************** 15*/ /* program number 3: cm(), center manifold reduction for ordinary 16 differential equations. see page 32 in "perturbation methods, 17 bifurcation theory and computer algebra". */ 18 19 20 21cm():=( 22 23/* input problem */ 24n:read("enter no. of eqs."), 25m:read("enter dimension of center manifold"), 26print("the d.e.'s must be arranged so that the first",m,"eqs."), 27print("represent the center manifold. i.e. all associated"), 28print("eigenvalues are zero or have zero real parts."), 29for i:1 thru n do 30 x[i]:read("enter symbol for variable no.",i), 31l:read("enter order of truncation"), 32for i:1 thru n do ( 33 print("enter rhs of eq.",i), 34 print("d",x[i],"/dt ="), 35 g[i]:read()), 36/* set up d.e.'s */ 37for i:1 thru n do 38 depends(x[i],t), 39for i:1 thru n do 40 (eq[i]:diff(x[i],t)=g[i], 41 print(eq[i])), 42 43/* form power series */ 44sub:makelist(k[i],i,1,m), 45var:product(x[i]^k[i],i,1,m), 46unk:[], 47for p:m+1 thru n do( 48 temp:a[p,sub]*var, 49 for i:1 thru m do 50 temp:sum(ev(temp,k[i]=j),j,0,l), 51 temp2:taylor(temp,makelist(x[i],i,1,m),0,l), 52 /* remove constant and linear terms */ 53 temp3:temp2-part(temp2,1)-part(temp2,2), 54 soln[p]:expand(temp3), 55 /* prepare list of unknowns */ 56 setxto1:makelist(x[i]=1,i,1,m), 57 /* turn sum into a list */ 58 unkn[p]:subst("[","+",ev(temp3,setxto1)), 59 unk:append(unk,unkn[p])), 60sol:makelist(x[p]=soln[p],p,m+1,n), 61 62/* substitute into d.e.'s */ 63cmde:makelist(eq[i],i,1,m), 64rest:makelist(lhs(eq[i])-rhs(eq[i]),i,m+1,n), 65temp4:ev(rest,sol,diff), 66temp5:ev(temp4,cmde,diff), 67temp6:ev(temp5,sol), 68temp7:taylor(temp6,makelist(x[i],i,1,m),0,l), 69 70/* collect terms */ 71counter:1, 72/* make list of terms */ 73terms:subst("[","+",soln[n]), 74terms:ev(terms,a[dummy,sub]:=1), 75for i:1 thru n-m do ( 76 exp[i]:expand(part(temp7,i)), 77 for j:1 thru length(terms) do( 78 cond[counter]:ratcoef(exp[i],part(terms,j)), 79 counter:counter+1)), 80conds:makelist(cond[i],i,1,counter-1), 81conds:ev(conds,makelist(x[i]=0,i,1,m)), 82 83/* solve for center manifold */ 84acoeffs:solve(conds,unk), 85centermanifold:ev(sol,acoeffs), 86print("center manifold:"), 87print(centermanifold), 88 89/* get flow on cm */ 90cmde2:ev(cmde,centermanifold), 91print("flow on the c.m.:"), 92print(cmde2))$ 93