Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

A EoM in Schwarzschild geometry: geodesic v Hamilton formalism

  1. Jun 5, 2017 #1
    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. jcsd
  3. Jun 5, 2017 #2

    Paul Colby

    User Avatar
    Gold Member

    Do the two approaches produce the same equations of motion? This could be directly checked in Mathematica.
     
  4. Jun 5, 2017 #3

    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!
     
  5. Jun 5, 2017 #4

    Paul Colby

    User Avatar
    Gold Member

    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}}##
     
  6. Jun 5, 2017 #5

    pervect

    User Avatar
    Staff Emeritus
    Science Advisor

    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>>.
     
  7. Jun 5, 2017 #6

    Paul Colby

    User Avatar
    Gold Member

    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.
     
  8. Jun 5, 2017 #7

    pervect

    User Avatar
    Staff Emeritus
    Science Advisor

    ##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.
     
  9. Jun 5, 2017 #8
    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.
     
  10. Jun 5, 2017 #9

    Paul Colby

    User Avatar
    Gold Member

    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.
     
  11. Jun 6, 2017 #10
    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.
     
  12. Jun 6, 2017 #11

    Paul Colby

    User Avatar
    Gold Member

    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
  13. Jun 6, 2017 #12

    Paul Colby

    User Avatar
    Gold Member

    Okay, I get a stiff system warning at ##s=10.8## which isn't too surprising.
     
  14. Jun 6, 2017 #13
    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.
     
  15. Jun 6, 2017 #14

    Paul Colby

    User Avatar
    Gold Member

    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)##.
     
  16. Jun 6, 2017 #15

    pervect

    User Avatar
    Staff Emeritus
    Science Advisor

    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).
     
  17. Jun 6, 2017 #16
    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.
     
  18. Jun 7, 2017 #17

    Paul Colby

    User Avatar
    Gold Member

    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...?
     
  19. Jun 7, 2017 #18

    Paul Colby

    User Avatar
    Gold Member

    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.
     
  20. Jun 7, 2017 #19

    Paul Colby

    User Avatar
    Gold Member

    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?
     
  21. Jun 7, 2017 #20

    Paul Colby

    User Avatar
    Gold Member

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

     
     
  22. Jun 7, 2017 #21

    vanhees71

    User Avatar
    Science Advisor
    Gold Member
    2017 Award

    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}.$$
    Taking the derivative leads to
    $$\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
  23. Jun 7, 2017 #22

    Paul Colby

    User Avatar
    Gold Member

    Totally true. I use this argument all the time in discussions with quality assurance about code bugs. Never seems to work.
     
  24. Jun 7, 2017 #23

    pervect

    User Avatar
    Staff Emeritus
    Science Advisor


    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
  25. Jun 7, 2017 #24

    PeterDonis

    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.)
     
  26. Jun 7, 2017 #25

    pervect

    User Avatar
    Staff Emeritus
    Science Advisor

    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
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted