Legendre Transform of a function of a 3-vector

Question

I am trying to implement the Legendre transform of a function in Mathematica with the purpose of calculating the Hamiltonian of a system starting from the Lagrangian.

I have found this page which explained how to perform the transform on a function of a single variable. But I need to evaluate the transform of a function of a vector in 3D.

Here it is my clumsy attempt at it:

Off[Solve::ifun];
Off[InverseFunction::ifun];
Clear[legendreTransform];
legendreTransform[f_, v_, p_] := 
  Module[{d1 = D[f, v[[1]]], d2 = D[f, v[[2]]], d3 = D[f, v[[3]]],transform}, 
    transform /; (transform = (-f + v[[1]]*d1 + v[[2]]*d2 + v[[3]]*d3) 
    /.Solve[{p[[1]] == d1 && p[[2]] == d2 && p[[3]] == d3}, v[[1]], v[[2]], v[[3]]}] // Last;
    FreeQ[transform, Solve] && FreeQ[transform, InverseFunction])] /; ! (v === p); 
vel = {vel1, vel2, vel3};
mom = {mom1, mom2, mom3};
legendreTransform[Exp[Total[vel]], vel, mom]

This piece of code is not working. I was not able to troubleshoot it also because I don't fully understand the code in the page I linked in the first place. I looked up every word in the documentation, but I still cannot understand code flow.

If you could give my any hint or suggestion about where to look I would be grateful.

EDIT (SOLUTION):

The VariationalMethods path seemed me more appealing at first and in fact I used it to do my calculations. But it produced slightly wrong results in my case. It took me some time to realize it but there was a problem with a relative sign.

Eventually I wrote this function that takes the best of the two approaches (the direct and the alternate answer) and it is working as intended in my case:

legendreTransform[l_, q_, p_] := 
  Module[{mom = Flatten[{p}], pos = Flatten[{q}], h}, 
    First[h /.Quiet[Solve[h == D[pos,t].VariationalD[l[pos], D[pos, t], t] - l[pos] && 
    mom == VariationalD[l[pos], D[pos, t], t], Append[D[pos, t], h]], 
  {Solve::incnst, Solve::ifun}]]]

More insight on what I am doing: I am doing calculations in primordial cosmology and, as you may know, in General Relativity the Hamiltonian is null in value (but not in form). This is sometimes called "frozen formalism" because with a null Hamiltonian is trickier to define the evolution of the system in time. Maybe that is the reason the standard methods in VariationalMethods are failing me. More precisely, there is the wrong sign between the kinetical term and the potential term (the spatial curvature term) in the Hamiltonian produced by Legendre transform of the Lagrangian (I am working in the ADM formalism). This at first was not a big issue since all I did was introducing the correct sign by hand.

Answer 1

Direct answer:

The function you were trying to transform in the question is not a scalar and can therefore not be used as an example.

Here is a version of legendreTransform that works for one or multiple variables, with only the bare minimum of code so as to make it more readable (hopefully):

legendreTransform[f_, v_, p_] := 
 Module[{vel = Flatten[{v}], mom = Flatten[{p}], h},
  First[h /. 
    Quiet[Solve[h == vel.Grad[f, vel] - f && mom == Grad[f, vel], 
      Append[vel, h]], {Solve::incnst, Solve::ifun}]]]

It leaves the hard work to Solve, with the relevant warnings turned off only inside the function, not globally. The arguments for the variables can be symbols (in one dimension) or lists of symbols (in two or more dimensions). Both cases are reduced to the same code by wrapping even a 1D variable in a List (using Flatten) so that it can be used in Grad. Of course the names for the components of velocity and momentum must be different, but are otherwise arbitrary. For the sign of the result, I used the convention of Wikipedia, which agrees with the one in the question.

The first example reproduces the one-dimensional test case:

legendreTransform[Exp[v], v, p]

$p (\log (p)-1)$

The second example is a toy Lagrangian in three dimensions:

legendreTransform[
   1/2 (v1^2 + v2^2 + v3^2) - (v1 - v2)/2, {v1, v2, v3}, {p1, p2, p3}]

$$\frac{1}{4} \left(2 \text{p1}^2+2 \text{p1}+2 \text{p2}^2-2 \text{p2}+2\text{p3}^2+1\right)$$

Alternative answer:

Given that you're interested in going from a Lagrangian to a Hamiltonian, you could use a generalization of this answer:

Needs["VariationalMethods`"]

Clear[pos, mom, r, p, t, h, V]
pos[t_] := Table[r[i][t], {i, 3}]
mom[t_] := Table[p[i][t], {i, 3}]

L = 1/2 m  pos'[t].pos'[t] - V[pos[t]];

Expand[h /. First@Solve[
  mom[t] == VariationalD[L, pos'[t], t] && h == FirstIntegral[t] /. 
   FirstIntegrals[L, pos[t], t], Append[pos'[t], h]]]

This uses the VariationalMethods package and therefore I had to define positions pos and momenta mom as functions of the time t.

Then I introduced a standard Lagrangian L from which you get to the Hamiltonian by Legendre transformation. But instead of doing the transformation by hand, I let FirstIntegrals find the Hamiltonian in terms of the positions and velocities. Then I just have to substitute the momenta for the velocities, using the fact that momenta are the variational derivative of L with respect to velocities pos'[t].

The substitution is done in Solve. I ask it to solve for the unknown Hamiltonian h, and also for the velocities. Although I then discard the velocities, this strategy makes sure that Solve finds a unique answer in terms of the variables I want to keep, i.e., the positions and momenta.