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