•  Trajectory Near a Black Hole

  Trajectory Near a Black Hole:
an application of Maple  to Lagrangian mechanics  

Frank Wang and J. F. Ogilvie  @ 2002 March

     The lagrangian formulation of mechanics has great advantages in practical use.  Mathematics involved in this subject is the calculus of variations, which is to find a function that minimizes an integral.  Calculations of this type require finding the derivative of a function with respect to another function.  There are many procedures in Maple  to support such differentiation, and there are many integrated procedures to derive the equations of motion directly.  These procedures are tremendously convenient in practical use, but the shortcomings are that one needs to invoke external procedures for Maple  to perform such calculation, and that perhaps a novice user tends not to comprehend fully the underlying physics.  We develop a method in Maple  to solve problems of the calculus of variations.  In our approach, we use basic commands subs , diff , and dsolve  only, without programming or invoking an external library.  Although our method is a pedagogic approach that might involve longer steps, it is a straightforward attack on this problem, and practically all problems in classical mechanics can be solved once the lagrangian is found.  In most real physics problems, there are no analytic solutions to differential equations.  We particularly emphasize forming plots based on a numerical method; by varying physical parameters and initial conditions and observing the graphic outputs, one can develop a sense of intuition of underlying physics. We introduce an example in general relativity, to find the trajectory of a particle near a black hole, which corresponds to nothing other than the shortest connection between two points in a curved space.

subject    calculus of variations, lagrangian, Schwarzschild metric, geodesic equation

 Trajectory Near a Black Hole

    John Wheleer summarized general relativity with the statement "geometry tells matter how to move, matter tells geometry how to curve".  This statement signifies that the presence of mass will curve the space, and the trajectory of a massive particle in a curved space is the shortest connection between two points.  We have solved the problem of finding the shortest connection between two points in a plane, and for a curved space it is exactly the same mathematics of applying the Euler-Lagrange equation; the differential equations therefrom are called geodesic equations.  

One solution of general relativity coresponding to the exterior of a spherically symmetric gravitational source  is the Schwarzschild metric,

which indicates a curved space.

     Similar to the example of finding the shortest plane curve connecting two points, we need to minimize the function

 

It can be proved that, if at one instant  is , it retains this value.  We can then rewrite the metric as

 

We let  be the independent variable; ,  and  become the dependent variables.  Therefore,

Identifying the function  as

.

because it is equivalent to find the minimum of  as .  

>

restart:

Defining

>	f := 1/2*(-(1 - 2*M/r(s))*diff(t(s), s)^2     + 1/(1 - 2*M/r(s))*diff(r(s), s)^2 + r(s)^2*diff(phi(s),s)^2);

and making the substitutions,

>	f1 := subs({t(s)=var1, diff(t(s),s)=var2, r(s)=var3, diff(r(s),s)=var4,             phi(s)=var5, diff(phi(s),s)=var6}, f);

we see that there is no explicit  term; there is a constant, which we denote as .

>	Eq11 := diff(f1, var2);

>	Eq12 := diff(f1, var1);

>	Eq13 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var4=diff(r(s),s), var5=phi(s), var6=diff(phi(s),s)}, Eq11);

>	Eq14 := Eq13 = -a;

Applying the Euler-Lagrange equation to the  coordinate yields

>	Eq21 := diff(f1, var4);

>	Eq22 := diff(f1, var3);

>	Eq23 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var4=diff(r(s),s), var5=phi(s), var6=diff(phi(s),s)}, Eq21);

>	Eq24 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s),         var4=diff(r(s),s), var5=phi(s), var6=diff(phi(s),s)}, Eq22);

>	Eq25 := diff(Eq23, s);

>	Eq26 := Eq25 - Eq24 = 0;

Also there is no explicit  term, so we denote another constant

>	Eq31 := diff(f1, var6);

>	Eq32 := diff(f1, var5);

>	Eq33 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s),         var4=diff(r(s),s), var5=phi(s), var6=diff(phi(s),s)}, Eq31);

>	Eq34 := Eq33 = h;

We list three equations of motion.

>

Eq14;

>

Eq26;

>

Eq34;

We make some rearrangements to decouple the equations.

>	Eq53 := isolate(Eq14, diff(t(s),s));

>	Eq54 := isolate(Eq34, diff(phi(s),s));

>	Eq55 := subs({Eq53, Eq54}, Eq26);

We list the   and   equations after rearrangement.

>

Eq54;

>

Eq55;

We solve these differential equations numerically.

>

M := 1;

>	Eq77 := r(0) = 11;

>	Eq78 := D(r)(0) = 0;

>	Eq79 := phi(0) = 0;

>	Eq80 := D(phi)(0) = 0.0295;

The constants  and  that we assigned earlier are related to energy and angular momentum; they can be evaluated from initial conditions.  

>	h := rhs(Eq77)^2*rhs(Eq80);

>	a := sqrt((1 - 2*M/rhs(Eq77))*(1 + h^2/rhs(Eq77)^2));

We make a plot of "effective potential", which is unrelated to our approach.  We list it so that the reader can compare with a convential approach such as in Misner, Thorne and Wheeler's Gravitation  (1973), p. 660.

>	Vsq := sqrt((1 - 2*M/r)*(1 + h^2/r^2)):

>	plot({Vsq, a}, r=2..25, a-0.01*a..a+0.01*a);

>	ini := Eq77, Eq78, Eq79;

>	Eq91 := dsolve({Eq54, Eq55, ini}, {r(s), phi(s)}, numeric,                output=listprocedure);

>	with(plots):

Warning, the name changecoords has been redefined

We can plot the trajectory of a massive particle in a strong gravitional field.

>	polarplot([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), s=0..500],           scaling=constrained);

This trajectory corresponds to nothing other than the shortest connection between two points in a curved space.  To visualize a curved space, it is convenient to use an embedding diagram.  This diagram is the slice of the equatorial plane, with a depth  to accommodate the actual length.  The embedding formula for this problem is

For details see ibid p. 615.  

>	p1 := plot3d([x*cos(t), x*sin(t), sqrt(8*(x-2))], x=3..14,              t=0..2*Pi, style=hidden):

>	p2 := spacecurve([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]),                   sqrt(8*(rhs(Eq91(s)[3])-2)), s=0..500], coords=cylindrical, color=black, numpoints=400):

>	display([p1, p2]);

We animate the motion by displaying plot frames in sequence.

>	N := 50; nstep := 10;

>

for i from 0 by nstep to N*nstep do   pl[i] := polarplot([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), s=i..i+nstep]):   pl3d[i] := spacecurve([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]),                   sqrt(8*(rhs(Eq91(s)[3])-2)), s=i..i+nstep], coords=cylindrical, color=black):   plsum[i] := display(seq(pl[j*nstep], j=0..i/nstep)):   pl3dsum[i] := display({p1, seq(pl3d[j*nstep], j=0..i/nstep)}): end do:

>	display([seq(plsum[i*nstep], i=0..N)], insequence=true,          scaling=constrained);

>	display([seq(pl3dsum[i*nstep], i=0..N)], insequence=true);

Maple  has a tensor  package that can calculate geometric properties in curved space and derive the geodesic equations in an integrated command.  

>

>

References

1. H. Goldstein, Classical Mechanics , 2nd ed., Reading, Mass: Addison-Wesley, 1980

2. C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation , San Francisco: Freeman, 1973