# A EoM in Schwarzschild geometry: geodesic v Hamilton formalism

Tags:
1. Jun 5, 2017

### Matter_Matters

Hi there guys,

Currently writing and comparing two separate Mathematica scripts which can be found here and also here. The first one I've slightly modified to suit my needs and the second one is meant to reproduce the same results.

Both scripts are attempting to simulate the trajectory of a timelike geodesic using
1. the geodesic equations of motion.
2. The Hamiltonian formalism for equations of motion.
Now, the results should be equivelent. However, my problem is as follows:

The Hamiltonian produces a timelike geodesic that falls radially inward directly to the event horizon, where the geodesic equation approach produces a beautiful precessing orbit. Now, the approach using the geodesic equation is correct, however, can anyone see why this is different in the second link above? I have no idea why they are not producing the same results!!

2. Jun 5, 2017

### Paul Colby

Do the two approaches produce the same equations of motion? This could be directly checked in Mathematica.

3. Jun 5, 2017

### Matter_Matters

Yes the output should be complete identical however it's not! I have been messing around for a couple of hours with no success unfortunately!

4. Jun 5, 2017

### Paul Colby

Wait, same equations, same solver, same boundary conditions yield a different answer? Hum, I doubtful. Something must not be identical. My suggestion was to reduce the Hamiltonian equations of motion down to the geodesic ones and test for difference before solving them. I suspect that different equations result and the hamiltonian derived ones are not equivalent to the geodesic ones. Your problem is very likely a programming or problem setup one not one and not with the physics.

in the hamiltonian approach the momentum is given by $p = \dot q$ which is not $p = \frac{\partial L}{\partial\dot{q}}$

5. Jun 5, 2017

### pervect

Staff Emeritus
As I recall, to derive the Hamiltonian approach, we can start out with some Lagrangian $L(q,\dot{q})$ and we DEFINE $p = \partial L / \partial \dot{q}$. The we compute the energy function $h(q,\dot{q})$, and use it to find the Hamiltonian H(p,q) with a change of variables. (Alternatively, we can define H as the Legendre transformation of L). So I'm not sure what point you're making. We write one of Hamilton's equations as $dq/dt = \partial H / \partial p$, but this isn't the same as what you wrote above. I'm not sure where you got $p = \dot{q}$, this doesn't look right. Wiki seems to bear me out, <<link>>.

6. Jun 5, 2017

### Paul Colby

In the link provided in the OP (the second "here"), $p = \dot{q}$ which is very much not $\frac{\partial L}{\partial\dot{q}}$ at least by my reading of the Mathematica code provided. My point was perhaps this is an error.

7. Jun 5, 2017

### pervect

Staff Emeritus
$p = \frac{\partial S}{\partial q}$ could make sense in the Hamilton-Jacobi approach (henceforth the HJ approach), but that doesn't seem to be what the script is saying, if I'm interpreting it correctly it seems indeed saying that $p = dq/ds$. This equation isn't making any sense to me, sorry.

I'd suggest looking up the equations of motion in a text, rather than trying to reverse engineer them from a pair of Mathematica scripts, one of which isn't working correctly.

Attempting to follow my own advice, I do see a fair amount of discussion of the Hamilton-Jacobi equations of motion in my text (MTW), but nothing that mentions the "Hamiltionian" method.

I'm not sufficiently familiar with either the HJ method or Mathematica to be of any more help, unfortunately - if indeed the OP is interested in the HJ method ( something I'm not totally convinced of either).

I'm fairly familiar with the geodesic method, but that seems to be working correctly and there aren't any questions about it.

8. Jun 5, 2017

### Matter_Matters

Hmm so I've got the following code. However, the output is still not corresponding to the geodesic equations approach:

[NOPARSE]
Code (Text):

(*Symplectic schemes*)
Clear[m];
q = {t[s], r[s], θ[s], ϕ[s]};
p = {pt[s], pr[s], pθ[s], pϕ[s]};
n = Length[q];
tt = 1 - 2 m/r[s];
rr = -1/tt;
θθ = -r[s]^2;
ϕϕ = -(r[s] Sin[θ[s]])^2;
metric = {{tt, 0, 0, 0}, {0, rr, 0, 0}, {0, 0, θθ, 0}, {0, 0, 0, ϕϕ}};
inversemetric = Simplify[Inverse[metric]];
ivs = {1.3, 0, 0, 0.088}; ics = {0, 6.5, π/2, 0};
m = 1;

(*hamiltonian=1/2 (D[q,s].inversemetric.D[q,s])-1;*)

hamiltonian = 1/2 (p.inversemetric.p) - 1;        (* define H in terms of p, q *)
pdot = Table[-D[hamiltonian, q[[i]]], {i, 1, n}];
qdot = Table[D[hamiltonian, p[[i]]], {i, 1, n}];

eqs1 = {{D[q, s] == qdot, D[p, s] == pdot},       (* adjust equations to new defs *)
{(q /. s -> 0) == ics, (p /. s -> 0) == ivs}};

invariants = hamiltonian;
time = {s, 0, 100};                               (* change T to s to match q etc. *)
solee = NDSolve[eqs1, q, time, Method -> "ExplicitRungeKutta",
StartingStepSize -> 0.1]
[/NOPARSE]
It's physical output is a particle that falls almost directly in to the event horizon of the Schwarzschild black hole. I have no idea why it does this though.

9. Jun 5, 2017

### Paul Colby

when $p$ is eliminated from the Hamiltonian system one must obtain the geodesic equation up to a possible parameter change. Once this is verified symbolically (not as a numerical solution but in terms of symbols) the answers will be the same. If they aren't symbolically equivalent then there is an error.

10. Jun 6, 2017

### Matter_Matters

I believe the code that I posted above should be doing the job. The EoM are given by Hamilton's equations and I've edited the code in line with the above comments. Still though, I'm encountering the singularity almost at the beginning of the simulation which I can't seem to understand.

11. Jun 6, 2017

### Paul Colby

Well, debugging others code is a good method to learn a language. Your inverse metric is poorly behaved at $\theta = 0$. Okay you're at $\theta = \pi/2$

Last edited: Jun 6, 2017
12. Jun 6, 2017

### Paul Colby

Okay, I get a stiff system warning at $s=10.8$ which isn't too surprising.

13. Jun 6, 2017

### Matter_Matters

Paul, I don't think this physically makes sense which is how I know there is an issue somewhere. It should have a stable orbit with those IC's.

14. Jun 6, 2017

### Paul Colby

Well, the $\phi(s)$ changes linearly with $s$ which is clearly wrong since angular momentum should lead to an ever increasing orbital rate. Also, there are constraints on the on the initial world velocity which seems broken. I did some rummaging around in the geodesic approach (your first "here" link) looking for the metric you are using. I couldn't find it. I agree with the form of the hamiltonian used but if the metric is non-physical all bets are off. Again, I suggest you forget solving and compare the equations after removing the $p$ dependence with the geodesic ones. My bet is this will fail. The reason I think the hamiltonian ones are broken because the $r(s)$ curve doesn't seem to be effected by changing $p\phi(0)$.

15. Jun 6, 2017

### pervect

Staff Emeritus
If you're just interested in the answer, see for instance see https://www.fourmilab.ch/gravitation/orbits/.- Basically, two of the geodesic equatoins for the Shwarzschild metric can be manipulated into the form.

$$E(\tau)= \left(1-\frac{2M}{r} \right) \frac{dt}{d\tau} \quad \frac{dE}{d\tau}=0 \quad L(\tau) = r^2 \, \frac{d\phi}{d\tau} \quad \frac{dL}{d\tau}=0$$

If you expand dE and dL using the chain rule, you basically get two of the geodesic equations. You need one more equation to solve for the equations of motion, which can be taken as setting the magnitude of the 4 velocity to -1 (with a -+++ metric signature).

16. Jun 6, 2017

### Matter_Matters

Thanks for the information. I'm trying to write it in terms of hamiltonian's specifically to later use symplectic integration schemes. However, this is proving to be quite the task! Well for my coding capabilities anyways.

17. Jun 7, 2017

### Paul Colby

Code (Text):

q = {t[s],r[s],θ[s],ϕ[s]};
p = {pt[s], pr[s], pθ[s], pϕ[s]};
tt = 1-2 M /r[s];
rr = -1/tt;
θθ = -r[s]^2;
ϕϕ = -(r[s] Sin[θ[s]])^2;
metric = {{tt,0,0,0},{0,rr,0,0},{0,0,θθ,0},{0,0,0,ϕϕ}};
inversemetric = Simplify[Inverse[metric]];
ivs = {1.3,0,0,0.088}; ics = {0,6.5, π/2, 0};
hamiltonian = 1/2 (p . inversemetric . p - 1);
M=1;
hamiltonian /. Join[Table[q[[n]]-> ics[[n]],{n,1,4}],Table[p[[n]]-> ivs[[n]],{n,1,4}]]
0.720464

The way the Hamiltonian is written its value at the initial conditions should be 0, not 0.720464. It's not a physical particle so...?

18. Jun 7, 2017

### Paul Colby

For what it is worth, fixing the initial canonical momenta to be $p^\mu p_\mu = 1$ does little to change the straight into the hole issue. It is interesting that $p_\phi$ and $\theta$ are all constant which is the correct physics. I don't know what to make of $p_t$ being constant? This is clearly implied by the equations.

19. Jun 7, 2017

### Paul Colby

Decreasing the mass is amusing. Stable orbits result but they are highly elliptical orbits. I get a stable orbit for $M=0.01$ but it really zings around the mass at closest approach. It's possible straight into the hole is well deserved for your initial conditions?

20. Jun 7, 2017

### Paul Colby

For completeness here is a "working" code
Code (Text):
Clear["Global`*"]
q = {t[s],r[s],θ[s],ϕ[s]};
p = {pt[s], pr[s], pθ[s], pϕ[s]};
tt = 1-2 M /r[s];
rr = -1/tt;
θθ = -r[s]^2;
ϕϕ = -(r[s] Sin[θ[s]])^2;
metric = {{tt,0,0,0},{0,rr,0,0},{0,0,θθ,0},{0,0,0,ϕϕ}};
inversemetric = Simplify[Inverse[metric]];
ivs = {1,0,0,0.088}; ics = {0,6.5, π/2, 0};
hamiltonian = 1/2 (p . inversemetric . p - 1);
M=0.01
0.01
norm = Sqrt[p . inversemetric . p /.
Join[Table[q[[n]]-> ics[[n]],{n,1,4}],
Table[p[[n]]-> ivs[[n]],{n,1,4}]]]
1.00145
ivs = ivs / norm;
ivs
{0.998552,0.,0.,0.0878725}
hamiltonian /.
Join[Table[q[[n]]-> ics[[n]],{n,1,4}],
Table[p[[n]]-> ivs[[n]],{n,1,4}]]
1.06604*10^-16
qdot = Table[D[hamiltonian,p[[n]]],{n,1,4}];
pdot = Table[-D[hamiltonian,q[[n]]],{n,1,4}];
qequs = Table[D[q[[n]],s]==qdot[[n]],{n,1,4}];
qbcs = Table[(q[[n]]/. s-> 0) == ics[[n]],{n,1,4}];
pbcs = Table[(p[[n]]/. s-> 0) == ivs[[n]],{n,1,4}];
pequs = Table[D[p[[n]],s]==pdot[[n]],{n,1,4}];
eqs1 = Join[qequs,pequs,qbcs, pbcs];
time = {s,0,1000};
pbcs
{pt(0)==0.998552,pr(0)==0.,pθ(0)==0.,pϕ(0)==0.0878725}
sol = NDSolve[eqs1,Join[p,q],time,StartingStepSize->0.01];