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];

21. Jun 7, 2017

vanhees71

Let's start in clear formulae not some mediocre computer code. The most simple method is to use the quadratic form of the action principle, i.e.,
$$L=\frac{1}{2} g_{\mu \nu} \dot{x}^{\mu} \dot{x}^{\nu}.$$
It automatically leads to the parametrization with an affine parameter since
$$H=p_{\mu} \dot{x}^{\mu}-L=L$$
is conserved, i.e.,
$$g_{\mu \nu} \dot{x}^{\mu} \dot{x}^{\nu}=\text{const}.$$
A massive particle has a time-like trajectory, so that you can choose the overall scale of the parameter such that
$$\text{const}=1,$$
and then the affine parameter is the proper time $\tau$ of the particle.

Further you have
$$p_{\mu}=\frac{\partial L}{\partial \dot{x}^{\mu}}=g_{\mu \nu} \dot{x}^{\nu},$$
which implies that the Hamiltonian in terms of the canonical momenta (as is needed for the Hamilton canonical equations) is
$$H=\frac{1}{2} g^{\mu \nu} p_{\mu} p_{\nu}.$$
Now let's check the equations of motion

(a) Lagrange formalism

$$\dot{p}_{\mu}=\frac{\mathrm{d}}{\mathrm{d} \lambda} [g_{\mu \nu} \dot{x}^{\nu}]=g_{\mu \nu} \dot{x}^{\nu} + \partial_{\alpha} g_{\mu \nu} \dot{x}^{\alpha} \dot{x}^{\alpha} = \frac{\partial L}{\partial x^{\mu}}=\frac{1}{2} \partial_{\mu} g_{\alpha \nu} \dot{x}^{\alpha} \dot{x}^{\nu},$$
or in the more familiar form
$$g_{\mu \nu} \ddot{x}^{\nu}+\Gamma_{\mu \alpha \nu} \dot{x}^{\alpha} \dot{x}^{\nu},$$
where
$$\Gamma_{\mu \alpha \nu}=\frac{1}{2} (\partial_{\alpha} g_{\mu \nu} + \partial_{\nu} g_{\mu \alpha}-\partial_{\mu} g_{\alpha \nu})=g_{\mu \rho} {\Gamma^{\rho}}_{\alpha \nu},$$
with the usual Christoffel symbols for the Lorentzian spacetime manifold.

(b) Hamilton formalism

Here the equations of motion are given by the Hamilton canonical equations,
$$\dot{x}^{\mu}=\frac{\partial H}{\partial p_{\mu}}=g^{\mu \nu} p_{\nu}, \quad \dot{p}_{\mu}=-\frac{\partial H}{\partial x^{\mu}}=-\frac{1}{2} \partial_{\mu} g^{\alpha \nu} p_{\alpha} p_{\nu}.$$
From the first equation you get
$$p_{\mu} = g_{\mu \nu} \dot{x}^{\nu}.$$
$$\dot{p}_{\mu}=g_{\mu \nu} \ddot{x}^{\nu}+\partial_{\alpha} g_{\mu \nu} \dot{x}^{\mu} \dot{x}^{\alpha}=-\frac{1}{2} (\partial_{\mu} g^{\alpha \nu}) g_{\alpha \rho} g_{\nu \sigma} \dot{x}^{\rho} \dot{x}^{\sigma}=+\frac{1}{2} g^{\alpha \nu} \partial_{\mu} g_{\alpha \rho} g_{\nu \sigma} \dot{x}^{\rho} \dot{x}^{\sigma}=\frac{1}{2}\partial_{\mu} g_{\alpha \rho} \dot{x}^{\alpha} \dot{x}^{\rho} .$$
This is obviously the same equation we derived above within the Lagrange formalism. That's no surprise, because the Lagrange and Hamilton formalisms derive the equations of motion from the same action principle and thus lead the same equations of motion.

Last edited: Jun 7, 2017
22. Jun 7, 2017

Paul Colby

Totally true. I use this argument all the time in discussions with quality assurance about code bugs. Never seems to work.

23. Jun 7, 2017

pervect

Staff Emeritus

Possibly I'm confused, but something seems wrong here.

Let's take the flat space case when $g_{\mu\nu} = \eta_{\mu\nu}$ = diag(-1,1,1,1), so we have a -+++ metric signature.

We know that the relativistic Lagrangian for a free particle in the above flat space-time should be $L = -m\,c^2 \sqrt{1-\beta^2}$, where $\beta = v/c$. If we use geometric units so that c=1, we can write this as:

$L(t,\textbf{x},\frac{d \textbf{x}}{dt}) = -m \, \sqrt{ -g_{\mu\nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt}}$

where $\textbf{x}$ is a spatial vector, and the implied sum is only over the spatial terms. I am using the non-covariant formalism of the Lagrangian, but when I look up the covariant formalism of the Lagrangian in Goldstein, I still see square roots.

I don't see how to reconcile this with what you wrote. I've seen in the past Einbein forrmalisms that eliminate the square root from my version of the Lagrangian above by introducing another variable, but if that's what's going on, a bit more explanation would be quite helpful.

Another way of putting this -a change in the action $\Delta S = L \Delta t$ should be proportional to change in proper time $\Delta \tau$, so that extremizing the action of a free particle is the same as extremizing proper time, and in flat space time this means that L is proportional to $\frac{d\tau}{dt} = 1/\gamma = \sqrt{1-\beta^2}$. And we know that even in curved space-time, test particles move along worldlines that extremize proper time. But I don't quite see how that applies to what you (vanhees71) wrote.

Last edited: Jun 7, 2017
24. Jun 7, 2017

Staff: Mentor

I think you left out a $1$ under the square root.

Yes, you just do $g_{\mu \nu} U^\mu U^\nu$ under the square root where $U^\mu$ is now the 4-velocity vector, i.e., a 4-vector, not a spatial 3-vector. (And of course there will be no $1$ under the square root this time.)

25. Jun 7, 2017

pervect

Staff Emeritus
Good points - I edited my post to just strike out that portion of my response - and fix up a few other things that bothered me on re-reading. I still think there's a square root missing from vanhees71 version of the Lagrangian.

Last edited: Jun 7, 2017