10

I want to go from Lagrangian description to Hamiltonian one. Using the example given by Mathematica I do something like:

genCoords = {x[t]};
T = 1/2 m x'[t]^2;
V = 1/2 K x[t]^2;
Q = -c x'[t];
L = T-V;

Now I want the Hamiltonian which is given by the relation : $$p = \frac{\partial L}{\partial \dot{q}}$$ and $$H=\dot{q}\frac{\partial L}{\partial \dot{q}}-L$$

The result I want is H as a function of $q$ and $p$ and I don't understand how to calculate the set {p[t]} nor to replace it inside the expression of the Hamiltonian.

VividD
  • 3,660
  • 4
  • 26
  • 42
sailx
  • 203
  • 1
  • 6

2 Answers2

16

Here is how you would do it using the standard add-on package VariationalMethods, which is meant for calculations like this:

Clear[m, k, c, x, t];
T = 1/2 m x'[t]^2;
V = 1/2 k x[t]^2;
L = T - V;

Needs["VariationalMethods`"]

hamiltonianEq = 
 h == FirstIntegral[t] /. Last[FirstIntegrals[L, {x[t]}, t]]

(* ==> h == 1/2 (k x[t]^2 + m Derivative[1][x][t]^2) *)

momentumEq = p[x] == VariationalD[L, x'[t], t]

(* ==> p[x] == m Derivative[1][x][t] *)

h /. 
 First[Solve[Eliminate[{hamiltonianEq, momentumEq}, x'[t]], h]]

(* ==> (p[x]^2 + k m x[t]^2)/(2 m) *)

Here, I defined two equations defining the Hamiltonian and momentum using the capabilities of the VariationalMethods package, and then used Eliminate to combine these equations into the final result.

The result is in simplified form, but you can also split it into two terms (which is more standard), by adding Apart@ in front of it.

For time-dependent Hamiltonians:

The above works for Hamiltonians that are equal to the conserved energy, as is the case in the question. However, it's easy to modify the definition of hamiltonianEq above so that it also works in general, following the definition given in the question. You would replace the first equation by

hamiltonianEq = h == x'[t] VariationalD[L, x'[t], t] - L

The rest stays the same.

Jens
  • 97,245
  • 7
  • 213
  • 499
10
genCoords = {x[t]};
ke = 1/2 m x'[t]^2;
v = 1/2 k x[t]^2;
q = -c x'[t];
l = ke - v;

Solve for x'[t] in terms of p[t]:

rule = First@Solve[p[t] == D[l, x'[t]], x'[t]]
(* {x'[t] -> p[t]/m} *)

and then replace this expression into the Hamiltonian:

x'[t] D[l, x'[t]] - l /. rule
(* p[t]^2/(2 m) + 1/2 k x[t]^2 *)

Note that I have used lower-case symbols for all of your variables. This is to avoid conflicts with built-in Mathematica symbols like K, which you used. This is good Mathematica practice.

march
  • 23,399
  • 2
  • 44
  • 100
  • 2
    This is the simplest answer, but it hasn't kept me from adding a more complicated alternative which may be of interest in more general situations. Especially for people who don't remember the definitions and don't want to look them up... – Jens Apr 19 '16 at 17:43
  • 2
    @Jens, ...of which there are a frightening lot. :o – J. M.'s missing motivation Apr 20 '16 at 05:26