# How to formulate the Hamiltonian equations of motion? Part II

I've recently asked a question on an issue I was facing with numerically integrating Hamiltonian equations of motion. I got a great answer.

Following on from this, I wanted to write a very similar code using alternate definitions in variational mechanics.

That is using the following:

1. The Lagrangian is given by $L= \frac{1}{2} g_{\mu\nu} \dot{x}^\mu(s) \dot{x}^\nu(s)$
2. The generalised coordinates are given by $x^\mu(s)$ with velocity given by $\dot{x}^\mu(s)$.
3. The conjugate momentum is given by $p = \partial L\;/ \;\partial \dot{x}^\mu(s)$
4. The Hamiltonian is given by $H=p\dot{x} - L$ or $H=\frac{1}{2} g^{\mu\nu} p_\mu p_\nu$.
5. Hamilton's equations are given by $\dot{p} = -g^{\mu\nu}\partial_q H, \; \dot{q} = g^{\mu\nu}\partial_p H$

Now, with the help of @jjc385 I think it is almost there. however, I still can't seem to finalise the equations of motion. I know they are getting there because I can compare them with the geodesic equations of motion for the same problem. (There is a reason why I'm being awkward and proceeding this way). The updated code thanks to @jjc385 is now given by:

SetAttributes[m, Constant];
q = {t[s], r[s], \[Theta][s], \[Phi][s]};
vel = D[q, s];
n = Length[q];
tt = 1 - 2 m/r[s];
rr = -1/tt;
\[Theta]\[Theta] = -r[s]^2;
\[Phi]\[Phi] = -(r[s] Sin[\[Theta][s]])^2;
metric = {{tt, 0, 0, 0}, {0, rr, 0, 0}, {0, 0, \[Theta]\[Theta],
0}, {0, 0, 0, \[Phi]\[Phi]}};
inversemetric = Simplify[Inverse[metric]];

p = FullSimplify[Table[D[L, vel[[i]]], {i, 1, n}]];

H = FullSimplify[1/2*(p.inversemetric.p)];
pSym = Symbol@*(StringJoin["p" <> #] &)@*ToString@*Head /@ q;
Solve[pSym == p, vel]
Hnew = H /. Flatten@%;

ivsnew = {1.3, 0, 0, 0.088}; ics = {0, 6.5, \[Pi]/2, 0};
m = 1;

pdot = FullSimplify[inversemetric.Table[-D[H, q[[i]]], {i, 1, n}]]
qdot = inversemetric.Table[D[Hnew, pSym[[i]]], {i, 1, n}]
eqs1 = {{D[q, s] == qdot,
D[p, s] == pdot}, {(q /. s -> 0) == ics, (p /. s -> 0) ==
ivsnew}};
time = {s, 0, 750};
solee = NDSolve[eqs1, q, time, Method -> "ExplicitRungeKutta"];


Any suggestions?

I know the equations of motion are incorrect but I don't know how to make them correct!

• You should be able to tell that pdot=... works without a problem, but qdot=... does not. As Mathematica is trying to tell you with the error messages, when you try to differentiate H with respect to p, p evaluates to a long expression -- Mathematica can't hope to differentiate with respect to that. – jjc385 Jun 7 '17 at 0:28
• @jjc385 I can see that. However, do you have any idea of how to fix that? – Rumplestillskin Jun 7 '17 at 0:33
• @jjc385 I wish that was the issue. That was just a typo. – Rumplestillskin Jun 7 '17 at 3:49

First, I think p=... has a typo, and should read

p = FullSimplify[Table[D[L, vel[[i]]], {i, 1, n}]]


One way to use Hamilton's equations of motions is to solve for the (symbolic) conjugate momenta in terms of q[[i]]:

pSym = Symbol@*(StringJoin["p" <> #] &)@*ToString@*Head /@ q

{pt, pr, pθ, pϕ}

Solve[pSym == p], vel]

{{ t'[s] -> pt * ... , r'[s] -> pr * ..., ... }}


Now you can write H as an explicit function of the symbolic conjugate momenta:

Hnew = H /. Flatten@%


Finally, you can do

qdot = inversemetric.Table[D[Hnew, pSym[[i]] ], {i, 1, n}]


which will be a function of the elements of vel and pSym. Hopefully the equation solving will proceed as you wish.

Edit1: You might find the last line simpler using either of the following, which all give the same output as the original:

qdot = inversemetric.Table[ D[Hnew, pi], {pi, pSym} ]
qdot = inversemetric.( D[Hnew,#]& /@ pSym )


__

When using Hamilton's equations of motion, your independent variables are $p$ and $q$ -- $\dot{q}$ should not appear at all.

First, it helps to write the components of pSym as functions of s, like you did for q :

pSym = #[s] & /@ (Symbol@*(StringJoin["p" <> #] &)@*ToString@*Head /@ q)

{pt[s], pr[s], pθ[s], pϕ[s]}


You should find pdot in terms of $p$ rather than $\dot{q}$ :

pdot = FullSimplify[inversemetric.Table[-D[Hnew, q[[i]]], {i, 1, n}]]


Then you can integrate the equations of motion:

eqs = { {D[pSym, s] == pdot, D[q, s] == qdot},
{(q /. s -> 0) == ics, (pSym /. s -> 0) == ivsnew}  };
NSolve[ eqs, time, Method -> "ExplicitRungeKutta" ]


I'm not 100% confident with the output, but at least it's giving something that seems somewhat reasonable.

Update: The full code I'm using:

SetAttributes[m, Constant];
q = {t[s], r[s], \[Theta][s], \[Phi][s]};
vel = D[q, s];
n = Length[q];
tt = 1 - 2 m/r[s];
rr = -1/tt;
\[Theta]\[Theta] = -r[s]^2;
\[Phi]\[Phi] = -(r[s] Sin[\[Theta][s]])^2;
metric = {{tt, 0, 0, 0}, {0, rr, 0, 0},
{0, 0, \[Theta]\[Theta], 0}, {0, 0, 0, \[Phi]\[Phi]}};
inversemetric = Simplify[Inverse[metric]];
L = FullSimplify[1/2 vel.metric.vel];

p = FullSimplify[Table[D[L, vel[[i]]], {i, 1, n}]];

H = FullSimplify[1/2*(p.inversemetric.p)]
pSym = #[s] & /@ (Symbol@*(StringJoin["p" <> #] &)@*ToString@*Head /@ q)
Solve[pSym == p, vel]
Hnew = H /. Flatten@%

ivsnew = {1.3, 0, 0, 0.088}; ics = {0, 6.5, \[Pi]/2, 0};
m = 1;

pdot = FullSimplify[inversemetric.Table[-D[Hnew, q[[i]]], {i, 1, n}]]
qdot = inversemetric.Table[D[Hnew, pSym[[i]]], {i, 1, n}]
time = {s, 0, 750};
eqs = {{D[pSym, s] == pdot,
D[q, s] == qdot}, {(q /. s -> 0) == ics, (pSym /. s -> 0) ==
ivsnew}};
NDSolve[eqs, pSym, time, Method -> "ExplicitRungeKutta"]

• It is almost perfect with the help of that great answer. However, I am still coming up short! I can see the equations of motion are beginning to take the same form as if I were to approach the problem using the geodesic equations of motion. I've updated the question hopefully for the final time. I'd love if you could take a look. Your answer helped a lot!! – Rumplestillskin Jun 7 '17 at 1:51
• @Rumplestillskin See my edit. Now NDSolve is producing output. – jjc385 Jun 7 '17 at 4:46
• I don't know why but the output is still wrong!!! This is far more challenging that I first though. Thanks again for the output. The r[s] variable should be period with a stable orbit!! I haven't got the faintest clue why this isn't working! – Rumplestillskin Jun 7 '17 at 6:17
• @Rumplestillskin (Ignore the now-deleted comment about the code not working -- I was accidentally using NSolve rather than NDSolve.) My most recent edit was just to fix a typo in eqs, and post my full code. – jjc385 Jun 7 '17 at 16:06
• @Rumplestillskin Note that one can find explicit expressions for $\dot{q}$ and $\dot{p}$ very easily by hand: $\dot{q}^\nu = p_\mu g^{\mu\nu}$ and $\dot{p}_\mu = -\frac{1}{2} p_\alpha p_\beta \partial_\mu g^{\alpha\beta}$. I calculated pdot and qdot with these expressions directly, and I seem to be getting the same numerical result. Anyway, I don't think I have time to work on this further for now. Good luck. – jjc385 Jun 7 '17 at 16:37