How to define the Poisson bracket in Mathematica

Question

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]]

Answer 1

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

Answer 2

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

Answer 3

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