1/* Filename lie.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 9: lie(), lie transformations for hamiltonian systems. 17 see page 144 in "perturbation methods, bifurcation theory and 18 computer algebra". */ 19 20 21lie():=block( 22 23/* input problem ? */ 24choice1:read("do you want to input a new problem (y/n) ?"), 25 26if choice1=n then go(jump1), 27 28/* input problem */ 29n:read("enter number of degrees of freedom"), 30for ii:1 thru n do ( 31 q[ii]:read("enter symbol for q[",ii,"]"), 32 p[ii]:read("enter symbol for p[",ii,"]")), 33kill(w), 34print("the hamiltonian depends on the q's, p's, t and e (small parameter)"), 35print("the e=0 problem must be of the form:"), 36print("h =",sum((p[ii]^2+w[ii]^2*q[ii]^2)/2,ii,1,n)), 37horiginal:read("enter the hamiltonian"), 38print("h =",horiginal), 39 40/* transform to action-angle variables */ 41/* find the w[ii]'s */ 42h0:ev(horiginal,e=0), 43for ii:1 thru n do 44 w[ii]:sqrt(diff(h0,q[ii],2)), 45print("the action-angle variables are j's for action, phi's for angle"), 46for ii:1 thru n do 47 tr[ii]:[q[ii]=sqrt(2*j[ii]/w[ii])*sin(phi[ii]), 48 p[ii]=sqrt(2*j[ii]*w[ii])*cos(phi[ii])], 49tran:makelist(tr[ii],ii,1,n), 50h:ev(horiginal,tran,assume_pos:true,infeval), 51h:trigsimp(h), 52h:expand(trigreduce(expand(h))), 53print("h =",h), 54 55jump1, 56 57/* input truncation order */ 58ntrunc:read("enter highest order terms in e to be kept"), 59for ii:0 thru ntrunc do 60 h[ii]:ratcoef(h,e,ii), 61 62/* lie transforms */ 63/* near identity transformation from (j,phi)'s to (i,psi)'s */ 64 65/* update variables */ 66for ii:1 thru n do( 67 p[ii]:i[ii], 68 q[ii]:psi[ii]), 69 70/* replace j and phi by i and psi in h's */ 71for ii:0 thru ntrunc do 72 h[ii]:ev(h[ii],makelist(j[iii]=i[iii],iii,1,n), 73 makelist(phi[iii]=psi[iii],iii,1,n)), 74 75k[0]:h[0], 76 77/* declare wgen[i] to be a fn of t, q's and p's */ 78kill(wgen), 79depends(wgen,[t]), 80for ii:1 thru n do 81 depends(wgen,[q[ii],p[ii]]), 82 83/* e=0 problem is of form sum(w[ii]*i[ii]) */ 84/* choose wgen[ii] to kill as much as possible in eq(ii) */ 85/* equate k[ii] to unremovable terms */ 86 87/* define pattern matching rules to isolate args of trig terms */ 88matchdeclare([dummy1,dummy2],true), 89defrule(cosarg,dummy1*cos(dummy2),dummy2), 90defrule(sinarg,dummy1*sin(dummy2),dummy2), 91 92for ii:1 thru ntrunc do( 93 94 eqn[ii]:expand(trigreduce(expand(eq(ii)))), 95 temp:expand(ev(eqn[ii],wgen[ii]=0)), 96/* change sum to a list */ 97 temp1:args(temp), 98/* remove constant terms */ 99 temp2:map(trigidentify,temp1), 100/* isolate arguments of trig terms */ 101 arg1:apply1(temp2,cosarg,sinarg), 102/* remove duplicate arguments */ 103 arg2:setify(arg1), 104/* remove resonant arguments */ 105 arg3:sublist(arg2,notresp), 106 107/* choose wgen to eliminate nonresonant terms */ 108 leng:length(arg3), 109 wgentemp1:0, 110 for jj:1 thru leng do( 111 wgentemp2:aaa*cos(part(arg3,jj))+bbb*sin(part(arg3,jj)), 112 temp4:ev(eqn[ii],wgen[ii]=wgentemp2,diff), 113 temp5:solve([ratcoef(temp4,cos(part(arg3,jj))), 114 ratcoef(temp4,sin(part(arg3,jj)))],[aaa,bbb]), 115 wgentemp1:wgentemp1+ev(wgentemp2,temp5)), 116 wgen[ii]:wgentemp1, 117 print("wgen[",ii,"] ="), 118 print(wgen[ii]), 119 k[ii]:expand(ev(eqn[ii],diff)), 120 k[ii]:expand(ratsimp(k[ii])), 121 print("the transformed hamiltonian k[",ii,"] ="), 122 print(k[ii])), 123 124kamiltonian:sum(k[ii]*e^ii,ii,0,ntrunc), 125print ("the transformed hamiltonian is "), 126print ("k =",kamiltonian), 127 128choice2:read("do you want to see the near identity transformation (y/n) ?"), 129if choice2=n then go(end), 130 131/* the near identity transformation */ 132for ii:1 thru n do( 133 jtrans[ii]:sum(s(iii,p[ii])*e^iii,iii,0,ntrunc), 134 phitrans[ii]:sum(s(iii,q[ii])*e^iii,iii,0,ntrunc)), 135for ii:1 thru n do( 136 print(j[ii],"=",jtrans[ii]), 137 print(phi[ii],"=",phitrans[ii])), 138 139end, 140kamiltonian)$ 141 142 143poisson(f,g):= 144 sum(diff(f,q[ii])*diff(g,p[ii])-diff(f,p[ii])*diff(g,q[ii]),ii,1,n)$ 145 146l(ii,f):=poisson(wgen[ii],f)$ 147 148s(ii,f):=(if ii=0 then f else sum(l(ii-m,s(m,f)),m,0,ii-1)/ii)$ 149 150eq(ii):=(h[ii]+(diff(wgen[ii],t)+poisson(wgen[ii],h[0]))/ii 151 +sum(l(ii-m,k[m])+m*s(ii-m,h[m]),m,1,ii-1)/ii)$ 152 153lzap(any):=diff(any,t)+poisson(any,h[0])$ 154 155trigidentify(exp):=if freeof(sin,exp) and freeof(cos,exp) then 0 else exp$ 156 157notresp(exp):=not is(lzap(exp) = 0)$ 158 159setify(list):=( 160 set:[list[1]], 161 for i:2 thru length(list) do( 162 if not member(list[i],set) then set:cons(list[i],set)), 163 set)$ 164