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.

Suppose that you have a Hamiltonian function

$$ H = H_0 + H_1 + H_2 + H_3 + \cdots + H_n, $$and the generator of the Lie series transformation

$$ W = W_1 + W_2 + W_3 + \cdots + 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.**

The normalization algorithm is very similar except for the fact that now we
don´t know the generating function `W`

and we will have to compute it along the
way. Suppose again that we have a Hamiltonian function

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]
```