Message Boards Message Boards

0
|
4240 Views
|
2 Replies
|
1 Total Likes
View groups...
Share
Share this post:

BHC formula : convert polynomial representation into a nested commutators

Hi,

I'm trying to convert polynomial representation of  z=log(exp(x)exp(y)) into nested commutators representation.
Polynomial representation can be obtained by Goldberg’s method.
(The code below is taken from Michael Weyraucha, Daniel Scholz, Computing the Baker–Campbell–Hausdorff series and the Zassenhaus product)
 g[1] = 1;
 
 g[s_] := g[s] = Expand[1/s*D[t*(t - 1)*g[s - 1], t]];
 c[w_] := c[w] = Module[
     {m, m1, m2, k},
     m = Length[w];
     m1 = Floor[m/2];
     m2 = Floor[(m - 1)/2];
     Integrate[t^m1*(t - 1)^m2*Product[g[w[[k]]], {k, m}], {t, 0, 1}]
     ];
BCH[n_Integer, alph_List] := Module[
   {p},
   p = Flatten[Permutations /@ IntegerPartitions[n], 1];
   Plus @@ (c[Sort[#]]*(words[#, alph] -
           (-1)^n*words[#, Reverse[alph]]) & /@ p) // Expand];
words[p_List, alph_List] := StringJoin @ (ConstantArray @@@
     Partition[Riffle[p, alph, {1, 2*Length[p], 2}], 2]);
For example:
In[2458]:= BCH[3, {"x", "y"}]
Out[2458]= ("xxy")/12 - ("xyx")/6 + ("xyy")/12 + ("yxx")/12 -("yxy")/6 + ("yyx")/12
The output is a polynomial of the form: P(x,y) = Sum a(x1,x2,..,xk) x1 x2 ... xk  with xi = x or y
In the same article (eq.14) the map that converts  polynomial representation into commutators is given:
F(P) = Sum a(x,y,x3,x4,...,xk)/(nx(x3,x4,...,xk)+1) [[...[[ [x,y],x3],x4],...,],xk],
here nx(x3,x4,...,xk) is a number of x in a monomial after xy, for e.g. n(xy)=1, n(xyx)=2, ...

Here is the algorithm I need to code:
1) for given n find BCH[n, {"x", "y"}]
2) remove all monomials (words) not starting with xy
3) for xy... monomials find a(x1,x2,..,xk) and nx(x3,x4,...,xk)
4) build  [[...[[ [x,y],x3],x4],...,],xk]

Example for z3:
1) ("xxy")/12 - ("xyx")/6 + ("xyy")/12 + ("yxx")/12 - ("yxy")/6 +("yyx")/122)  - ("xyx")/6 + ("xyy")/12
3)  - ("xyx")/6:    a(x,y,x)=-1/6 and n(x)=1
        ("xyy")/12:  a(x,y,y)=1/12 and n(y)=0
4) a(x,y,x)/(n(x)+1) [[x,y],x] + a(x,y,y)/(n(y)+1) [[x,y],y] = -1/12[[x,y],x]+1/12[[x,y],y]


Any tips on how to implement this map in Mathematica?
Thanks in advance,
I.M.
POSTED BY: Ivan Morozov
2 Replies
Here is what I've managed to do
 (* generate poly *)
 poly = BCH[3, {"x", "y"}] ;
 (* create list1 and list2 for a(x1,x2,...,xk) and x1x2...xk *)
 list1 = list2 = Table[0, {i, 1, Length[poly]}] ;
 Do[
   {
    list1[[i]] = poly[[i]] [[1]],
    list2[[i]] = poly[[i]] [[2]]
    },
  {i, 1, Length[poly]}
  ];
(* select words starting with xy *)
tmp = StringCases[list2, StartOfString ~~ "xy" ~~ ___] ;
(* find xy words positions and take corr. a's *)
A = list1[[Take[Position[tmp, x_String] // Flatten, {1, -1, 2}]]];
(* shape arr. for words *)
WORDS = tmp // Flatten;
(* build n(x3x4...)+1 array *)
B = StringCount[WORDS, "x"];
(* build the answer *)
answer = Table[0, {i, 1, Length[WORDS]}];
Do[
{
  answer[[i]] = A[[i]]/B[[i]] Apply[NonCommutativeMultiply,
     ToExpression[StringCases[WORDS[[i]], Repeated[_, 1]]], {0}]
  }, {i, 1, Length[WORDS]}
]
Total[answer]
I've used NonCommutativeMultiply to shape the answer into the form -(1/12) x ** y ** x + x ** y ** y/12,
but I still need to build correspondence between ** and commutators, i.e. to change x**y**x into [[[x,y],x],...] or ((x <my operation> y) < my operation>) < my operation> ... )

Can I redefine ** to be non associative or  define non associative (acting from left to right) <my operation> and replace ** with it, or replace ** with  <my operation> but correctly group the answer with parentheses?

I.M.
POSTED BY: Ivan Morozov
Here is a complete (not optimized) solution to the problem stated.I've checked terms up to z5, higher terms also seem to be correct from numerical experiments.
 (* M.Weyrauch,D.Scholz/Computer Physics Communications 180 (2009) \
 1558\[Dash]1565 *)
 (* Goldberg's method *)
 (* polynomial generator *)
 g[1] = 1;
 g[s_] := g[s] = Expand[1/s*D[t*(t - 1)*g[s - 1], t]];
 c[w_] := c[w] = Module[
     {m, m1, m2, k},
     m = Length[w];
    m1 = Floor[m/2];
    m2 = Floor[(m - 1)/2];
    Integrate[t^m1*(t - 1)^m2*Product[g[w[[k]]], {k, m}], {t, 0, 1}]
     ];
BCH[n_Integer, alph_List] := Module[
   {p},
   p = Flatten[Permutations /@ IntegerPartitions[n], 1];
   Plus @@ (c[Sort[#]]*(words[#, alph] -
           (-1)^n*words[#, Reverse[alph]]) & /@ p) // Expand];
words[p_List, alph_List] := StringJoin @ (ConstantArray @@@
     Partition[Riffle[p, alph, {1, 2*Length[p], 2}], 2]);
(* polynomial converter *)
CircleTimes[a_, b_, c__] = CircleTimes[CircleTimes[a, b], c];
PolyToCom[n_Integer] := Module[
   {
    POLY, LISTA, LISTB, TMP, A, B, WORDS, ANSWER, i
    },
   POLY = Apply[List, BCH[n, {"x", "y"}]];
   LISTA = POLY[[All, 1]];
   LISTB = POLY[[All, 2]];
   TMP = StringCases[LISTB, StartOfString ~~ "xy" ~~ ___] ;
   WORDS = TMP // Flatten;
   A = LISTA[[Take[Position[TMP, x_String] // Flatten, {1, -1, 2}]]];
   B = StringCount[WORDS, "x"];
   ANSWER = A;
   Do[
    {
     ANSWER[[i]] =
      A[[i]]/B[[i]] Apply[CircleTimes,
        ToExpression[StringCases[WORDS[[i]], Repeated[_, 1]]], {0}]
     }, {i, 1, Length[WORDS]}
    ];
   Total[ANSWER]
   ];
(* EG *)
PolyToCom[2]
PolyToCom[3]
PolyToCom[4]
PolyToCom[5]

I.M.
POSTED BY: Ivan Morozov
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract