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};
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