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