/* PROGRAM NUMBER 9: LIE(), LIE TRANSFORMATIONS FOR HAMILTONIAN SYSTEMS.
   SEE PAGE 144 IN "PERTURBATION METHODS, BIFURCATION THEORY AND
   COMPUTER ALGEBRA". */


LIE():=BLOCK(
 
/* INPUT PROBLEM ? */
CHOICE1:READ("DO YOU WANT TO INPUT A NEW PROBLEM (Y/N) ?"),
 
IF CHOICE1=N THEN GO(JUMP1),
 
/* INPUT PROBLEM */
N:READ("ENTER NUMBER OF DEGREES OF FREEDOM"),
FOR II:1 THRU N DO (
   Q[II]:READ("ENTER SYMBOL FOR Q[",II,"]"),
   P[II]:READ("ENTER SYMBOL FOR P[",II,"]")),
KILL(W),
PRINT("THE HAMILTONIAN DEPENDS ON THE Q'S, P'S, T AND E (SMALL PARAMETER)"),
PRINT("THE E=0 PROBLEM MUST BE OF THE FORM:"),
PRINT("H =",SUM((P[II]^2+W[II]^2*Q[II]^2)/2,II,1,N)),
HORIGINAL:READ("ENTER THE HAMILTONIAN"),
PRINT("H =",HORIGINAL),
 
/* TRANSFORM TO ACTION-ANGLE VARIABLES */
/* FIND THE W[II]'S */
H0:EV(HORIGINAL,E=0),
FOR II:1 THRU N DO
   W[II]:SQRT(DIFF(H0,Q[II],2)),
PRINT("THE ACTION-ANGLE VARIABLES ARE J'S FOR ACTION, PHI'S FOR ANGLE"), 
FOR II:1 THRU N DO
   TR[II]:[Q[II]=SQRT(2*J[II]/W[II])*SIN(PHI[II]),
          P[II]=SQRT(2*J[II]*W[II])*COS(PHI[II])],
TRAN:MAKELIST(TR[II],II,1,N),
H:EV(HORIGINAL,TRAN,ASSUME_POS:TRUE,INFEVAL),
H:TRIGSIMP(H),
H:EXPAND(TRIGREDUCE(EXPAND(H))),
PRINT("H =",H),
 
JUMP1,
 
/* INPUT TRUNCATION ORDER */
NTRUNC:READ("ENTER HIGHEST ORDER TERMS IN E TO BE KEPT"),
FOR II:0 THRU NTRUNC DO
   H[II]:RATCOEF(H,E,II),
 
/* LIE TRANSFORMS */
/* NEAR IDENTITY TRANSFORMATION FROM (J,PHI)'S TO (I,PSI)'S */
 
/* UPDATE VARIABLES */
FOR II:1 THRU N DO(
   P[II]:I[II],
   Q[II]:PSI[II]),
 
/* REPLACE J AND PHI BY I AND PSI IN H'S */
FOR II:0 THRU NTRUNC DO
   H[II]:EV(H[II],MAKELIST(J[III]=I[III],III,1,N),
                MAKELIST(PHI[III]=PSI[III],III,1,N)),
 
K[0]:H[0],
 
/* DECLARE WGEN[I] TO BE A FN OF T, Q'S AND P'S */
KILL(WGEN),
DEPENDS(WGEN,[T]),
FOR II:1 THRU N DO
  DEPENDS(WGEN,[Q[II],P[II]]),
 
/* E=0 PROBLEM IS OF FORM SUM(W[II]*I[II]) */
/* CHOOSE WGEN[II] TO KILL AS MUCH AS POSSIBLE IN EQ(II) */
/* EQUATE K[II] TO UNREMOVABLE TERMS */
 
/* DEFINE PATTERN MATCHING RULES TO ISOLATE ARGS OF TRIG TERMS */
MATCHDECLARE([DUMMY1,DUMMY2],TRUE),
DEFRULE(COSARG,DUMMY1*COS(DUMMY2),DUMMY2),
DEFRULE(SINARG,DUMMY1*SIN(DUMMY2),DUMMY2),
 
FOR II:1 THRU NTRUNC DO(
 
   EQN[II]:EXPAND(TRIGREDUCE(EXPAND(EQ(II)))),
   TEMP:EXPAND(EV(EQN[II],WGEN[II]=0)),
/* CHANGE SUM TO A LIST */
   TEMP1:ARGS(TEMP),
/* REMOVE CONSTANT TERMS */
   TEMP2:MAP(TRIGIDENTIFY,TEMP1),
/* ISOLATE ARGUMENTS OF TRIG TERMS */
   ARG1:APPLY1(TEMP2,COSARG,SINARG),
/* REMOVE DUPLICATE ARGUMENTS */
   ARG2:SETIFY(ARG1),
/* REMOVE RESONANT ARGUMENTS */
   ARG3:SUBLIST(ARG2,NOTRESP),
 
/* CHOOSE WGEN TO ELIMINATE NONRESONANT TERMS */
   LENG:LENGTH(ARG3),
   WGENTEMP1:0,
   FOR JJ:1 THRU LENG DO(
      WGENTEMP2:AAA*COS(PART(ARG3,JJ))+BBB*SIN(PART(ARG3,JJ)),
      TEMP4:EV(EQN[II],WGEN[II]=WGENTEMP2,DIFF),
      TEMP5:SOLVE([RATCOEF(TEMP4,COS(PART(ARG3,JJ))),
                   RATCOEF(TEMP4,SIN(PART(ARG3,JJ)))],[AAA,BBB]),
      WGENTEMP1:WGENTEMP1+EV(WGENTEMP2,TEMP5)),
   WGEN[II]:WGENTEMP1,
   PRINT("WGEN[",II,"] ="),
   PRINT(WGEN[II]),
   K[II]:EXPAND(EV(EQN[II],DIFF)),
   K[II]:EXPAND(RATSIMP(K[II])),
   PRINT("THE TRANSFORMED HAMILTONIAN K[",II,"] ="),
   PRINT(K[II])),   
 
KAMILTONIAN:SUM(K[II]*E^II,II,0,NTRUNC),
PRINT ("THE TRANSFORMED HAMILTONIAN IS "),
PRINT ("K =",KAMILTONIAN),
 
CHOICE2:READ("DO YOU WANT TO SEE THE NEAR IDENTITY TRANSFORMATION (Y/N) ?"),
IF CHOICE2=N THEN GO(END),
 
/* THE NEAR IDENTITY TRANSFORMATION */
FOR II:1 THRU N DO(
   JTRANS[II]:SUM(S(III,P[II])*E^III,III,0,NTRUNC),
   PHITRANS[II]:SUM(S(III,Q[II])*E^III,III,0,NTRUNC)),
FOR II:1 THRU N DO(
   PRINT(J[II],"=",JTRANS[II]),
   PRINT(PHI[II],"=",PHITRANS[II])),
 
END,
KAMILTONIAN)$
 
 
POISSON(F,G):=
          SUM(DIFF(F,Q[II])*DIFF(G,P[II])-DIFF(F,P[II])*DIFF(G,Q[II]),II,1,N)$
 
L(II,F):=POISSON(WGEN[II],F)$
 
S(II,F):=(IF II=0 THEN F ELSE SUM(L(II-M,S(M,F)),M,0,II-1)/II)$
 
EQ(II):=(H[II]+(DIFF(WGEN[II],T)+POISSON(WGEN[II],H[0]))/II
        +SUM(L(II-M,K[M])+M*S(II-M,H[M]),M,1,II-1)/II)$
 
LZAP(ANY):=DIFF(ANY,T)+POISSON(ANY,H[0])$

TRIGIDENTIFY(EXP):=IF FREEOF(SIN,EXP) AND FREEOF(COS,EXP) THEN 0 ELSE EXP$

NOTRESP(EXP):=NOT IS(LZAP(EXP) = 0)$

SETIFY(LIST):=(
   SET:[LIST[1]],
   FOR I:2 THRU LENGTH(LIST) DO(
      IF NOT MEMBER(LIST[I],SET) THEN SET:CONS(LIST[I],SET)),
   SET)$