# Computing normal forms with Mathematica

2009-10-01

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 = H0 + H1 + H2 + H3 + ... + Hn,

and the generator of the Lie series transformation

W = W1 + W2 + W3 + ... + Wk.

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

H = H0 + H1 + H2 + H3 + ... + Hn,

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