# Addendum to “Computing normal forms with Mathematica”

2016-11-14

Several years ago I had written here about a Mathematica implementation of the algorithm for computing normal forms for Hamiltonian systems. That post was written for people who are already familiar how the algorithm works and could figure out what functions like `rangePart`

are `inverseAd`

are doing. To make that implementation more accessible I give here an example of how to use it to compute the normal form for the Hénon-Heiles Hamiltonian given by
`rangePart`

and `inverseAd`

below.

First, we define the Hamiltonian as a list of homogeneous polynomials of successively higher degree.

Hpq = {1/2 (p1^2 + p2^2) + 1/2 (q1^2 + q2^2), q1^2 q2 - 1/3 q2^3};

Now, let

We define thew new coordinates by

pq = {p1, p2, q1, q2}; xy = {x1, x2, y1, y2}; xy2pq = {x1 -> (q1 - I p1)/Sqrt[2], x2 -> (q2 - I p2)/Sqrt[2], y1 -> (p1 - I q1)/Sqrt[2], y2 -> (p2 - I q2)/Sqrt[2]}; pq2xy = Simplify[Solve[MapThread[Equal, {xy, xy /. xy2pq}], pq][[1]]];

where `xy2pq`

gives the transformation from `pq2xy`

gives the inverse transformation. One can check that the given transformation is symplectic. In particular, in

PBxy[F_, G_] := Sum[ D[F, xy[[i]]] D[G, xy[[i + 2]]] - D[G, xy[[i]]] D[F, xy[[i + 2]]], {i, 2}]

The Hamiltonian in the new coordinates becomes

Hxy = Expand[Hpq /. pq2xy]

The quadratic part is

Simplify[PBxy[x1^a1 x2^a2 y1^b1 y2^b2, Hxy[[1]]]]

gives

I (a1 + a2 - b1 - b2) x1^a1 x2^a2 y1^b1 y2^b2

Let

ctAd = Function[{m}, I {1, 1, -1, -1}.Exponent[m, xy]];

give the number `rangePart`

and `inverseAd`

by

rangePartxy[F_] := If[F === 0, 0, Select[F, ctAd[#] != 0 &]]; inverseAdxy[F_] := If[F === 0, 0, Map[#/ctAd[#] &, F]];

The idea here is that `rangePartxy`

keeps the monomials `F`

that can be obtained through the Poisson bracket of some monomial `inverseAdxy`

gives the function `rangePartxy`

. Also note that these functions give the correct result only if their argument `F`

has been `Expand`

ed. The function `LieSeriesNormalization`

given in the earlier post takes care of that.

Now, everything is in place for the normalization. We can compute

{Nxy, Wxy} = LieSeriesNormalization[ Expand[Hxy], 4, PBxy, inverseAdxy, rangePartxy, Verbose -> False];

`Nxy`

is the resulting normal form, and `Wxy`

is the generating function. The number `4`

is the normalization order; here it means that we normalize the Hamiltonian up to terms of degree 6, since we start counting at degree 2 (the lowest degree terms in `Hxy`

).

We can go back to coordinates

Npq = Expand[Nxy /. xy2pq]

and the normal form can be obtained as one function through

Total[Npq]

The result I finally obtain from `Total[Npq] // InputForm`

is

p1^2/2 + (5*p1^4)/48 - (101*p1^6)/3456 + p2^2/2 + (5*p1^2*p2^2)/24 + (907*p1^4*p2^2)/1152 + (5*p2^4)/48 - (773*p1^2*p2^4)/1152 + (235*p2^6)/3456 + q1^2/2 + (5*p1^2*q1^2)/24 - (101*p1^4*q1^2)/1152 - (3*p2^2*q1^2)/8 + (65*p1^2*p2^2*q1^2)/64 - (73*p2^4*q1^2)/1152 + (5*q1^4)/48 - (101*p1^2*q1^4)/1152 + (263*p2^2*q1^4)/1152 - (101*q1^6)/3456 + (7*p1*p2*q1*q2)/6 + (161*p1^3*p2*q1*q2)/144 - (175*p1*p2^3*q1*q2)/144 + (161*p1*p2*q1^3*q2)/144 + q2^2/2 - (3*p1^2*q2^2)/8 + (263*p1^4*q2^2)/1152 + (5*p2^2*q2^2)/24 - (47*p1^2*p2^2*q2^2)/64 + (235*p2^4*q2^2)/1152 + (5*q1^2*q2^2)/24 + (65*p1^2*q1^2*q2^2)/64 - (47*p2^2*q1^2*q2^2)/64 + (907*q1^4*q2^2)/1152 - (175*p1*p2*q1*q2^3)/144 + (5*q2^4)/48 - (73*p1^2*q2^4)/1152 + (235*p2^2*q2^4)/1152 - (773*q1^2*q2^4)/1152 + (235*q2^6)/3456