Application Center - Maplesoft

App Preview:

Trajectory Near a Black Hole: an application of Lagrangian mechanics

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



  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 M  is the Schwarzschild metric,

ds^2 = -(1-2*M/r)*dt^2+(1-2*M/r)^(-1)*dr^2+r^2*(d*theta^2+sin(theta)^2*d*phi^2)

which indicates a curved space.

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

delta*int(1,s) = 0  

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

ds^2 = -(1-2*M/r)*dt^2+(1-2*M/r)^(-1)*dr^2+r^2*d*phi^2  

We let s  be the independent variable; t(s) , r(s)  and phi(s)  become the dependent variables.  Therefore,

ds^2 = [-(1-2*M/r)*(dt/ds)^2+(1-2*M/r)^(-1)*(dr/ds)^2+r^2*(d*phi/ds)^2] ds^2

Identifying the function f  as

2*f = -(1-2*M/r)*(dt/ds)^2+(1-2*M/r)^(-1)*(dr/ds)^2+r^2*(d*phi/ds)^2 .

because it is equivalent to find the minimum of 1/2 ds^2  as ds .  

>    restart:

Defining f

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

f := -1/2*(1-2*M/r(s))*diff(t(s),s)^2+1/2*1/(1-2*M/r(s))*diff(r(s),s)^2+1/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);

f1 := -1/2*(1-2*M/var3)*var2^2+1/2/(1-2*M/var3)*var4^2+1/2*var3^2*var6^2

we see that there is no explicit t(s)  term; there is a constant, which we denote as -a .

>    Eq11 := diff(f1, var2);

Eq11 := -(1-2*M/var3)*var2

>    Eq12 := diff(f1, var1);

Eq12 := 0

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

Eq13 := -(1-2*M/r(s))*diff(t(s),s)

>    Eq14 := Eq13 = -a;

Eq14 := -(1-2*M/r(s))*diff(t(s),s) = -a

Applying the Euler-Lagrange equation to the r  coordinate yields

>    Eq21 := diff(f1, var4);

Eq21 := 1/(1-2*M/var3)*var4

>    Eq22 := diff(f1, var3);

Eq22 := -M/var3^2*var2^2-1/(1-2*M/var3)^2*var4^2*M/var3^2+var3*var6^2

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

Eq23 := 1/(1-2*M/r(s))*diff(r(s),s)

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

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

>    Eq25 := diff(Eq23, s);

Eq25 := -2/(1-2*M/r(s))^2*diff(r(s),s)^2*M/r(s)^2+1/(1-2*M/r(s))*diff(r(s),`$`(s,2))

>    Eq26 := Eq25 - Eq24 = 0;

Eq26 := -1/(1-2*M/r(s))^2*diff(r(s),s)^2*M/r(s)^2+1/(1-2*M/r(s))*diff(r(s),`$`(s,2))+M/r(s)^2*diff(t(s),s)^2-r(s)*diff(phi(s),s)^2 = 0

Also there is no explicit phi(t)  term, so we denote another constant h

>    Eq31 := diff(f1, var6);

Eq31 := var3^2*var6

>    Eq32 := diff(f1, var5);

Eq32 := 0

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

Eq33 := r(s)^2*diff(phi(s),s)

>    Eq34 := Eq33 = h;

Eq34 := r(s)^2*diff(phi(s),s) = h

We list three equations of motion.

>    Eq14;

-(1-2*M/r(s))*diff(t(s),s) = -a

>    Eq26;

-1/(1-2*M/r(s))^2*diff(r(s),s)^2*M/r(s)^2+1/(1-2*M/r(s))*diff(r(s),`$`(s,2))+M/r(s)^2*diff(t(s),s)^2-r(s)*diff(phi(s),s)^2 = 0

>    Eq34;

r(s)^2*diff(phi(s),s) = h

We make some rearrangements to decouple the equations.

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

Eq53 := diff(t(s),s) = -a/(-1+2*M/r(s))

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

Eq54 := diff(phi(s),s) = h/r(s)^2

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

Eq55 := -1/(1-2*M/r(s))^2*diff(r(s),s)^2*M/r(s)^2+1/(1-2*M/r(s))*diff(r(s),`$`(s,2))+M/r(s)^2*a^2/(-1+2*M/r(s))^2-1/r(s)^3*h^2 = 0

We list the r   and phi   equations after rearrangement.

>    Eq54;

diff(phi(s),s) = h/r(s)^2

>    Eq55;

-1/(1-2*M/r(s))^2*diff(r(s),s)^2*M/r(s)^2+1/(1-2*M/r(s))*diff(r(s),`$`(s,2))+M/r(s)^2*a^2/(-1+2*M/r(s))^2-1/r(s)^3*h^2 = 0

We solve these differential equations numerically.

>    M := 1;

M := 1

>    Eq77 := r(0) = 11;

Eq77 := r(0) = 11

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

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

>    Eq79 := phi(0) = 0;

Eq79 := phi(0) = 0

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

Eq80 := D(phi)(0) = .295e-1

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

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

h := 3.5695

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

a := .9509661236

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

[Maple Plot]

>    ini := Eq77, Eq78, Eq79;

ini := r(0) = 11, D(r)(0) = 0, phi(0) = 0

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

Eq91 := [s = proc (s) option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; s end proc, phi(s) = proc (s) local res, solnproc, outpoint, ndsol; option `Copyright (c) 2000 by ...

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

[Maple Plot]

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 z  to accommodate the actual length.  The embedding formula for this problem is

z = sqrt(8*M*(r-2*M))

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

[Maple Plot]

We animate the motion by displaying plot frames in sequence.

>    N := 50; nstep := 10;

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,

[Maple Plot]

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

[Maple Plot]

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




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