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
28nf():= block(
29
30/* new variable names? */
31test : read ("do you want to enter new variable names (y/n)?"),
32if test = n then go(jump),
33n : read ("how many eqs"),
34for i thru n do (x[i] : read ("symbol for old x[",i,"]")),
35for i thru n do( y[i] : read ("symbol for new x[",i,"]")),
36for i thru n do depends([x[i],y[i]],t),
37kill(flag), /* flag used in gen */
38
39jump,
40
41/* new d.e.'s? */
42print ("do you want to enter new d.e.'s (y/n)?"),
43test:read(),
44if test = n then go(loop),
45for 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])),
48kill(rhs2),
49rhs2[i,j] := rhs[i],
50rhs3:genmatrix(rhs2,n,1),
51
52loop,
53
54/* near-identity transformation */
55print("input near-identity transformation
56(use prev for previous transformation)"),
57for 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 */
64trans : makelist (x[i]=y[i]+tr[i],i,1,n),
65m : read("enter truncation order (highest order terms to be kept)"),
66
67
68/* transform the d.e.'s */
69temp2 :ev(rhs3,trans),
70
71/* solve for the transformed derivatives */
72kill(jacob),
73jacob[i,j]:=diff(tr[i],y[j]),
74jacob2:genmatrix(jacob,n,n),
75temp3:sum((-1)^i*jacob2^^i,i,0,m-1).temp2,
76
77/* taylor expand the resulting eqs */
78newrhs : taylor(temp3,makelist(y[i],i,1,n),0,m),
79newdes:makelist(diff(y[i],t)=newrhs[i,1],i,1,n),
80for i:1 thru n do
81          print(part(newdes,i)),
82
83/* enter another transformation? */
84branch:read("do you want to enter another transformation (y/n)"),
85if branch = y then go(loop),
86newdes)$
87
88
89gen(nn):=(
90if 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),
99ev(tempgen3[nn],rowdummy=row)) $
100
101decompose():=(
102kill(c),
103if not numberp(flag[m]) then gen(m),
104temp8:subst("[","+",tempgen3[m]),
105terms:ev(temp8,a[dummy,sub]:=1),
106coeffs:ev(temp8,a[dummy,sub]:=c[dummy,sub],makelist(y[i]=1,i,1,n)),
107for 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,
110j)))))$
111
112vars(nn):=(
113temp5:sum(ev(tempgen3[nn]),rowdummy,1,n),
114temp6:subst("[","+",temp5),
115temp7: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)))$
159