# Schwarzschild Orbits in Cartesian coordinates

## Is Carl going to find the Schwarzschild orbits in Cartesian DE form?

• Total voters
14
Homework Helper
My Java applet gravity simulator http://www.gaugegravity.com/testapplet/SweetGravity.html [Broken]
draws beautiful orbits, however the GR simulation is very badly broken as one can tell when comparing it with Newton at long distances. The source code is at:
http://www.gaugegravity.com/testapplet/SwGrav_Top.java [Broken]

Test bodies are required to all lie in a plane. The Newtonian physics is produced by writing the x and y accelerations in terms of the phase space parameters of the test body, that is, in terms of its (x,y) position and its (Vx,Vy) velocity. With Newton, the acceleration has no dependence on velocity of course.

The GR "physics" is produced by a complicated calculation that follows the artificial potential that is discussed with respect to the Schwarzschild orbits in MTW and in various places all over the web. These methods use energy and angular momentum as constants of the motion to describe the orbits but they are less than optimal when used in an integration form.

I've seen a computer method for computing GR orbits that should work better, but I am hesitating to implement it because it is based on polar coordinates. This means that the calculation has to use trig functions to convert back and forth and these are slow.

So I'm wondering how hard it is to put the Schwarzschild orbits into phase space form in Cartesian coordinates. Here's the basic plan:

(1) Write the Schwarzschild metric in Cartesian coordinates.

(2) Write the proper length of a path as an integral over coordinate time.

(3) Vary the path and use the Euler-Lagarange equation to determine a pair of 2nd order differential equations that the orbits solve.

(4) Find $$d^2x/dt^2$$ and $$d^2y/dt^2$$ from the two 2nd order DEs.

Now for the poll: Do you think that this can be done reasonably easily? Part of the reason I am working on it is that I want to get a better idea on how these gravity theories differ one from another. I will post calculations as I make them, if convenient.

Carl

Last edited by a moderator:

## Answers and Replies

Related Special and General Relativity News on Phys.org
Homework Helper
Step (1), write the Schwarzschild metric in Cartesian coordinates.

$$x = r \cos(\phi)$$
$$y = r \sin(\phi)$$

$$dx = \cos(\phi) dr - r\sin(\phi) d\phi$$
$$dy = \sin(\phi) dr + r\cos(\phi) d\phi$$

$$dr = \cos(\phi) dx + \sin(\phi) dy$$
$$rd\phi = -\sin(\phi) dx + \cos(\phi) dy$$

$$dr^2 = \cos^2(\phi) dx^2 + \sin^2(\phi) dy^2 + \sin(2\phi) dx\;dy$$
$$r^2d\phi^2 = \sin^2(\phi) dx^2 + \cos^2(\phi) dy^2 - \sin(2\phi) dx\;dy$$

$$dr^2 = (x^2dx^2 + y^2dy^2 + 2xy dx\;dy)/r^2$$
$$r^2d\phi^2 = (y^2dx^2 + x^2dy^2 - 2xy dx\;dy)/r^2$$

$$ds^2 = -\frac{(r-2)}{r}dt^2 + \frac{r}{r-2}dr^2 + r^2d\phi^2$$

$$ds^2 = -\frac{(r-2)}{r}dt^2 + \frac{1}{r(r-2)}(x^2dx^2 + y^2dy^2 + 2xy dx\;dy) + \frac{1}{r^2}(y^2dx^2 + x^2dy^2 - 2xy dx\;dy)$$

$$ds^2 = -\frac{(r-2)}{r}dt^2 + \frac{1}{r^2(r-2)}((rx^2+ ry^2 -2y^2)dx^2 + (ry^2+rx^2-2x^2)dy^2 + (2rxy -2rxy +4rxy) dx\;dy)$$

$$ds^2 = -\frac{(r-2)}{r}dt^2 + \frac{1}{r^2(r-2)}((r^3 -2y^2)dx^2 + (r^3-2x^2)dy^2 + (4rxy) dx\;dy)$$

$$ds^2 = -\frac{(r-2)}{r}dt^2 + \frac{(r^3-2y^2)}{(r^3-2r^2)}dx^2 + \frac{(r^3-2x^2)}{(r^3-2r^2)}dy^2 +\frac{4rxy}{(r^3-2r^2)}dx\;dy$$

Last edited:
Jorrie
Gold Member
CarlB said:
...
Now for the poll: Do you think that this can be done reasonably easily? Part of the reason I am working on it is that I want to get a better idea on how these gravity theories differ one from another. I will post calculations as I make them, if convenient.
Carl
Interesting excercise Carl, but how do you think having the orbital equation in Cartesian form (as against polar form) will help you get a "better idea on how these gravity theories differ one from another."?

pervect
Staff Emeritus
There isn't any such thing as "Cartesian coordinates" for a black hole. Cartesian coordinates work only for flat space-times. The space-time around a black hole is not flat. In particular, r is not a radial coordinate.

The closest you can come to "cartesian" coordiantes are isotropic coordinates.

Howver, from experience, I can tell you it will be far far easier to get your program working for Schwarzschild coordinates than isotropic cooridinates.

As I mentioned before, you can just use the geodesic equations to express your set of equations as a second order system that look pretty much like the second order set you have already for Newton's equations.

This is the "geodesic" equation

http://en.wikipedia.org/wiki/Solving_the_geodesic_equations

(Wiki somewhat confusingly uses "t" instead of "tau" for the variables.

The necessary Christoffel symbols are simplest for the schwarzschild metric, they are

Code:
                        r                 m
CC   [r r] = - -----------
(r - 2 m) r

theta
CC       [theta r] = 1/r

phi
CC     [phi r] = 1/r

t               m
CC   [t r] = -----------
(r - 2 m) r

theta
CC       [r theta] = 1/r

r
CC   [theta theta] = -r + 2 m

phi                cos(theta)
CC     [phi theta] = ----------
sin(theta)

phi
CC     [r phi] = 1/r

phi                cos(theta)
CC     [theta phi] = ----------
sin(theta)

r                                   2
CC   [phi phi] = -(r - 2 m) sin(theta)

theta
CC       [phi phi] = -sin(theta) cos(theta)

t               m
CC   [r t] = -----------
(r - 2 m) r

r          (r - 2 m) m
CC   [t t] = -----------
3
r
As I also mentioned, you should be able to obtain the same set of geodesic equations from the conserved quantites, as well, by replacing "conserve quantity = constant" with "d/dtau (conserved quantity) = 0"

Last edited:
pervect
Staff Emeritus
Let's try to make things really, really, really simple.

theta will always be 90 degrees. Therfore d theta / dtau will be zero.

We then have for our variables

r(tau), t(tau), phi(tau).

The second order geodesic equations of motion (the geodesic equations) will therfore be

$$\frac{d^2 r}{d tau^2} + \Gamma^r{}_{rr} \left( \frac{dr}{dtau} \right) ^2 + \Gamma^r{}_{tt} \left( \frac{dt}{d tau} \right)^2 + \Gamma^r{}_{\phi\phi} \left( \frac{d \phi}{d tau} \right)^2 = 0$$

$$\frac{d^2 t}{d tau^2} + \left( \Gamma^t{}_{tr} + \Gamma^t{}_{rt} \right) \left( \frac{dr}{d tau} \right) \left( \frac{dt}{dtau} \right) = 0$$

$$\frac{d^2 \phi}{d tau^2} + \left( \Gamma^{\phi}{}_{r\phi} + \Gamma^{\phi}{}_{\phi r} \right) \left( \frac{dr}{d tau} \right) \left( \frac{d \phi} {d tau} \right) = 0$$

See the previous post for the values for the Christoffel symbols.

Note for instance that since r^2 ${d \phi}/ {d tau}$ = constant, when differentiated, gives by the chain rule

$$2 r \frac{dr}{d\tau} \frac{d \phi}{d tau} + r^2 \frac{d^2 \phi} {d tau^2} = 0$$

and that this is equivalent to the last equation above when you substitute in the correct value for the Christoffel symbol.

The conseved energy when differentiated will give one of the other geodesic equations (the one for time).

The last geodesic equation is a consequence of the others plus the equation for the metric.

Last edited:
Jorrie
Gold Member
Cartesian coordinates

pervect said:
There isn't any such thing as "Cartesian coordinates" for a black hole. Cartesian coordinates work only for flat space-times. The space-time around a black hole is not flat. In particular, r is not a radial coordinate.

The closest you can come to "Cartesian" coordinates are isotropic coordinates.
Are you not a bit harsh on CarlB? I understood that he simply wants to express Schwarzschild coordinates in terms of x,y,z rather than in terms of r,phi,theta. Nothing fundamentally wrong with that, I think! Maybe his (Cartesian) choice of words was wrong, but let’s not get too semantic!

Nevertheless, I agree that we should rather stick to the tried and tested procedure of treating black holes.

Homework Helper
pervect, one of the things I'm trying to get away from is differential equations written in proper time. It makes for a more complicated simulation and also makes it a lot easier for the simulation to fail in various ways. So I'm continuing with the process of writing the equations in Cartesian coordinates. I will allow you to continue calling them isotropic coordinates if that's what you want to do.

Carl

pervect
Staff Emeritus
Carl (and Jorrie, too) I think you've totally missed the point.

The 'r' in the Scwarzschild coordinate system is not the same 'r' in the isotropic coordinate system.

In fact, one can see that r = r'(1+2m/r')^2 from the webpage, where r' is the isotropic radial coordinate and r is the Schwarzschild r coordinate.

Neither of these 'r' coordinates can be interpreted as 'radial distance'. The isotropic coordinates, though, define a set of coordinates where light cones will apear to be circles.

Thus assuming that r^2 = x^2 + y^2 does not make any physical sense. This is assumed implicitly when one writes x = r cos(phi) y = r sin(phi), one assumes that x^2 + y^2 = r^2. One can think of it as a purely formal mathematical expression, I suppose.

As such a formal exercise, unless Carl has made an algebraic error he's got a new and totally non-standard metric that doesn't even have a name in terms of new and toatally non-standard coordiantes that don't have names either (except what Carl has christened them). They definitely shouldn't be confused with Cartesian coordinates, which is what Carl has chosen to call his variables.

One can easily see that Carl's metric is not the same as the isotropic metric either by comparing the line elements.

Aside from the fact that it is non-standard, and thus will be hard to check solutions for, even deriving the correct solutions is going to be a giant pain. Re-writing a simple metric in a more complicated form is not going to get anywyere as far as generating a solution for the metric.

If the problem is with the factor of 'tau' in the standard equations, just eliminate that tau. This shouldn't be too hard since

(1-2m/r) dt/dtau = E

The most annoying thing is going to be this

If we write r as a function of t : r = rr(t) then

r(tau) = rr(t(tau))

and dr/dtau = drr/dt dt/dtau

but by the chain rule

$$\frac{d^2 r} {d tau^2} = \frac{d^2 rr}{d t^2} \left( \frac{dt}{d tau} \right) ^2 + \frac{dr}{dt} \frac{d^2 t}{d tau^2}$$

Still, a little bit of algebra should eliminate tau as a variable if this is really what is desired.

Solving for phi shouldn't be terriby difficult, since

r^2 dphi/dtau =L

therfore

r^2 dphi/dt * dt/dtau = L

And we've already noted that dt/dtau is easy to eliminate

pervect
Staff Emeritus
Unless I've made a mistake (quite possible, even with computer assistance) the equations with proper time eliminated for the Schwarzschld metric are just

r(t), phi(t)

d phi/dt = K*(1-2*m/r) / r^2

K is a constant of motion, equal to "L/E" in the old notation. Knowing dphi/dt and r at one point will determine the value of K

The second order equation for r(t) is

$$\frac{d^2 r}{dt^2} = \frac{2\,m}{r^2 \left( 1 - \frac{2m}{r} \right) } \left( \frac{dr}{dt} \right) ^2 - \left( 1-{\frac {2 m}{r}} \right) \frac{m}{r^2}+ \left( 1-{\frac {2m}{r}} \right) ^{3}\frac{{K}^{2}}{r^3}$$

It looks sensible at first glance, dr/dt is forced to zero, however I have not tested it.

Compare this to the newtonian equations in polar coordinates (r, phi)

r^2 * d phi/dt = K

where K = L/m

d^2 r / dt^2 = -m/r^2 + r (dphi/dt)^2 = -m/r^2 + K^2 / r^3

we can also see that the Schwarzschild equations reduce to the Newtonian equations when dr/dt << 1 and r >> 2m as they should, another sign that they are probably right.

Last edited:
pervect
Staff Emeritus
The Christoffel symbols for the Painleve metric are:

Code:
                                   1/2      1/2
T          2    (M/r)    M
CC   [T T] = ---------------
2
r

r          (r - 2 M) M
CC   [T T] = -----------
3
r

T           M
CC   [r T] = ----
2
r

1/2      1/2
r            2    (M/r)    M
CC   [r T] = - ---------------
2
r

T           M
CC   [T r] = ----
2
r

1/2      1/2
r            2    (M/r)    M
CC   [T r] = - ---------------
2
r

1/2
T             2    M
CC   [r r] = -------------
1/2  2
2 (M/r)    r

r             M
CC   [r r] = - ----
2
r

theta
CC       [theta r] = 1/r

phi
CC     [phi r] = 1/r

theta
CC       [r theta] = 1/r

T                    1/2      1/2
CC   [theta theta] = -2    (M/r)    r

r
CC   [theta theta] = -r + 2 M

phi                cos(theta)
CC     [phi theta] = ----------
sin(theta)

phi
CC     [r phi] = 1/r

phi                cos(theta)
CC     [theta phi] = ----------
sin(theta)

T                1/2      1/2             2
CC   [phi phi] = -2    (M/r)    r sin(theta)

r                                   2
CC   [phi phi] = -(r - 2 M) sin(theta)

theta
CC       [phi phi] = -sin(theta) cos(theta)
a similar analysis using the geodesic equations and eliminating tau should be possible.

Homework Helper
pervect said:
Carl (and Jorrie, too) I think you've totally missed the point.

The 'r' in the Scwarzschild coordinate system is not the same 'r' in the isotropic coordinate system.

In fact, one can see that r = r'(1+2m/r')^2 from the webpage, where r' is the isotropic radial coordinate and r is the Schwarzschild r coordinate.

Neither of these 'r' coordinates can be interpreted as 'radial distance'. The isotropic coordinates, though, define a set of coordinates where light cones will apear to be circles.

Thus assuming that r^2 = x^2 + y^2 does not make any physical sense. This is assumed implicitly when one writes x = r cos(phi) y = r sin(phi), one assumes that x^2 + y^2 = r^2. One can think of it as a purely formal mathematical expression, I suppose.

As such a formal exercise, unless Carl has made an algebraic error he's got a new and totally non-standard metric that doesn't even have a name in terms of new and toatally non-standard coordiantes that don't have names either (except what Carl has christened them). They definitely shouldn't be confused with Cartesian coordinates, which is what Carl has chosen to call his variables.

Aside from the fact that it is non-standard, and thus will be hard to check solutions for, even deriving the correct solutions is going to be a giant pain. Re-writing a simple metric in a more complicated form is not going to get anywyere as far as generating a solution for the metric.

It is my intention to work only in the Schwarzschild metric, but with the polar coordinates replaced by Cartesian coordinates. I realize that you think this is unphysical, but here in the United States computers are built with screens that have rectangular arrays of pixels, and so any computation I make must be eventually so converted. Now that you've explained what they are, I don't mean to use "isotropic" coordinates at all. It's just extra confusion you've brought in, not me. I thought it was you calling my coordinates "isotropic". Recall, it was I who "allowed" you to call them that if you like.

Nor am I the only person in the world who calls these things I'm using "Cartesian". The United States is a very primitive place, filled with horses, saloons and cowboys. Other obscure, amateurish authors here use the same confusing and incorrect terminology. For example see:

Numerical Relativity for Inspiraling Binaries in Co-Rotating Coordinates: Test Bed for Lapse and Shift Equations
Kip S. Thorne, 1998
http://arxiv.org/abs/gr-qc/9808024

By the way, I should mention that while I only took one graduate level General Relativity class at the University of California, Irvine, I was the top student, the only student in the class receiving an "A+". I don't think that taking one class makes you an expert in a subject, but you can be sure that I understand the principles of general relativity at least a little better than some.

 But I see that it is not only on this side of the Atlantic that one finds the use of Cartesian coordinates in describing black holes:

Shortcut Method of Solution of Geodesic Equations for Schwarzschild Black Hole
Jean-Alain Marck,
D´epartement d’Astrophysique Relativiste et de Cosmologie
Observatoire de Paris - Section de Meudon
F-92195 Meudon Cedex, France

Classical and Quantum Gravity 13 (1996) 393-402

$$r^2 = x^2 + y^2 + z^2\;\;\;\;\;\;\;(6)$$

http://arxiv.org/abs/gr-qc/9505010
[/edit]

Carl

Last edited:
Homework Helper
pervect said:
Unless I've made a mistake (quite possible, even with computer assistance) the equations with proper time eliminated for the Schwarzschld metric are just:
This is exactly what I need and much thanks. To convert it into Cartesian form, I need only split the radial acceleration into x and y accelerations. I should have a new simulation running later today, but Seattle is sooooooo hot (about 90 F) that everyone has left for the beach and very little work is getting done. Maybe I'll wait until it begins to cool down.

I was in mathematics grad school at the University of Washington 10 years ago and got a student priced copy of Mathematica. I played around with it, but never worked at it long enough to get around to messing with the symbolic calculus. Since then, most of my work has been in Clifford algebra and for that I wrote Java programs. But after this fiasco I'm convinced that I need to get some sort of a symbolic math package.

Carl

Last edited:
pervect
Staff Emeritus
Symbolic math really rocks. Couldn't live without it.

I've also added horizon crossing solutions to the Painleve thread (but only in terms of conseved quantites) I haven't eliminated dt/dtau, it appears to be a royal pain to do so. dt/dtau will range in values from 1 to E for the horizon crossing solutions.

I note that non-horizon crossing alternate solutions exist, but I haven't explored them very far.

I've also computed d^2 r / dtau^2 only in the special case where dr/dtau = 0.

Last edited:
Homework Helper
pervect said:
$$\frac{d^2 r}{dt^2} = \frac{2\,m}{r^2 \left( 1 - \frac{2m}{r} \right) } \left( \frac{dr}{dt} \right) ^2 - \left( 1-{\frac {2 m}{r}} \right) \frac{m}{r^2}+ \left( 1-{\frac {2m}{r}} \right) ^{3}\frac{{K}^{2}}{r^3}$$

Unfortunately, I couldn't get it to work. This could always be a programming error of course. It acts like there is too much energy in the angular kinetic energy term, that is the K^2 term. In fact, the simulation doesn't match Newton at large r because of this term. To get it to work, I have to divide by an extra factor of r.

Take a look at the units for large r, in the first and third terms of the above and see if there is a transcription error, most likely in the last term of the above, but at least conceivably in the formula for $$K = \dot{\phi} r^2 / (1-2m/r)$$ which can be used to simplify the third term.

Having a differential equation in r and a constant of the motion in phi is unsatisfying to me. I realize that I can put the Newtonian equation into the same form, but I am one of those crazy guys who thinks that the equations of motion are closer to the physics than the constants of the motion. Yes, I believe that the constants of the motion are nothing more than convenient methods we use for simplifying the original physical problem, which to me are always defined in terms of positions and velocities. So it looks like I will continue my process of deriving the Newton style equations of motion for the Schwarzschild metric.

Carl

Homework Helper
Step (2), write the proper length of a path as an integral over coordinate time.

Here's where I begin:

CarlB said:
$$ds^2 = \frac{(r-2)}{r}dt^2 - \frac{(r^3-2y^2)}{(r^3-2r^2)}dx^2 - \frac{(r^3-2x^2)}{(r^3-2r^2)}dy^2 -\frac{4rxy}{(r^3-2r^2)}dx\;dy$$
Note that in the above I've swapped over to the other coast metric in order to make s time like instead of spacelike. Sorry for starting the problem with the spacelike signature, it's a habit from Euclidean relativity where s is interpreted as a little curled up spatial dimension. Dividing by (dt)^2 gives:

$$\left(\frac{ds}{dt}\right)^2 = \frac{(r-2)}{r} - \frac{(r^3-2y^2)}{(r^3-2r^2)}\dot{x}^2 - \frac{(r^3-2x^2)}{(r^3-2r^2)}\dot{y}^2 -\frac{4rxy}{(r^3-2r^2)}\dot{x}\dot{y}$$

where the dots indicate derivatives with respect to coordinate time t rather than the more usual proper time. Thus the proper time along the path is:

$$S = \int_{t_1}^{t_2} \frac{ds}{dt}\; dt$$

$$= \int_{t_1}^{t_2} \sqrt{\frac{(r-2)}{r} - \frac{(r^3-2y^2)}{(r^3-2r^2)}\dot{x}^2 - \frac{(r^3-2x^2)}{(r^3-2r^2)}\dot{y}^2 -\frac{4rxy}{(r^3-2r^2)}\dot{x}\dot{y} }\;\; dt$$

Last edited:
Homework Helper
The next step is to vary the path by the Euler-Lagrange equations. Discretion being the better part of valor, I will wait until I have a symbolic mathematics engine before I attempt that slope. Basically, I need to calculate the partial derivatives with respect to the coordinates and velocities. The equations of motion are then the time derivatives of the partials with respect to velocities set equal to the partials with respect to positions. That is, for example,

$$\frac{d}{dt} \left(\frac{\partial ds/dt}{\partial x}\right) = \frac{\partial ds/dt}{\partial \dot{x}}$$

Carl

Okay, now I've got Maxima running on my laptop. Pretty cool. And cheap.

Last edited:
pervect
Staff Emeritus
If you use Maxima to calculate the Christoffel symbols, watch out for the reversed ordering of the indices. I fell afoul of that issue in GrTensor. Maxima does it to, I believe (I've only used that briefly).

If you have mathematica I think grtensor version M will work with it. Grtensor is better than Maxima in my experience, and Maxima is much better than nothing at all.

I've recalculated the formula for the Schwarzschild metric, and gotten a different answer. The main difference is that I removed some manual copying steps, collecting the values of the Christoffel symbols automatically instead of transcribing them by hand. The new formula seems to conserve energy and momentum in the numerical test case I ran, though some variance can be noted near the horizon. This variance is improved by tightning up the tolearances of the integration routine, though, so I think it's a numerical issue.

The new formula are:

$$\frac{d^2 r}{d t^2} = \frac {3 m{{\it vr}}^{2}}{ \left( r-2\,m \right) r} + \left( r-2\,m \right) \left( {{\it vphi}}^{2}-{\frac {m}{{r}^{3}}} \right)$$

$$\frac{d^2 phi}{d t^2} = -\frac {2 {\it vr}\,{\it vphi}\, \left( r -3\,m \right) }{ \left( r-2\,m \right) r}$$

vr = dr/dt and vphi = dphi/dt

Last edited:
pervect
Staff Emeritus
Here's what the solution looks like at the code level, minus some manual simplifications of the formula.

manual simplifications are checked by computing simplify(manual - auto) and making sure it is zero.

CC(up,dn,dn) calculates the Christoffel symbols - it's a macro that uses grtensor's built-in routine and changes the order to the standard order.

>grcalc(CC(up,dn,dn));
> grdisplay(_);

> tmp := grarray(CC(up,dn,dn));
> vlist := [rdot,0,phidot,tdot];

> eq1 := 0;
> eq2 := 0;
> eq3 := 0;
> eq4 := 0;

> # sum geodesic equations
> for j from 1 by 1 to 4 do
> for k from 1 by 1 to 4 do
> eq1 := eq1 + tmp[1,j,k]*vlist[j]*vlist[k];
> eq2 := eq2 + tmp[2,j,k]*vlist[j]*vlist[k];
> eq3 := eq3 + tmp[3,j,k]*vlist[j]*vlist[k];
> eq4 := eq4 + tmp[4,j,k]*vlist[j]*vlist[k];
> end do;
> end do;
>
> eq1 := subs(sin(theta)=1,eq1);
> eq2 := subs(cos(theta)=0,eq2);
> eq3 := subs(sin(theta)=1,eq3);
> eq4 := subs(sin(theta)=1,eq4);

> # we have t(tau), r(tau), phi(tau)
> # tdot = dt/dtau, rdot = dr/dtau, phidot = dphi/dtau
> # we want vr(t) = dr/dt, vphi(t) = dphi/dt
> # vr * tdot = rdot; vphi * tdot = phidot
> # rdotdot = ar*tdot^2+vr*tdotdot
> # but tdotdot = -eq4, thus

> # equation for r, ar is d^2 r / dt^2
> eq1a := ar*tdot^2 -vr*eq4 + eq1;
> eq1b := subs(rdot = vr*tdot, phidot = vphi*tdot, eq1a);
> eq1c:= simplify(solve(eq1b,ar));
> collect(eq1c,vr);

> # phi equation
> eq3a := aphi*tdot^2 -vphi*eq4 + eq3;
> eq3b := subs(rdot = vr*tdot, phidot = vphi*tdot, eq3a);
> eq3c := simplify(solve(eq3b,aphi));
> collect(eq3c,vr);

pervect
Staff Emeritus
CarlB said:
The next step is to vary the path by the Euler-Lagrange equations. Discretion being the better part of valor, I will wait until I have a symbolic mathematics engine before I attempt that slope.
The standard textbook approach is to write the Lagrangian as

$$L = g_{\mu \nu} \frac{d x^\mu}{d \lambda} \frac{d x^\nu}{d \lambda} d \lambda$$

(Summation is implied).

Applyin the Euler-Lagrange equations to this Lagrangian will yield the geodesic equations.

Of course this approach yields equations parameterized in terms of proper time - i.e. the standard geodesic equations.

AFAIK the easiest way to eliminate proper time is to use the approach I outlined earlier. It is not immediately obvious why one doesn't need to solve for "tdot", but it seems to factor itself out automatically, so one need not solve for it.

pervect
Staff Emeritus
Finally, the painleve metric equations with proper time eliminated. They also seemed to conserve energy and angular momentum in the one test case I ran

$$\frac{d^2 r}{d t^2} = \sqrt {{\frac {2 M}{r}}}{\it vr}\, \left( {\frac {{{\it vr }}^{2}}{2 r}}+3\,{\frac {M}{{r}^{2}}}-r{{\it vphi}}^{2} \right) +3\,{ \frac {M{{\it vr}}^{2}}{{r}^{2}}}+ \left( r-2\,M \right) \left( {{ \it vphi}}^{2}-{\frac {M}{{r}^{3}}} \right)$$

$$\frac{d^2 phi}{d t^2} = {\it vphi}\, \left( \sqrt {{\frac {2 M}{r}}} \left( {\frac {{{\it vr}}^{2}}{2 r}}-r{{\it vphi}}^{2}+{\frac {M}{{r}^{2}}} \right) -2\,{\frac {{\it vr}\, \left( r-M \right) }{{r}^{2}}} \right)$$

One can see that while the Painleve equations are simple in the newtonian limit when vr=vphi is almost zero they are not as simple as the Schwarzschild equations in the fully relativistic case.

This could be anticipated by the "count the number of Christoffel symbols" method.

Last edited:
pervect
Staff Emeritus
Finally, since I have the automaton working, here's the sort of result the above procedure gets when applied to the psueudo-cartesian metric you suggested earlier.

d^2x/dt^2 :=
-(-6*x^3*vy^2*y^4+6*x^5*vy^2*y^2+4*vr*y^7*vy-8*x^5-6*y*vx*vy*r*x^6-4*x
^3*vx^2*y^6-8*x^3*y^2+2*vr*r*y^3*vy*x^4+4*x*vy^2*z^6*y^2+2*vr*r*x^5*vx
*y^2-6*vr*r*x^3*vx*z^4-8*x^3*z^2-8*x^5*y^2+2*x^3*y^4-16*x^5*z^2+2*x^9*
vy^2-6*x^7-12*x^3*y^2*z^2+8*x^7*vy^2+4*x^7*vy^2*y^2+8*x^7*vy^2*z^2+4*v
r*x^7*vx-4*x*vx^2*y^8-4*x^3*vx^2*y^4+4*x^5*vx^2*z^2+2*x^3*vx^2*z^4+4*x
^7*vx^2*y^2+4*x^5*vx^2*y^4-4*x^3*vy^2*y^6-2*x*vy^2*y^8+6*x^5*vy^2*z^2-
6*x^3*vy^2*z^4+12*x^5*vy^2*z^4+8*x^3*vy^2*z^6+12*y^3*vx*vy*x^4-8*y^5*v
x*vy*x^2+12*vr*x^5*vx*y^2+8*vr*x^5*vx*z^2+12*vr*x^3*vx*y^4+4*vr*x^3*vx
*z^4+12*vr*y^3*vy*x^4+12*vr*y^5*vy*x^2+8*vr*y^5*vy*z^2+4*vr*y^3*vy*z^4
-4*x*vx^2*y^6-12*x*vx^2*y^6*z^2+6*x^3*vx^2*z^2*y^2+2*x^5*vx^2*y^2+4*x^
5*vx^2*y^2*z^2-4*x*vy^2*y^6-4*x*vy^2*y^6*z^2-12*x^3*vy^2*y^2*z^2-12*x*
vy^2*y^4*z^2+12*x^5*vy^2*y^2*z^2+12*x^3*vy^2*z^4*y^2-4*x*vy^2*z^6-12*x
*vy^2*z^4*y^2+12*y^3*vx*vy*x^2*z^2+4*y^6*x-14*x^3*z^4-4*x*z^6+2*x^7*vx
^2+16*vr*x^3*vx*y^2*z^2+16*vr*y^3*vy*x^2*z^2+2*x*vy^2*z^8+4*y^4*x*z^2-
4*x*z^4*y^2+4*vr*x*vx*y^6+8*vr*x*vx*y^4*z^2+4*vr*x*vx*y^2*z^4+4*vr*y*v
y*x^6+4*x*vx^2*y^2*z^4-8*x^3*vx^2*y^4*z^2-12*x*vx^2*y^4*z^4-4*x^3*vx^2
*y^2*z^4-4*x*vx^2*y^2*z^6+20*y*vx*vy*x^6+8*y^5*vx*vy*z^2+16*y^3*vx*vy*
z^4+8*y*vx*vy*z^6+8*vr*y*vy*x^4*z^2+4*vr*y*vy*x^2*z^4+32*y*vx*vy*x^4*z
^2+20*y*vx*vy*x^2*z^4-8*vr*r*y^3*vy*x^2+13*x*vx^2*r*y^4*z^2-6*y*vx*vy*
r*z^6-16*y*vx*vy*r*x^4+9*x*vy^2*r*y^4*z^2+3*x*vy^2*r*y^2*z^4-6*vr*r*x^
5*vx*z^2+2*vr*r*x^3*vx*y^4-2*vr*r*x^7*vx-2*vr*r*y^7*vy-3*x^5*vx^2*r*z^
2+5*x^3*vx^2*r*y^4-3*x^3*vx^2*r*z^4-6*x^5*vx^2*r*y^2+10*x*vx^2*r*y^6-x
*vx^2*r*z^6+8*x^3*vx^2*r*y^2+8*x^3*vy^2*r*z^2-17*x^5*vy^2*r*z^2+4*x^3*
vy^2*r*y^4-10*x^3*vy^2*r*z^4-9*x^5*vy^2*r*y^2+5*x*vy^2*r*y^6-x*vy^2*r*
z^6+8*x^3*vy^2*r*y^2+2*r*x^3*z^2*y^2-y^4*r*x*z^2+x*r*z^4*y^2+3*r*x^5*z
^2+8*y^2*r*x^3-4*y^4*r*x-8*x^7*vy^2*r+16*r*x^3*z^2+4*x*r*z^4-x^7*vx^2*
r+3*r*x^3*z^4+x^5*r*y^2-y^4*r*x^3-y^6*r*x+x*r*z^6+12*x^5*r+x^7*r-4*x^3
*vx^2*r*y^2*z^2-8*vr*r*x^3*vx*y^2-2*vr*r*x*vx*y^6-6*vr*r*x*vx*y^4*z^2-
4*vr*r*x^3*vx*y^2*z^2-6*vr*r*x*vx*y^2*z^4-6*vr*r*y*vy*x^4*z^2-4*vr*r*y
^3*vy*x^2*z^2-6*vr*r*y*vy*x^2*z^4-2*vr*r*x*vx*z^6-2*vr*r*y*vy*x^6+2*vr
*r*y^5*vy*x^2-6*vr*r*y^5*vy*z^2-6*vr*r*y^3*vy*z^4-2*vr*r*y*vy*z^6+2*x*
vx^2*r*y^2*z^4-16*y*vx*vy*r*x^2*z^2-6*y*vx*vy*r*x^4*z^2+6*y^5*vx*vy*r*
x^2-6*y*vx*vy*r*x^2*z^4-6*y^5*vx*vy*r*z^2-12*y^3*vx*vy*r*z^4-6*x^3*vy^
2*r*y^2*z^2)/(5*y^8*z^2+y^8*x^2-2*y^4*x^6+44*x^4*y^2*z^2+44*x^2*y^4*z^
2+28*x^2*y^2*z^4+4*x^8+4*y^8+20*x^6*y^2+12*x^6*z^2+32*x^4*y^4+12*x^4*z
^4+20*x^2*y^6+4*x^2*z^6+12*y^6*z^2+12*y^4*z^4+4*y^2*z^6+y^10+x^10+z^10
+6*y^4*x^4*z^2+8*y^6*x^2*z^2+18*y^4*x^2*z^4+8*x^6*z^2*y^2+18*x^4*z^4*y
^2+16*x^2*z^6*y^2+x^8*y^2+10*x^6*z^4-2*y^6*x^4+10*x^4*z^6+5*x^8*z^2+5*
x^2*z^8+10*y^4*z^6+5*y^2*z^8+10*y^6*z^4-8*y^2*r*x^2*z^2-28*r*x^2*y^2*z
^4-26*r*x^2*y^4*z^2-26*r*x^4*y^2*z^2-8*y^2*r*x^4-8*y^4*r*x^2-14*r*x^6*
z^2-8*r*x^4*y^4-18*r*x^4*z^4-10*r*x^2*z^6-14*r*y^6*z^2-18*r*y^4*z^4-10\

I haven't checked this for accuracy, and I don't plan to - actually, it looks to me like the whole thing didn't paste correctly, it was probably too long for the paste buffer....

Homework Helper
pervect said:
Finally, since I have the automaton working, here's the sort of result the above procedure gets when applied to the psueudo-cartesian metric you suggested earlier.
It appears I've got MAXIMA running. To crunch out the equations of motion by varying a path, the square root is somewhat ugly. Consequently, it helps to substitute $$L^2$$ for L in the Euler-Lagrange equations. One begins by defining:

$$I = \left(\frac{ds}{dt}\right)^2$$

then the thing which is to be minimized is $$\int \sqrt{I} \;dt.$$ The usual Euler-Lagrange equations minimize $$\int L\;dt.$$ To convert them to the square root form, simply substitute and obtain:

$$\frac{d}{dt}\left(\frac{\partial I}{\partial \dot{x}}\right) = \frac{1}{2I}\frac{\partial I}{\partial \dot{x}} + \frac{\partial I}{\partial x},$$

and similarly for y. The equations work out a lot simpler if one abbreviates sqrt(x^2+y^2) as r. While I've computed the equations of motion once, I'm not going to post the results until tomorrow as the tools are new to me and I'd like to go through the exercise again. But I have no doubt that the Maxima software will save me much grief:

http://maxima.sourceforge.net/

Here it is the next day and I'm too busy to deal with this. One thing I'd like to mention, working out the equations of motion in coordinate time this way has the advantage that the resulting geodesics can be massless as well as massive.

Two of Einstein's principles were simplicity and the use of geometry. It seems that there are three different things that one could require to be simple. One could require the equations of motion, the energies / lagrangian, or one could require an overall principle be simple. What Einstein did was to require the overall principle be simple.

I'm not at all convinced that this is the correct choice. Maybe the Gravity Probe B will give us more information on this. Newton's gravitation is simple in all three parts, the equations of motion, the Lagrangian, and the overall principle (linearity). I know that the equations of motion have to be a part of reality, but the others might just be parts of our mathematics. So my inclination is to require that the equations of motion be simple. Accordingly, I'll be looking for equations of motion that are simple rather than Lagrangians or principles that are simple.

Carl

Last edited:
pervect
Staff Emeritus
CarlB said:
It appears I've got MAXIMA running. To crunch out the equations of motion by varying a path, the square root is somewhat ugly. Consequently, it helps to substitute $$L^2$$ for L in the Euler-Lagrange equations. One begins by defining:

$$I = \left(\frac{ds}{dt}\right)^2$$

then the thing which is to be minimized is $$\int \sqrt{I} \;dt.$$ The usual Euler-Lagrange equations minimize $$\int L\;dt.$$ To convert them to the square root form, simply substitute and obtain:

$$\frac{d}{dt}\left(\frac{\partial \sqrt{I}}{\partial \dot{x}}\right) = \frac{d}{d t} \left[ \frac{1}{2 \sqrt{I}} \left( \frac{\partial I}{\partial \dot{x}} \right) \right] = \frac{1}{2 \sqrt{I}} \frac{\partial I}{\partial x}$$
Good luck with this approach. (I think there was a small typo which I changed, see above).

Minimizing the action sqrt(I) [in some sign conventions, sqrt(-I)] gives you a parameter independent path. However, I think you'll find when you actually try this direct route is that it gets very messy.

Minimizing the action I (rather than the square root) is the usual route to the geodesic equation. While the path generated by the previous approach was independent of the parametrization of the path, the path generated by the approach using the action I (rather than its square root) requires that the path be parameterized by proper time or an affine transformation of proper time.

You might want to consult a textbook on GR for more verification if you are interested in the details.

I expect that you'll find that minimzing I, and then eliminating proper time as a separate step, will be computationally much more convenient than your direct approach using sqrt(I).

Last edited:
pervect
Staff Emeritus
BTW, Maxima's ctensor() package will calculate the geodesic equations for you:

For example (/* are comments */):

ct_coordsys(exteriorschwarzschild);
//* calculate and setup stuff such as inverse metric */
cmetric();
//* display some stuff */
//* coordinates */
ct_coords;
//* metric */
lg;
//* inverse metric */
ug;
//* the geodesic equations */
cgeodesic(true);

Homework Helper
pervect said:
Good luck with this approach. (I think there was a small typo which I changed, see above).
Thanks for pointing this out. I've run out of all possible other things worth doing at the moment, so time to run this calculation out. Let me begin where you left off:

$$\frac{d}{d t} \left[ \frac{1}{2 \sqrt{I}} \left( \frac{\partial I}{\partial \dot{x}} \right) \right] = \frac{1}{2 \sqrt{I}} \frac{\partial I}{\partial x}$$

Multiply by 2 to get rid of the halves. I want to eliminate the sqrt(I). Accordingly, use the product rule on the time derivative:

$$-\frac{1}{2 I^{3/2}}\left(\frac{\partial I}{\partial \dot{x}}\right)\frac{dI}{dt} + \frac{1}{I^{1/2}}\frac{d}{dt}\left(\frac{\partial I}{\partial \dot{x}}\right) = \frac{1}{I^{1/2}}\frac{\partial I}{\partial x}$$

Bring all the terms to the LHS and multiply by I sqrt(I) to get:

$$I\frac{d}{dt}\left(\frac{\partial I}{\partial \dot{x}}\right) -\frac{1}{2}\left(\frac{\partial I}{\partial \dot{x}}\right)\frac{dI}{dt} - I\frac{\partial I}{\partial x} = 0$$

Write out the derivatives with respect to time in terms of partial derivatives:

$$I\left(\frac{\partial^2 I}{\partial \dot{x}\partial x}\right)\dot{x} + I\left(\frac{\partial^2 I}{\partial \dot{x}^2}\right) \ddot{x} + I\left(\frac{\partial^2 I}{\partial \dot{x}\partial y}\right) \dot{y} + I\left(\frac{\partial^2 I}{\partial \dot{x}\partial \dot{y}}\right) \ddot{y}$$
$$-\frac{1}{2}\left(\frac{\partial I}{\partial \dot{x}}\right)\left(\frac{\partial I}{\partial x}\right)\dot{x} -\frac{1}{2}\left(\frac{\partial I}{\partial \dot{x}}\right)\left(\frac{\partial I}{\partial \dot{x}}\right)\ddot{x} -\frac{1}{2}\left(\frac{\partial I}{\partial \dot{x}}\right)\left(\frac{\partial I}{\partial y}\right)\dot{y} -\frac{1}{2}\left(\frac{\partial I}{\partial \dot{x}}\right)\left(\frac{\partial I}{\partial \dot{y}}\right)\ddot{y} -I\left(\frac{\partial I}{\partial x}\right) = 0$$

Collect terms according to type in the phase coordinates:

$$\begin{array}{rcl} &&\left[I\left(\frac{\partial^2 I}{\partial \dot{x}^2}\right) -\frac{1}{2}\left(\frac{\partial I}{\partial \dot{x}}\right)\left(\frac{\partial I}{\partial \dot{x}}\right)\right] \ddot{x} + \left[I\left(\frac{\partial^2 I}{\partial \dot{x}\partial \dot{y}}\right) -\frac{1}{2}\left(\frac{\partial I}{\partial \dot{x}}\right)\left(\frac{\partial I}{\partial \dot{y}}\right)\right] \ddot{y}\\ &=&\left[\frac{1}{2}\left(\frac{\partial I}{\partial \dot{x}}\right)\left(\frac{\partial I}{\partial x}\right) -I\left(\frac{\partial^2 I}{\partial \dot{x}\partial x}\right) \right]\dot{x} + \left[\frac{1}{2}\left(\frac{\partial I}{\partial \dot{x}}\right)\left(\frac{\partial I}{\partial y}\right) -I\left(\frac{\partial^2 I}{\partial \dot{x}\partial y}\right) \right] \dot{y} + I\left(\frac{\partial I}{\partial x}\right) \end{array}$$

A similar equation obtains for y:

$$\begin{array}{rcl} &&\left[I\left(\frac{\partial^2 I}{\partial \dot{y}\partial \dot{x}}\right) -\frac{1}{2}\left(\frac{\partial I}{\partial \dot{y}}\right)\left(\frac{\partial I}{\partial \dot{x}}\right)\right] \ddot{x} + \left[I\left(\frac{\partial^2 I}{\partial \dot{y}^2}\right) -\frac{1}{2}\left(\frac{\partial I}{\partial \dot{y}}\right)\left(\frac{\partial I}{\partial \dot{y}}\right)\right] \ddot{y}\\ &=&\left[\frac{1}{2}\left(\frac{\partial I}{\partial \dot{y}}\right)\left(\frac{\partial I}{\partial y}\right) -I\left(\frac{\partial^2 I}{\partial \dot{y}\partial y}\right) \right]\dot{y} + \left[\frac{1}{2}\left(\frac{\partial I}{\partial \dot{y}}\right)\left(\frac{\partial I}{\partial x}\right) -I\left(\frac{\partial^2 I}{\partial \dot{y}\partial x}\right) \right] \dot{x} + I\left(\frac{\partial I}{\partial y}\right)$$

Carl

Last edited: