Addendum to “Computing normal forms with Mathematica”

H = \frac{1}{2} (p_{1}^{2} + p_{2}^{2} + q_{1}^{2} + q_{2}^{2}) + q_{1}^{2} q_{2} - \frac{1}{3} 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};

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 $A F^{'} = 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