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

A How to Integrate the geodesic equations numerically?

  1. Mar 7, 2017 #1
    Hello there, I've been considering the geodesic equations of motion for a test particle in Schwarzschild geometry for some time now. Similar to what we can do with the Kepler problem I would like to be able to numerically integrate the equations of motion. I'm quite interested to see how coordinate time and proper time change over a simulation and also the trajectory itself. But the problem is.. I have no idea how to do that. Any suggestions?
  2. jcsd
  3. Mar 7, 2017 #2


    User Avatar
    Science Advisor

  4. Mar 7, 2017 #3


    User Avatar
    Staff Emeritus
    Science Advisor

    The Schwarzschild geodesic has a couple of conserved quantites that mean you don't have to numerically integrate the equations, you can solve them analytically by taking advantage of the existence of these constants of motion.

    Sean Caroll discusses this in his lecture notes on GR, though it's not an easy read. The applicable concept is called "Killing Vectors", a quote from the relevant section (found at https://www.preposterousuniverse.com/grnotes/) is:

    For another online reference, http://www.fourmilab.ch/gravitation/orbits/ might or might not be easier to follow - it is taken pretty much from MTW's treatment in "Gravitation", which would be a better source than the above (being more detailed and of better provenance) - if you can get a hold of it. This approach is a bit more physical, they give the conserved quantites names, like "energy at infinity" and "angular momentum". The resulting equations are rather similar to how we might solve for motion in the Newtonian case by taking advantage of energy and angular momentum conservation.

    It might help to review the "effective potential" method for the discussion in the second link. Orbits occur in a plane, so they're basically two dimensional. So one can introduce unit vectors ##\hat{r}## and ##\hat{\phi}## in the Newtonian case to describe the orbital plane, the former being in the radial direction, the other being at right angles to the radial vector in the orbital plane.

    By the conservation of angular momentum, one can write a relationship between the the rate of change of ##\phi## and the radial distance r. This gives the non-radial component of the velocity. Knowing the radial and transverse-to-radial components of the velocity, we can get the total kinetic energy, and using the fact that energy is conserved, we can solve for the radial and transverse components of the velocity at any point.

    The "effective potential" method is basically a matter of interpreting this solution as the motion of a 1 dimensional particle in a force field that gives the proper solution for the radial component of the motion. Many texts should discuss the effective potential method, including Goldstein's "Classical Mechanics".
  5. Mar 7, 2017 #4
    [To the OP]
    Yes, this is most definitely the way to go. There are some, er, "fascinating" technical issues you will need to deal with, but it's far better than just blasting away at the second order ODEs in the geodesic equation (particularly if you like to cross horizons!).

    Oh, couldn't resist. Here is my take on the fourmilab animations, with a spinning black hole (it has one solar mass, and is orbited by a test particle similar to Mercury). Also does Newtonian orbit for "comparison".
  6. Mar 10, 2017 #5
    Absolutely fantastic responses both very helpful and also a lot I didn't think about. Being a "trained mathematician" I was always satisfied to leave the geodesic equations as they are and have no other concerns i.e. I was completely content with just having a set of DE's and that was that. However, now that my needs require more involvement I find myself almost lost in the technicalities of things. Still though, absolutely fascinating to study. I do find it interesting that such topics are not heavily discussed in the literature, after all, nearly every student would know a member of the Runge-Kutta family. So extending this should be a topic of interest. Don't you agree?
  7. Mar 10, 2017 #6
    Since you are a trained mathematician (I most definitely am NOT), I will try to be brief ;) The seminal paper on geodesics IMO is this one by Wilkins. The first order equations for the Kerr metric are all there in plain sight in section II.

    One other topic I would bring up is the Symplectic Integrator, which is frequently used over RK4 for long term stability.

    And finally, a plug for my proudest "contribution" to the field, in case you like to mess with working code.
  8. Mar 10, 2017 #7


    User Avatar
    Staff Emeritus
    Science Advisor

    The end result is that if you consider the motion in an equatorial orbital plane (t,r, ##\phi##), the solution for writing the geodesic equation is three functions ##t(\tau), r(\tau), \phi(\tau)##.

    We can define ## E = (1-2m/r) dt / d\tau## and ##L = r^2 (d\phi/d\tau)## and re-write two of these equations as ##dE/d\tau=0, dL/d\tau=0##, using the chain rule to expand them. The last equation comes from the Schwarzschild metric itself, which has various forms. One is:

    $$c^2 \,{d \tau}^{2} = \left(1 - 2m/r \right) c^2 \,dt^2 - \left(1- 2m/r \right)^{-1} \,dr^2 - r^2 \left(d\theta^2 + \sin^2\theta \, d\phi^2\right)$$

    but the fourmilab reference assumes c=1, for instance, and by assuming an equatorial orbit we have previously set ##\theta=0##. Making these assumptions and dividing both sides by ##d\tau^2## should basically give the last equation.

    I think I posted something more detailed along these lines a long time ago, but I'm not finding it in search.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted