3
$\begingroup$

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!

$\endgroup$
  • 1
    $\begingroup$ 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[1], p[1] evaluates to a long expression -- Mathematica can't hope to differentiate with respect to that. $\endgroup$ – jjc385 Jun 7 '17 at 0:28
  • $\begingroup$ @jjc385 I can see that. However, do you have any idea of how to fix that? $\endgroup$ – Rumplestillskin Jun 7 '17 at 0:33
  • $\begingroup$ @jjc385 I wish that was the issue. That was just a typo. $\endgroup$ – Rumplestillskin Jun 7 '17 at 3:49
3
$\begingroup$

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 )

__

Edit2/ Significant Addition:

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"]
$\endgroup$
  • $\begingroup$ 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!! $\endgroup$ – Rumplestillskin Jun 7 '17 at 1:51
  • $\begingroup$ @Rumplestillskin See my edit. Now NDSolve is producing output. $\endgroup$ – jjc385 Jun 7 '17 at 4:46
  • $\begingroup$ 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! $\endgroup$ – Rumplestillskin Jun 7 '17 at 6:17
  • 1
    $\begingroup$ @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. $\endgroup$ – jjc385 Jun 7 '17 at 16:06
  • $\begingroup$ @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. $\endgroup$ – jjc385 Jun 7 '17 at 16:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.