# How to define the Poisson bracket in Mathematica

If $f,g$ are functions of the independent variables $\{q_1, q_2, ..., q_N, p_1, p_2, ..., p_N\}$, then the Poisson bracket is defined as:

\begin{equation*} [f, g] = \sum_{i=1}^N \left(\frac{\partial f}{\partial q_i}\frac{\partial g}{\partial p_i} - \frac{\partial f}{\partial p_i}\frac{\partial g}{\partial q_i}\right) \end{equation*}

I'd like to define Poisson bracket in Mathematica. Just so to prove I'm not lazy, I wrote the following snippet (I doubt it is correct, let alone slick).

ClearAll["Global*"];
pbra[f_, g_, q_, p_] :=
With[{c = Join[q, p]},
Sum[
D[f@@c, q[[k]]] D[g@@c, p[[k]]] - D[f@@c, p[[k]]] D[g@@c, q[[k]]],
{k, 1, Length@q}]]


And call it like (dummy example):

f[x_, y_, px_, py_] := x+y+px*py;
h[x_, y_, px_, py_] := x-y+px/py;
p[f, g, {x, y}, {px, py}]


How would you define Poisson bracket ? For example, I am wondering whether it would be more sound to define it in a way that the user calls it like:

f[x_, y_, px_, py_] := x+y+px*py;
h[x_, y_, px_, py_] := x-y+px/py;
p[f[x, y, px, py], g[x, y, px, py]]


This is how I have them defined in a very old notebook of mine. I believe I took the code - hopefully adapted, if not straigthforward copied - by Gerd Baumann's "Mathematica in Theoretical Physics". The first edition in one volume.

For a single pair of canonical coordinates:

PoissonBracket[a_, b_, q_Symbol, p_Symbol] :=
Simplify[D[a, q]  D[b, p] - D[b, q ] D[a, p] ]


For several coordinates

PoissonBracket[a_, b_, q_List, p_List] := Block[{pk, n},
n = Length[q];
If[n == Length[p],
pk = Simplify[Sum[
D[a, q[[j]]] D[b, p[[j]]] - D[b, q[[j]]] D[a, p[[j]]],
{j, 1, n}]
],
Print["Incompatible lengths"]]]


Here's an example

H = p^2/(2m) + k^2 q^2/2;
PoissonBracket[q, H, q, p]
PoissonBracket[p, H, q, p]


gives

p/m

-k^2 q

which translates in the motion equations

eqs = { q'[t] == p[t]/m,
p'[t] == -k^2 q[t] }

• Thank you @Peltio! I see that Gerd Baumann is the same person that wrote the article @rasher mentioned earlier. Do you suggest getting that book ? Commented Feb 7, 2014 at 5:17
• It's been a long time since I read it. I remember it had several very interesting approaches on how to solve physics problems with mma. I also remember that I hated its graphical layout (something that was probably due to the publisher and not the author) - I do remember very clearly Poincarè sections and phase portraits that were painful to watch for the dimension of the points and the thickness of the lines used... So, I suggest you try to get a look at it in a library before making up your mind. Also, there is a second edition in 2 vols, now. Commented Feb 7, 2014 at 5:33

Just another approach:

PoissonBracket[f_, g_, q_List, p_List] /; Length[q] == Length[p] :=
Fold[Plus, 0, MapThread[D[f, #1] D[g, #2] - D[f, #2] D[g, #1] &, {q, p}]]


Then:

H = p^2/(2 m) + k^2 q^2/2;

PoissonBracket[p, H, {q}, {p}]

(* -k^2 q *)


Edit

I think the following might be faster since Total is called once on the entire List

PoissonBracket2[f_, g_, q_List, p_List] /; Length[q] == Length[p] :=
Total @ Flatten @ Fold[List, 0, MapThread[D[f, #1] D[g, #2] - D[f, #2] D[g, #1] &, {q, p}]]


You could also do:

PoissonBracket[f_, g_, q_List, p_List] /; Length[q] == Length[p] :=
D[f, {q}].D[g, {p}] - D[f, {p}].D[g, {q}]


and if needed add the following to accommodate non-lists:

PoissonBracket[most__, q : Except[_List], p_] :=
PoissonBracket[most, {q}, p]

PoissonBracket[most__, p : Except[_List]] :=
PoissonBracket[most, {p}]


Example:

H = p^2/(2 m) + k^2 q^2/2;
PoissonBracket[p, H, q, p]
`

-k^2 q