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