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 [tex]d^2x/dt^2[/tex] and [tex]d^2y/dt^2[/tex] 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
Step (1), write the Schwarzschild metric in Cartesian coordinates. [tex]x = r \cos(\phi)[/tex] [tex]y = r \sin(\phi)[/tex] [tex]dx = \cos(\phi) dr - r\sin(\phi) d\phi[/tex] [tex]dy = \sin(\phi) dr + r\cos(\phi) d\phi[/tex] [tex]dr = \cos(\phi) dx + \sin(\phi) dy[/tex] [tex]rd\phi = -\sin(\phi) dx + \cos(\phi) dy[/tex] [tex]dr^2 = \cos^2(\phi) dx^2 + \sin^2(\phi) dy^2 + \sin(2\phi) dx\;dy[/tex] [tex]r^2d\phi^2 = \sin^2(\phi) dx^2 + \cos^2(\phi) dy^2 - \sin(2\phi) dx\;dy[/tex] [tex]dr^2 = (x^2dx^2 + y^2dy^2 + 2xy dx\;dy)/r^2[/tex] [tex]r^2d\phi^2 = (y^2dx^2 + x^2dy^2 - 2xy dx\;dy)/r^2[/tex] [tex]ds^2 = -\frac{(r-2)}{r}dt^2 + \frac{r}{r-2}dr^2 + r^2d\phi^2[/tex] [tex]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)[/tex] [tex]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)[/tex] [tex]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)[/tex] [tex]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[/tex]
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."?
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 (Text): 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"
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 [tex] \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 [/tex] [tex] \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 [/tex] [tex] \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 [/tex] See the previous post for the values for the Christoffel symbols. Note for instance that since r^2 [itex]{d \phi}/ {d tau}[/itex] = constant, when differentiated, gives by the chain rule [tex] 2 r \frac{dr}{d\tau} \frac{d \phi}{d tau} + r^2 \frac{d^2 \phi} {d tau^2} = 0 [/tex] 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.
Cartesian 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.
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
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 [tex] \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} [/tex] 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
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 [tex] \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} [/tex] 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.
The Christoffel symbols for the Painleve metric are: Code (Text): 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.
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 [tex]r^2 = x^2 + y^2 + z^2\;\;\;\;\;\;\;(6)[/tex] http://arxiv.org/abs/gr-qc/9505010 [/edit] Carl
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
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.
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 [tex]K = \dot{\phi} r^2 / (1-2m/r)[/tex] 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
Step (2), write the proper length of a path as an integral over coordinate time. Here's where I begin: 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: [tex]\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}[/tex] 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: [tex] S = \int_{t_1}^{t_2} \frac{ds}{dt}\; dt[/tex] [tex] = \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[/tex]
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, [tex]\frac{d}{dt} \left(\frac{\partial ds/dt}{\partial x}\right) = \frac{\partial ds/dt}{\partial \dot{x}}[/tex] Carl Okay, now I've got Maxima running on my laptop. Pretty cool. And cheap.
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: [tex] \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) [/tex] [tex] \frac{d^2 phi}{d t^2} = -\frac {2 {\it vr}\,{\it vphi}\, \left( r -3\,m \right) }{ \left( r-2\,m \right) r} [/tex] vr = dr/dt and vphi = dphi/dt
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);
The standard textbook approach is to write the Lagrangian as [tex] L = g_{\mu \nu} \frac{d x^\mu}{d \lambda} \frac{d x^\nu}{d \lambda} d \lambda [/tex] (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.
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 [tex] \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) [/tex] [tex] \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) [/tex] 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.