Schwarzschild Orbits in Cartesian coordinates

[CarlB]

My Java applet gravity simulator http://www.gaugegravity.com/testapplet/SweetGravity.html
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

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

 

[CarlB]

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$$

 

[Jorrie]

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

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

                        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"

 

[pervect]

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 therefore 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.

 

[Jorrie]

Cartesian coordinates

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.

 

[CarlB]

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]

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]

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.

 

[pervect]

The Christoffel symbols for the Painleve metric are:

                                   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.

 

[CarlB]

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.

[edit] 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

 

[CarlB]

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

 

[pervect]

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.

 

[CarlB]

$$\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

 

[CarlB]

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

Here's where I begin:

$$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$$

 

[CarlB]

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.

 

[pervect]

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

 

[pervect]

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.


>qload(schw);
>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]

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]

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.

 

[pervect]

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...

 

[CarlB]

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

 

[pervect]

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).

 

[pervect]

BTW, Maxima's ctensor() package will calculate the geodesic equations for you:

For example (/* are comments */):

load(ctensor);
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);

 

[CarlB]

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

 

[CarlB]

The above equations can be written as:

$$\begin{array}{rcl} A &=& 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) \\ B &=& 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) \\ C &=& \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) \\ D &=& 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) \\ E &=& 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) \\ F &=& \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)$$

so that

$$A \ddot{x} + B \ddot{y} = C$$
$$D \ddot{x} + E \ddot{y} = F$$

The solution is:

$$\begin{array}{rcl} \ddot{x} &=& (B*F - C*E)/(B*D-A*E),\\ \ddot{y} &=& (C*D - A*F)/(B*D-A*E). \end{array}$$

To compute A, B, C, D, E, and F, I need the following terms:
$$\begin{array}{cccc} \left(\frac{\partial I}{\partial x}\right) & \left(\frac{\partial I}{\partial \dot{x}}\right) & \left(\frac{\partial I}{\partial y}\right) & \left(\frac{\partial I}{\partial \dot{y}}\right) \\ \left(\frac{\partial^2 I}{\partial \dot{x}^2}\right) & \left(\frac{\partial^2 I}{\partial \dot{y}^2}\right) & \left(\frac{\partial^2 I}{\partial \dot{x} \partial \dot{y}}\right) & \\ \left(\frac{\partial^2 I}{\partial \dot{x} \partial x}\right) & \left(\frac{\partial^2 I}{\partial \dot{y} \partial y}\right) & \left(\frac{\partial^2 I}{\partial \dot{x} \partial y}\right) & \left(\frac{\partial^2 I}{\partial \dot{y} \partial x}\right) \end{array}$$

It remains to compute these from the value of I for the Schwarzschild and Painleve coordinates.

Carl

 

[CarlB]

As a test of this solution to the equations of motion, let me see what happens when I use it to derive Newtonian gravitation. I need to minimize the following function over the path:

$$L(x,y,\dot{x},\dot{y}) = T - V = \frac{\dot{x}^2 + \dot{y}^2}{2} + \frac{1}{r}$$

Writing this in $$\sqrt{I}$$ form, I have that

$$I = I_N = \left(\frac{\dot{x}^2 + \dot{y}^2}{2} + \frac{1}{r}\right)^2$$

The first partial derivatives are:

$$\frac{\partial I}{\partial x} = -\left(\dot{x}^2 + \dot{y}^2 + \frac{2}{r}\right)\frac{x}{r^3}$$
$$\frac{\partial I}{\partial y} = -\left(\dot{x}^2 + \dot{y}^2 + \frac{2}{r}\right)\frac{y}{r^3}$$
$$\frac{\partial I}{\partial \dot{x}} = 2\left(\dot{x}^2 + \dot{y}^2 + \frac{2}{r}\right)\dot{x}$$
$$\frac{\partial I}{\partial \dot{y}} = 2\left(\dot{x}^2 + \dot{y}^2 + \frac{2}{r}\right)\dot{y}$$

The second order partial derivative with respect to momenta are:
$$\frac{\partial^2 I}{\partial \dot{x}^2} = \left(6\dot{x}^2 + 2\dot{y}^2 + \frac{4}{r}\right)$$
$$\frac{\partial^2 I}{\partial \dot{y}^2} = \left(6\dot{y}^2 + 2\dot{x}^2 + \frac{4}{r}\right)$$
$$\frac{\partial^2 I}{\partial \dot{x} \partial \dot{y}} = 4 \dot{x}\dot{y}$$

And finally, the four mixed second order partial derivatives are:
$$\frac{\partial^2 I}{\partial \dot{x} \partial x} = -\frac{4 \dot{x} x}{r^3}$$
$$\frac{\partial^2 I}{\partial \dot{y} \partial y} = -\frac{4 \dot{y} y}{r^3}$$
$$\frac{\partial^2 I}{\partial \dot{x} \partial y} = -\frac{4 \dot{x} y}{r^3}$$
$$\frac{\partial^2 I}{\partial \dot{y} \partial x} = -\frac{4 \dot{y} x}{r^3}$$

It's clear that these are going to simplify wonderfully when I substitute them in for A, B, C, D, E, and F.

Carl

 

[pervect]

Are you still working on this? Given that I presented an answer, quite a long time ago, by using the standard geodesic equations and then eliminating proper time, I would have thought that you would have been able to confirm or deny the correctness of my answer by now.

Specifically, I'm talking about the answers in:
https://www.physicsforums.com/showpost.php?p=1046874&postcount=17
for the Schwarzschild metric and
https://www.physicsforums.com/showpost.php?p=1047041&postcount=20
for the Painleve metric.

The prorcess of eliminating proper time turns out to be quite simple, BTW, if approached in the right manner. It turns out to be a matter of just subtracting a multiple of one quadtratic form from another quadratic form, as you'll see if you actually attempt this method.

I'm reasonably sure the answers I gave in these two posts are correct, BTW. but I haven't done a lot of testing, just one simulation that seemed to conserve the quantities that should be conserved. (And there's always the possibility of typos).

Interestingly enough, there is an alternate method available that was mentioned in another thread, called the Einbein. See Self-adjoint's thread of that title for more info. This accomplishes the same result as using "proper time" by using a different auxillary variable which accomplishes much the same purpose.

I'll also standy by remarks that changing to "cartesian-like" x and y coordinates before solving the equations of motion is Not A Good Idea. See

https://www.physicsforums.com/showpost.php?p=1047056&postcount=21

for why.

So in concluion, I think I've done about all I can possibly do to help you with this problem, to the point of actually giving you a correct solution for the Schwarzschild and Painleve metrics.

The one thing more that I could do is to recommend some specific textbooks, which I'd be quite willing to do. However, as you are not a GR person and probably don't have GR textbooks in your library, I'm not sure that that would be all that helpful. If you have some questions about how the standard geodesic equations are derived, this issue should be covered in most GR textbooks.

I imagine that particle physicists might have their own approach to the problem (like the Einbein I mentioned) as well. In any event, I don't think there's much more I can do for you (except perhaps to give some textbook references, which won't help that much if you don't have the textbooks in question) having already solved the problem and having presented the solution.

 

[jbusc]

Not to take away from this very interesting discussion, but are the trig functions really that slow? They used to be in older versions of Java but are fairly optimized now in modern JVM's...

 

[CarlB]

Are you still working on this? Given that I presented an answer, quite a long time ago, by using the standard geodesic equations and then eliminating proper time, I would have thought that you would have been able to confirm or deny the correctness of my answer by now.

Business has been very busy these past few weeks and I haven't spent much time on it, but yes, I'm still working on this. Eventually I will get around to checking your equation.

Interestingly enough, there is an alternate method available that was mentioned in another thread, called the Einbein.

I've been following this with interest.

I'll also standy by remarks that changing to "cartesian-like" x and y coordinates before solving the equations of motion is Not A Good Idea.

You're wrong.

The one thing more that I could do is to recommend some specific textbooks, which I'd be quite willing to do. However, as you are not a GR person and probably don't have GR textbooks in your library, I'm not sure that that would be all that helpful.

Of course I have MTW which covers the problem in the mode you've given at some length. I took a one quarter (graduate) GR class at the Univ. of California Irvine. I had the highest grade in the class and earned the "A+" that were assigned to the top student in a class at that time. Probably the reason I did very well in physics grad school was because I had a masters degree in mathematics before I switched to physics.

Carl

 

[CarlB]

Not to take away from this very interesting discussion, but are the trig functions really that slow? They used to be in older versions of Java but are fairly optimized now in modern JVM's...

It's partly speed and partly elementary particles stuff that I am not ready to explain. Java allows you to write subroutines in other languages (i.e. assembler or C), so you can write blindingly fast integration subroutines.

Example, AMD-64 CPU with 64-bit arithmetic. The latency of a 64-bit floating multiply or add (FMUL or FADD) is 4 clocks, while the latency of the 64-bit floating cosine/sine (FSINCOS) is 104. For these numbers, see
http://www.amd.com/us-en/assets/content_type/white_papers_and_tech_docs/25112.PDF

Carl

 

[pervect]

You're wrong.

Carl

The geodesic equation for your metric isn't any ordinary rabbit. It's got teeth like - it can leap like - look at the bones man! Look at the pages of output from Maxima's built in cgeodesic() function!

I warned ye! I try and try and try and tell them, but does anybody ever listen to Tim? Nooooo. "It's just a little bunny rabit", they say. "Go and chop it's head off", they say.

Oh well, I'll go put my "no" vote in on the poll now...

 

[CarlB]

Looks like an extra factor of 2 appeared in a previous post, let me give the corrected values for the partial derivatives with the Newtonian gravitation (as computed in sqrt(I) form):

$$I = I_N = \left(\frac{\dot{x}^2 + \dot{y}^2}{2} + \frac{1}{r}\right)^2 = (\sqrt{I})^2$$

$$\frac{\partial I}{\partial x} = -\sqrt{I}\left(\frac{2x}{r^3}\right)$$
$$\frac{\partial I}{\partial y} = -\sqrt{I}\left(\frac{2y}{r^3}\right)$$
$$\frac{\partial I}{\partial \dot{x}} = \sqrt{I}(2\dot{x})$$
$$\frac{\partial I}{\partial \dot{y}} = \sqrt{I}(2\dot{y})$$

The second order partial derivative with respect to momenta are:
$$\frac{\partial^2 I}{\partial \dot{x}^2} = \left(3\dot{x}^2 + \dot{y}^2 + \frac{2}{r}\right)$$
$$\frac{\partial^2 I}{\partial \dot{y}^2} = \left(3\dot{y}^2 + \dot{x}^2 + \frac{2}{r}\right)$$
$$\frac{\partial^2 I}{\partial \dot{x} \partial \dot{y}} = 2 \dot{x}\dot{y}$$

And finally, the four mixed second order partial derivatives are:
$$\frac{\partial^2 I}{\partial \dot{x} \partial x} = -\frac{2 \dot{x} x}{r^3}$$
$$\frac{\partial^2 I}{\partial \dot{y} \partial y} = -\frac{2 \dot{y} y}{r^3}$$
$$\frac{\partial^2 I}{\partial \dot{x} \partial y} = -\frac{2 \dot{x} y}{r^3}$$
$$\frac{\partial^2 I}{\partial \dot{y} \partial x} = -\frac{2 \dot{y} x}{r^3}$$

 

[CarlB]

Now, substituting the above into post #26 gives:

$$\begin{array}{rcl}A &=& 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) \\ &=& I(3\dot{x}^2+\dot{y}^2+2/r)-0.5(2\dot{x}\sqrt{I})(2\dot{x}\sqrt{I}) \\ &=& I(\dot{x}^2+\dot{y}^2 +2/r)\end{array}$$

$$\begin{array}{rcl}B &=& 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) \\ &=& I (2\dot{x}\dot{y})-0.5(2\sqrt{I}\dot{x})(2\sqrt{I}\dot{y}) \\ &=& 0\end{array}$$

$$\begin{array}{rcl}C &=& \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) \\ &=& (0.5(2\sqrt{I}\dot{x})(-2\sqrt{I}x/r^3) + I2\dot{x}x/r^3)\dot{x} + (0.5(2\sqrt{I}\dot{x})(-2\sqrt{I}y/r^3) +I2\dot{x}y/r^3)\dot{y} + I\sqrt{I}(-2x/r^3) \\ &=& -I(\dot{x}^2+\dot{y}^2 +2/r )x/r^3 \end{array}$$

$$\begin{array}{rcl} D &=& 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) \\ &=& I(2\dot{x}\dot{y})-0.5(2\sqrt{I}\dot{y})(2\sqrt{I}\dot{x}) \\ &=& 0\end{array}$$

$$\begin{array}{rcl} E &=& 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) \\ &=& I(6\dot{y}^2+2\dot{x}^2+4/r) -0.5(4\sqrt{I}\dot{y})^2 \\&=& I(-2\dot{y}^2+2\dot{x}^2+4/r)\end{array}$$

$$\begin{array}{rcl} F &=& \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) \\ &=& (0.5(2\sqrt{I}\dot{y})(-2\sqrt{I}y/r^3) + I2\dot{y}y/r^3)\dot{y} + (0.5(2\sqrt{I}\dot{y})(-2\sqrt{I}x/r^3) + I2\dot{y}x/r^3)\dot{x} +I\sqrt{I}(-2y/r^3) \\ &=& -I(\dot{x}^2+\dot{y}^2+2/r)y/r^3\end{array}$$

Since A, B, C, D, E and F are each written in the above as multiples of I, we can factor this out to get A', B', C', D', E' and F':

$$\begin{array}{rcl} A' = E' &=& (\dot{x}^2+\dot{y}^2 +2/r) \\ B' = D &=& 0 \\ C' &=& -(\dot{x}^2+\dot{y}^2 +2/r )x/r^3 \\ F' &=& -(\dot{x}^2+\dot{y}^2 +2/r )y/r^3\end{array}$$


And the equations of motion are now simply:

$$\begin{array}{rcccl}\ddot{x} &=& C'/A' &=& -x/r^3,\\ \ddot{y} &=& F'/E' &=& -y/r^3,\end{array}$$

which are the correct equations of motion for Newton's gravitation. Thus the method works. Next, the Schwarzschild metric.

Carl

 

[CarlB]

Okay, it's time for me to push the Schwarzschild metric of post #15 through the machinery developed in post #26. I was going to do this the other night but PF crashed and I deleted the post.

Rewriting the "I" of post #15 (and removing a factor of r that unnaturally appeared in post #2) gives a particularly simple form:

$$I = 1-\frac{2}{r}-(\dot{x}^2+\dot{y}^2) -2\frac{(x\dot{x}+y\dot{y})^2}{r^3-r^2}$$

The first order partial derivatives are:

$$\left(\frac{\partial I}{\partial x}\right) = \frac{2x}{r^3}-\frac{4\dot{x}(x\dot{x}+y\dot{y})}{r^3-r^2} +2\frac{(x\dot{x}+y\dot{y})^2 (3rx-2x)}{(r^3-r^2)^2}$$

$$\left(\frac{\partial I}{\partial y}\right) = \frac{2y}{r^3}-\frac{4\dot{y}(x\dot{x}+y\dot{y})}{r^3-r^2} +2\frac{(x\dot{x}+y\dot{y})^2 (3ry-2y)}{(r^3-r^2)^2}$$

$$\left(\frac{\partial I}{\partial \dot{x}}\right) = -2\dot{x}-\frac{4x(x\dot{x}+y\dot{y})}{r^3-r^2}$$

$$\left(\frac{\partial I}{\partial \dot{y}}\right) = -2\dot{y}-\frac{4y(x\dot{x}+y\dot{y})}{r^3-r^2}$$

It doesn't look like the 2nd order partials will be particularly difficult. However, this busy old man is going to call it a night early.

Carl