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

$$ H = \frac12 (p_1^2 + p_2^2 + q_1^2 + q_2^2) + q_1^2 q_2 - \frac13 q_2^3. $$

It is important to notice that the quadratic part of the Hénon-Heiles Hamiltonian is a resonant oscillator so we will be computing here a resonant normal form. This affects the definition of 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 $A$ denote the operator given by $A[F] = \\{F,H_2\\}$ where $F$ is a smooth function of $(p,q)$ and $H_2$ is the quadratic part of $H$. Computations become easier if we can express $H_2$ in coordinates $(x,y)$ in which $A$ acts diagonally on monomials $m=x_1^{a_1} x_2^{a_2} y_1^{b_1} y_2^{b_2}$. That means that we want to have coordinates such that

$$ A[m] = c(m) m, $$

where $c(m)$ is a (complex) number that depends on $m$.

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 $(x,y)$ coordinates to $(p,q)$ and pq2xy gives the inverse transformation. One can check that the given transformation is symplectic. In particular, in $(x,y)$ coordinates we have the standard Poisson bracket

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 $H_2 = i(x_1 y_1 + x_2 y_2)$ and we check that $A$ acts diagonally on monomials, that 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 $c(m) = i(a_1+a_2-b_1-b_2)$ for the monomial $m = x_1^{a_1} x_2^{a_2} y_1^{b_1} y_2^{b_2}$. Then we can define 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 $m$ of F that can be obtained through the Poisson bracket of some monomial $m'$ with $H_2$, while inverseAdxy gives the function $F'$ for which $AF'=\\{F',H_2\\}=F$. For this to work we must have that $F$ is in the range of the operator $A$ and this is why it is first necessary to use rangePartxy. Also note that these functions give the correct result only if their argument F has been Expanded. 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 $(p,q)$ as follows,

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