Computing normal forms with Mathematica

I have been unhappy for some time now with the function that I had written in Mathematica for computing the normal form of a Hamiltonian system. I used the standard Lie series approach but I had written the code using a traditional approach with several Do loops. Except for the fact that this kind of code is completely contrary to Mathematica´s philosophy I found the code hard to read and inefficient to the degree at least that Mathematica can be efficient.

So, last Sunday I decided to revise this code. I worked out again from scratch the Lie series transformation algorithm and I quickly figured out that I could write a very simple, completely declarative, implementation of this algorithm in Mathematica.

H = H_{0} + H_{1} + H_{2} + H_{3} + \dots + H_{n},

W = W_{1} + W_{2} + W_{3} + \dots + W_{k} .

Then the following Mathematica function takes as arguments the lists {H0, H1, ..., Hn} and {W1, W2, ..., Wk}, and the Poisson bracket PB. If, for example, the functions H and W are expressed in canonical coordinates qp = {q1, ..., qm, p1, ..., pm} then PB can be defined as

PB[F_,G_] := Sum[D[F,qp[[i]]] D[G,qp[[i+m]]] - D[G,qp[[i]]] D[F,qp[[i+m]]], {i,m}];

The function returns the transformed Hamiltonian as a list {N0, N1, ..., Nk}. The transformation algorithm is

LieSeriesTransformation[Hlist_, Wlist_, pb_, 
  OptionsPattern[{MaxOrder -> 0, Verbose -> True}]] := Module[
  {order, verbose, H, W, L, N},

  If[Head[OptionValue[MaxOrder]] =!= Integer || 
    OptionValue[MaxOrder] < 0, 
   Message[LieSeriesTransformation::maxord]; Return[]];
  If[OptionValue[MaxOrder] == 0, order = Length[Wlist], 
   order = OptionValue[MaxOrder]];
  verbose = OptionValue[Verbose];
  If[verbose =!= True && verbose =!= False, 
   Message[LieSeriesTransformation::verbose, verbose]; Return[]];
    
  W[s_] := W[s] = If[s <= Length[Wlist], Wlist[[s]], 0];

  L[0, s_] := L[0, s] = If[s + 1 <= Length[Hlist], Hlist[[s + 1]], 0];

  L[k_, m_] := L[k, m] = Expand[Sum[1/k pb[L[k - 1, m - l], W[l + 1]], {l, 0, m}]];
  N[s_] := N[s] = Sum[L[n, s - n], {n, 0, s}];

  Table[
   If[verbose, Print[StringForm["Order `1`", s]]];
   N[s], {s, 0, order}]]

The first few lines in the function just take care of options passed to the function. Then L[0,s] and W[s] are set to the elements of the lists that we passed as arguments. The core of the algorithm are the following two lines where the rules are given for computing L[k,m] and the transformed Hamiltonian N[s]. Finally, a list of N[s] for all values of s up to the required order is created. It is only at this step that anything is actually computed. When we ask for, let´s say, N[4] the algorithm starts applying the rules for N[s] and L[k,m] computing along the way everything that is necessary. Pure beauty.

H = H_{0} + H_{1} + H_{2} + H_{3} + \dots + H_{n},

and we want to normalize it up to some prescribed order with respect to H0. The algorithm is

LieSeriesNormalization[Hlist_, order_, pb_, inverseAd_, rangePart_, 
     OptionsPattern[{Verbose -> True, 
    ReturnGeneratingFunction -> True}]] := 
   Module[{ N, W, adW, L, Nlist, Wlist, verbose, returngenfunc},
  
  verbose = OptionValue[Verbose];
  returngenfunc = OptionValue[ReturnGeneratingFunction];
  If[verbose =!= True && verbose =!= False,
   Message[LieSeriesNormalization::boolarg, "Verbose", verbose];
   Return[]];
  If[returngenfunc =!= True && returngenfunc =!= False,
   Message[LieSeriesNormalization::boolarg, 
    "ReturnGeneratingFunction", returngenfunc];
   Return[]];
  
  L[0, s_] := L[0, s] = If[s + 1 <= Length[Hlist], Hlist[[s + 1]], 0];
  L[k_, m_] := L[k, m] = Expand[Sum[pb[ L[k - 1, m - l], W[l + 1]]/k, {l, 0, m}]];

  N[0] = L[0, 0];
  N[s_] := N[s] = Module[{rp, temp},
     If[verbose, Print[StringForm["Order `1`", s]]];
     W[s] = 0;
     temp = Sum[L[n, s - n], {n, 0, s}];
     rp = rangePart[temp];
     W[s] = inverseAd[Expand[-rp]];
     L[1, s - 1] = Expand[L[1, s - 1] - rp];
     Expand[temp - rp]
     ];
  
  Nlist = Table[N[s], {s, 0, order}];
  If[verbose, Print[StringForm["Finished normalization"]]];
  If[returngenfunc,
   Wlist = Table[W[s], {s, 1, order}]; {Nlist, Wlist},
   Nlist]
  ]

The algorithm is again concise, but not as elegant as the previous one. The reason is that in order to compute N[s] we must first compute W[s]. The easiest way to do this is to set W[s] = 0, compute N[s] (which is called temp in the algorithm and does not in general commute with H0 so it is not the N[s] that we are after), and then find the W[s] that kills the part rp of N[s] that is in the range of PB[H0,_]. The remaining part temp - rp is the required N[s]. Some people may find this little ‘dance’ clever, but it is much less elegant than the previous algorithm for the Lie series transformation. I have not thought yet of an efficient, completely declarative, way to write this algorithm.

As a test, I used the last algorithm to compute the normal form for the perturbation of a three degree of freedom 1:1:2 resonant oscillator. The computation up to terms of order 12 (10 normalization steps) in canonical coordinates q, p took on my MacBook approximately 1 minute and up to order 14 (12 normalization steps) approximately 6 minutes. The FORTRAN crowd will probably find this extremely long but for my purposes it is more than satisfactory.

Finally, here is a very cool trick. If you want to see the result of the Lie series transformation in purely symbolic form run the following code after making sure that H0, H1, …, W1, W2, …, and PB are not already defined.

PB[a_Integer F_, G_] := a PB[F, G];
PB[a_Rational F_, G_] := a PB[F, G];
PB[F_ + G_, H_] := PB[F, H] + PB[G, H];
NH = LieSeriesTransformation[{H0, H1, H2, H3, H4}, {W1, W2, W3, W4}, PB]