share/contrib/rand/nf.mac

   1 /* Filename nf.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 */
  16 /*program number 5: nf(), normal form transformations. see page 62 in
  17   "perturbation methods, bifurcation theory and computer algebra". */
  18
  19
  20 /* this file contains nf(), a normal form utility function. it also contains
  21    these additional functions:
  22    gen(n) will generate a homogeneous order n transformation.
  23    decompose() isolates the coefficients of the new equations.
  24    vars(n) generates a list of unknown coefficients of degree n.
  25    hopfk(), for k=2,3,4,5,6,7 solves for the coefficients of a system of
  26    2 de's so as to put the eqs in hopf normal form */
  27
  28 nf():= block(
  29
  30 /* new variable names? */
  31 test : read ("do you want to enter new variable names (y/n)?"),
  32 if test = n then go(jump),
  33 n : read ("how many eqs"),
  34 for i thru n do (x[i] : read ("symbol for old x[",i,"]")),
  35 for i thru n do( y[i] : read ("symbol for new x[",i,"]")),
  36 for i thru n do depends([x[i],y[i]],t),
  37 kill(flag), /* flag used in gen */
  38
  39 jump,
  40
  41 /* new d.e.'s? */
  42 print ("do you want to enter new d.e.'s (y/n)?"),
  43 test:read(),
  44 if test = n then go(loop),
  45 for i thru n do
  46           (rhs[i]:read("enter rhs of eq. no.",i,",d",x[i],"/dt ="),
  47            print("d",x[i],"/dt =",rhs[i])),
  48 kill(rhs2),
  49 rhs2[i,j] := rhs[i],
  50 rhs3:genmatrix(rhs2,n,1),
  51
  52 loop,
  53
  54 /* near-identity transformation */
  55 print("input near-identity transformation
  56 (use prev for previous transformation)"),
  57 for i thru n do
  58           (row:i,
  59            prev :tr[i],
  60            tr[i] :read (x[i],"=",y[i],"+ ?"),
  61            print (x[i],"=",y[i]+tr[i])),
  62
  63 /* input truncation order */
  64 trans : makelist (x[i]=y[i]+tr[i],i,1,n),
  65 m : read("enter truncation order (highest order terms to be kept)"),
  66
  67
  68 /* transform the d.e.'s */
  69 temp2 :ev(rhs3,trans),
  70
  71 /* solve for the transformed derivatives */
  72 kill(jacob),
  73 jacob[i,j]:=diff(tr[i],y[j]),
  74 jacob2:genmatrix(jacob,n,n),
  75 temp3:sum((-1)^i*jacob2^^i,i,0,m-1).temp2,
  76
  77 /* taylor expand the resulting eqs */
  78 newrhs : taylor(temp3,makelist(y[i],i,1,n),0,m),
  79 newdes:makelist(diff(y[i],t)=newrhs[i,1],i,1,n),
  80 for i:1 thru n do
  81           print(part(newdes,i)),
  82
  83 /* enter another transformation? */
  84 branch:read("do you want to enter another transformation (y/n)"),
  85 if branch = y then go(loop),
  86 newdes)$
  87
  88
  89 gen(nn):=(
  90 if not numberp(flag[nn]) then (
  91                sub:makelist(k[i],i,1,n),
  92                var:product(y[i]^k[i],i,1,n),
  93                tempgen:a[rowdummy,sub]*var,
  94                for i:1 thru n do
  95                        tempgen:sum(ev(tempgen,k[i]=j),j,0,nn),
  96                tempgen2:last(taylor(tempgen,makelist(y[i],i,1,n),0,nn)),
  97                tempgen3[nn]:expand(tempgen2),
  98                flag[nn] : 1),
  99 ev(tempgen3[nn],rowdummy=row)) $
 100
 101 decompose():=(
 102 kill(c),
 103 if not numberp(flag[m]) then gen(m),
 104 temp8:subst("[","+",tempgen3[m]),
 105 terms:ev(temp8,a[dummy,sub]:=1),
 106 coeffs:ev(temp8,a[dummy,sub]:=c[dummy,sub],makelist(y[i]=1,i,1,n)),
 107 for i:1 thru n do(
 108      for j:1 thru length(terms) do(
 109           ev(part(coeffs,j),rowdummy=i)::ratcoef(expand(newrhs[i,1]),part(terms,
 110 j)))))$
 111
 112 vars(nn):=(
 113 temp5:sum(ev(tempgen3[nn]),rowdummy,1,n),
 114 temp6:subst("[","+",temp5),
 115 temp7:ev(temp6,makelist(y[i]=1,i,1,n)))$
 116
 117
 118  hopf2():=(decompose(),
 119            solve([c[1,[2,0]],c[1,[1,1]],c[1,[0,2]],c[2,[2,0]],
 120                   c[2,[1,1]],c[2,[0,2]]],
 121                  vars(2)))$
 122
 123  hopf3():=(decompose(),
 124            solve([c[1,[3,0]]=c[1,[1,2]],c[1,[3,0]]=c[2,[2,1]],
 125                   c[1,[3,0]]=c[2,[0,3]],c[1,[0,3]]=c[1,[2,1]],
 126                   c[1,[0,3]]=-c[2,[3,0]],c[1,[0,3]]=-c[2,[1,2]]],
 127                  vars(3)))$
 128
 129  hopf4():=(decompose(),
 130            solve([c[1,[4,0]],c[1,[3,1]],c[1,[2,2]],c[1,[1,3]],
 131                   c[1,[0,4]],c[2,[4,0]],c[2,[3,1]],c[2,[2,2]],
 132                   c[2,[1,3]],c[2,[0,4]]],
 133                  vars(4)))$
 134
 135  hopf5():=(decompose(),
 136            solve([c[1,[5,0]]=c[1,[3,2]]/2,c[1,[5,0]]=c[1,[1,4]],
 137                   c[1,[5,0]]=c[2,[4,1]],c[1,[5,0]]=c[2,[2,3]]/2,
 138                   c[1,[5,0]]=c[2,[0,5]],c[2,[5,0]]=c[2,[3,2]]/2,
 139                   c[2,[5,0]]=c[2,[1,4]],c[2,[5,0]]=-c[1,[4,1]],
 140                   c[2,[5,0]]=-c[1,[2,3]]/2,c[2,[5,0]]=-c[1,[0,5]]],
 141                  vars(5)))$
 142
 143  hopf6():=(decompose(),
 144            solve([c[1,[6,0]],c[1,[5,1]],c[1,[4,2]],c[1,[3,3]],
 145                   c[1,[2,4]],c[1,[1,5]],c[1,[0,6]],c[2,[6,0]],
 146                   c[2,[5,1]],c[2,[4,2]],c[2,[3,3]],c[2,[2,4]],
 147                   c[2,[1,5]],c[2,[0,6]]],
 148                  vars(6)))$
 149
 150  hopf7():=(decompose(),
 151            solve([c[1,[7,0]]=c[1,[5,2]]/3,c[1,[7,0]]=c[1,[3,4]]/3,
 152                   c[1,[7,0]]=c[1,[1,6]],c[1,[7,0]]=c[2,[6,1]],
 153                   c[1,[7,0]]=c[2,[4,3]]/3,c[1,[7,0]]=c[2,[2,5]]/3,
 154                   c[1,[7,0]]=c[2,[0,7]],c[2,[7,0]]=c[2,[5,2]]/3,
 155                   c[2,[7,0]]=c[2,[3,4]]/3,c[2,[7,0]]=c[2,[1,6]],
 156                   c[2,[7,0]]=-c[1,[6,1]],c[2,[7,0]]=-c[1,[4,3]]/3,
 157                   c[2,[7,0]]=-c[1,[2,5]]/3,c[2,[7,0]]=-c[1,[0,7]]],
 158                  vars(7)))$