BalckHole_Trajectory.mws
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 EulerLagrange 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
.
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);

Applying the EulerLagrange 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);

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

We list three equations of motion.
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.
We solve these differential equations numerically.
> 
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, a0.01*a..a+0.01*a);

> 
ini := Eq77, Eq78, Eq79;

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

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*(x2))], 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):

We animate the motion by displaying plot frames in sequence.
> 
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: AddisonWesley, 1980
2. C. W. Misner, K. S. Thorne, and J. A. Wheeler,
Gravitation
, San Francisco: Freeman, 1973