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

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 `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 $(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
```