share/contrib/rand/cm.mac

   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
  21 cm():=(
  22
  23 /* input problem */
  24 n:read("enter no. of eqs."),
  25 m:read("enter dimension of center manifold"),
  26 print("the d.e.'s must be arranged so that the first",m,"eqs."),
  27 print("represent the center manifold.  i.e. all associated"),
  28 print("eigenvalues are zero or have zero real parts."),
  29 for i:1 thru n do
  30    x[i]:read("enter symbol for variable no.",i),
  31 l:read("enter order of truncation"),
  32 for 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 */
  37 for i:1 thru n do
  38    depends(x[i],t),
  39 for i:1 thru n do
  40    (eq[i]:diff(x[i],t)=g[i],
  41     print(eq[i])),
  42
  43 /* form power series */
  44 sub:makelist(k[i],i,1,m),
  45 var:product(x[i]^k[i],i,1,m),
  46 unk:[],
  47 for 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])),
  60 sol:makelist(x[p]=soln[p],p,m+1,n),
  61
  62 /* substitute into d.e.'s */
  63 cmde:makelist(eq[i],i,1,m),
  64 rest:makelist(lhs(eq[i])-rhs(eq[i]),i,m+1,n),
  65 temp4:ev(rest,sol,diff),
  66 temp5:ev(temp4,cmde,diff),
  67 temp6:ev(temp5,sol),
  68 temp7:taylor(temp6,makelist(x[i],i,1,m),0,l),
  69
  70 /* collect terms */
  71 counter:1,
  72 /* make list of terms */
  73 terms:subst("[","+",soln[n]),
  74 terms:ev(terms,a[dummy,sub]:=1),
  75 for 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)),
  80 conds:makelist(cond[i],i,1,counter-1),
  81 conds:ev(conds,makelist(x[i]=0,i,1,m)),
  82
  83 /* solve for center manifold */
  84 acoeffs:solve(conds,unk),
  85 centermanifold:ev(sol,acoeffs),
  86 print("center manifold:"),
  87 print(centermanifold),
  88
  89 /* get flow on cm */
  90 cmde2:ev(cmde,centermanifold),
  91 print("flow on the c.m.:"),
  92 print(cmde2))$