# A Maple procedure for computing the Campbell-Hausdorff formula # ============================================================== # by B. Keller (November 1998) # # # The procedure camph(o) computes the Campbell-Hausdorff formula up to # order o. # # The result is given as a polynomial in adX, adY and s. # The variable s only serves to keep track of the degree of the # monomials and to order them. Put s=1. The result is to be interpreted # as the operator H(X,?) which satisfies # # exp(H(X,Y))=exp(X)exp(Y) # # up to order o. # ################################################################## # # Here is some sample output # #> read `campha.txt`; #> campha(1); # s + 1/2 adX # #> campha(2); # 2 # s + 1/2 s adX - 1/12 adY adX + 1/12 adX adX # #> campha(3); # 3 2 # s + 1/2 s adX - 1/12 s adY adX + 1/12 s adX adX - 1/24 adX adY adX # #> campha(4); # 4 3 2 2 #s + 1/2 s adX - 1/12 s adY adX + 1/12 s adX adX - 1/24 s adX adY adX # # + 1/360 adY adX adX adX - 1/240 adX adY adX adX # # - 1/360 adY adY adX adX + 1/720 adY adY adY adX # # - 1/360 adX adY adY adX - 1/720 adX adX adX adX + 1/90 adY adX adY adX # # - 1/240 adX adX adY adX #> campha(5); #bytes used=23024124, alloc=1834672, time=31.53 #bytes used=24024532, alloc=1834672, time=32.87 # 5 4 3 3 2 #s + 1/2 s adX - 1/12 s adY adX + 1/12 s adX adX - 1/24 s adX adY adX # # + 1/360 s adY adX adX adX - 1/240 s adX adY adX adX # # - 1/360 s adY adY adX adX + 1/720 s adY adY adY adX # # - 1/360 s adX adY adY adX - 1/720 s adX adX adX adX # # + 1/90 s adY adX adY adX - 1/240 s adX adX adY adX # # + 1/180 adX adY adX adY adX + 1/720 adX adY adX adX adX # # - 1/480 adX adX adY adX adX - 1/720 adX adY adY adX adX # # - 1/720 adX adX adY adY adX + 1/1440 adX adY adY adY adX # # + 1/720 adX adX adX adY adX # elcomp:=proc(A,B) local F,P,Q,Q1,Q2,R,S; # # Elementary composition # Composition of two products containing each one occurrence # of a U[x,y,z ...] at most # Example : elcomp(2*t*U[1], 3*t^2*U[1,2])= 6*t^3*U[1,1,2] # The variables U[x,y,z, ...] are reserved. # if (not has(A,U)) or (not has(B,U)) then S:=A*B else P:=expand(A*subs(U=V,B)); F:=remove(has,P,{U,V}); Q:=select(has,P,{U,V}); Q1:=op(1,Q); Q2:=op(2,Q); if op(0,Q1)=`U` then Q:=U[op(Q1),op(Q2)] else Q:=U[op(Q2),op(Q1)] fi; if nops(Q)2 then T:=T+P fi; fi; od; T; end; eltransl:=proc(S) # # translates U[1,2,2,1] to adX adY adY adX # local i,w; w:=``; for i from 1 to nops(S) do if op(i,S)=1 then w:=cat(w,` `, `adX`) else w:=cat(w,` `, `adY`) fi; od; w; end; translate:=proc(S) # # translates all the terms U[1,2,2,1] in S to ad X adY adY adX # local i, T, P, L, Q; T:=S; for i from 1 to nops(S) do P:=op(i,S); if has(P,U) then Q:=select(has,P,U); L:=eltransl(Q); T:=subs(Q=L,T); fi; od; T; end; myorder:=proc(a,b); if degree(a,s)