Relativistic Orbits with Maple

Frank Wang, fw@phys.columbia.edu

Using Maple  to derive equations of motion from the Lagrangian, to solve differential equations numerically, and to form graphs based on numerical solutions, we present trajectories of a particle under gravity based on Newtonian theory and general relativity with the Schwarzschild metric and the Kerr metric.  

The Lagrangian for the Kerr solution in general relativity, which describes a spherical object of mass M    and spin angular momentum J   , is

L = -(1-2*M/r)*diff(t(tau),tau)^2/2-2*a*M/r*diff(t(tau),tau)*diff(phi(tau),tau)+r^2/(2*Delta)*diff(r(tau),tau)^2+(r^2+a^2+2*M*a^2/r)*diff(phi(tau),tau)^2/2    

where

a = J/M    and Delta = r^2-2*M*r+a^2   .

If a = 0   , namely the object does not rotate, the Kerr solution becomes the Schwarzschild solution.  

In the limit of c = infinity    and t = tau   , the Schwarzschild solution further reduces to the Lagrangian for Newtonian mechanics, i.e., the difference between kinetic energy and potential energy.  

To find equations of motion, we simply employ the Euler-Lagrange equation

d/dt     diff(L,`q'`[i])     -   diff(L,q[i])     =  0

where q    is a coordinate and q*`'`    denotes its derivative with respect to time:   q*`'`    =   dq/dt  (conventional notation for this quantity as q    with a dot above this letter).

For additional information, see F. Wang, "Relativistic orbits with computer algebra," Am . J . Phys . 72 , 1040--1044 (2004), and references therein.  Contact the author for reprint.     

By varying the initial conditions, one can discover many interesting orbits in relativity.  I thank Professor Roman Koniuk (koniuk@yorku.ca) for suggesting the last three sections.  

>    restart: with(plots): with(plottools): 

Warning, the name changecoords has been redefined




Warning, the name arrow has been redefined

Newtonian mechanics

First define the kinetic energy, the potential energy, and the Lagrangian:

>    T := 1/2*(diff(r(t),t)^2 + r(t)^2*diff(phi(t),t)^2);

T := 1/2*diff(r(t),t)^2+1/2*r(t)^2*diff(phi(t),t)^2

>    V := -G*M/r(t);

V := -G*M/r(t)

>    L := T - V;

L := 1/2*diff(r(t),t)^2+1/2*r(t)^2*diff(phi(t),t)^2+G*M/r(t)

>    L1 := subs({r(t)=var1, diff(r(t),t)=var2, phi(t)=var3, diff(phi(t),t)=var4}, L):

Employ the Euler-Lagrange equation for the Newtonian Lagrangian (Epr11 through Eq26)

>    Epr11 := diff(L1, var4):

>    Epr12 := diff(L1, var3):

>    Epr13 := subs({var1=r(t),var2=diff(r(t),t),var3=phi(t),var4=diff(phi(t),t)},Epr11):

>    Eq14 := Epr13 = l;

Eq14 := r(t)^2*diff(phi(t),t) = l

>    Epr21 := diff(L1, var2):

>    Epr22 := diff(L1, var1):

>    Epr23 := subs({var1=r(t), var2=diff(r(t),t), var3=phi(t), var4=diff(phi(t),t)}, Epr21):

>    Epr24 := subs({var1=r(t), var2=diff(r(t),t), var3=phi(t), var4=diff(phi(t),t)}, Epr22):

>    Epr25 := diff(Epr23,t):

>    Eq26 := Epr25 - Epr24 = 0;

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

Eq31 and Eq32 decouple the differential equations.

>    Eq31 := isolate(Eq14, diff(phi(t),t));

Eq31 := diff(phi(t),t) = l/r(t)^2

>    Eq32 := eval(Eq26, Eq31);

Eq32 := diff(r(t),`$`(t,2))-1/r(t)^3*l^2+G*M/r(t)^2 = 0

Provide the initial values, and use Maple  to solve these differential equations numerically and to form plots based on numerical solutions.  

>    G:=1; M:=1;

G := 1

M := 1

>    Eq41 := r(0) = 26;

Eq41 := r(0) = 26

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

Eq42 := D(r)(0) = 0

>    Eq43 := phi(0) = 0;

Eq43 := phi(0) = 0

>    Eq44 := D(phi)(0) = 0.0071;

Eq44 := D(phi)(0) = .71e-2

>    En := eval(T + V, {r(t)=rhs(Eq41), diff(r(t),t)=rhs(Eq42), diff(phi(t),t)=rhs(Eq44)});

En := -.2142295846e-1

>    l := eval(lhs(Eq14), {r(t)=rhs(Eq41), diff(phi(t),t)=rhs(Eq44)});

l := 4.7996

The eccentricity is not needed; it is listed as a reference.

>    epsilon := sqrt(1 + 2*En*l^2/((G*M)^2));  

epsilon := .1139938402

>    ini1 := Eq41, Eq42, Eq43;

ini1 := r(0) = 26, D(r)(0) = 0, phi(0) = 0

>    Eq51:=dsolve({Eq31, Eq32, ini1}, {r(t), phi(t)}, numeric, output=listprocedure);

Eq51 := [t = proc (t) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc := ...
Eq51 := [t = proc (t) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc := ...

>    polarplot([rhs(Eq51(t)[3]), rhs(Eq51(t)[2]), t=0..720], scaling=constrained, legend="Newtonian");

[Maple Plot]

>    pnewton := polarplot([rhs(Eq51(t)[3]), rhs(Eq51(t)[2]), t=0..720], scaling=constrained, color=black, legend="Newtonian", linestyle=3, thickness=2):

Schwarzschild metric

>    G := 'G': M := 'M':

The procedure for the Schwarzschild metric is almost identical to the Newtonian mechanics; we merely change the Lagrangian.

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

f := -1/2*(c^2-2*G*M/r(tau))*diff(t(tau),tau)^2+1/2*1/(1-2*G*M/c^2/r(tau))*diff(r(tau),tau)^2+1/2*r(tau)^2*diff(phi(tau),tau)^2

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

Employ the Euler-Lagrange equations (Epr11 through Eq34)

>    Epr11 := diff(f1, var2):

>    Epr12 := diff(f1, var1):

>    Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):

>    Eq14 := Epr13 = -a;

Eq14 := -(c^2-2*G*M/r(tau))*diff(t(tau),tau) = -a

>    Epr21 := diff(f1, var4):

>    Epr22 := diff(f1, var3):

>    Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):

>    Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):

>    Epr25 := diff(Epr23, tau):

>    Eq26 := Epr25 - Epr24 = 0;

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

>    Epr31 := diff(f1, var6):

>    Epr32 := diff(f1, var5):

>    Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):

>    Eq34 := Epr33 = h;

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

Decouple the differential equations in the following three lines.

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

Eq53 := diff(t(tau),tau) = -a/(-c^2+2*G*M/r(tau))

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

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

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

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

Provide the initial conditions, and plot the numerical solutions.

>    M := 1; G := 1; c := 1;

M := 1

G := 1

c := 1

>    Eq77 := r(0) = 26;

Eq77 := r(0) = 26

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

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

>    Eq79 := phi(0) = 0;

Eq79 := phi(0) = 0

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

Eq80 := D(phi)(0) = .71e-2

>    h := eval(lhs(Eq34), {r(tau)=rhs(Eq77), diff(phi(tau),tau)=rhs(Eq80)});

h := 4.7996

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

a := .9770019258

The effective potential illuminates the turning points.  

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

>    plot([Vsq, a], r=2..40, 0.973..1.12, legend=["effective potential", "energy"]);

[Maple Plot]

>    fsolve(Vsq=a, r, 10..20);

15.46812439

>    fsolve(Vsq=a, r, 20..30);

25.99999991

>    ini := Eq77, Eq78, Eq79;

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

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

Eq91 := [tau = proc (tau) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc...
Eq91 := [tau = proc (tau) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc...

>    psch1 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2300],
          scaling=constrained, axesfont=[TIMES, ROMAN, 12], legend="Schwarzschild"):

>    psch2 := disk([15.47*cos(3.7797), 15.47*sin(3.7797)], 0.5, color=black):

>    psch3 := disk([15.47*cos(11.3399), 15.47*sin(11.3399)], 0.5, color=black):

>    psch4 := disk([15.47*cos(18.9008), 15.47*sin(18.9008)], 0.5, color=black):

>    psch5 := disk([15.47*cos(26.4221), 15.47*sin(26.4221)], 0.5, color=black):

>    pns := disk([0,0],5.8, color=orange):

>    display([psch1, psch2, psch3, psch4, psch5, pns, pnewton]);

[Maple Plot]

>    pschwarz := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000],
          scaling=constrained, color=blue, legend=Schwarzschild):

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

Kerr metric

>    M := 'M': l := 'l':

The same procedure also applies to the Kerr metric; we only change the Lagrangian:

>    a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2;

a := J/M

Delta := r(tau)^2-2*M*r(tau)+J^2/M^2

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

f := -1/2*(1-2*M/r(tau))*diff(t(tau),tau)^2+1/2*r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2+1/2*(r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)^2-2*J/r(tau)*diff(t(tau),tau)*diff(phi...
f := -1/2*(1-2*M/r(tau))*diff(t(tau),tau)^2+1/2*r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2+1/2*(r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)^2-2*J/r(tau)*diff(t(tau),tau)*diff(phi...

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

Employ the Euler-Lagrange equations (Epr11 through Eq34)

>    Epr11 := diff(f1, var2):

>    Epr12 := diff(f1, var1):

>    Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):

>    Eq14 := Epr13 = -E;

Eq14 := -(1-2*M/r(tau))*diff(t(tau),tau)-2*J/r(tau)*diff(phi(tau),tau) = -E

>    Epr21 := diff(f1, var4):

>    Epr22 := diff(f1, var3):

>    Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):

>    Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):

>    Epr25 := diff(Epr23, tau):

>    Eq26 := Epr25 - Epr24 = 0;

Eq26 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq26 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq26 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...

>    Epr31 := diff(f1, var6):

>    Epr32 := diff(f1, var5):

>    Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):

>    Eq34 := Epr33 = l;

Eq34 := (r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)-2*J/r(tau)*diff(t(tau),tau) = l

A similar decoupling manipulation (Eq53 through Eq55)

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

Eq53 := diff(t(tau),tau) = (-E+2*J/r(tau)*diff(phi(tau),tau))/(-1+2*M/r(tau))

>    Eq53p := subs(Eq53, Eq34);

Eq53p := (r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)-2*J/r(tau)*(-E+2*J/r(tau)*diff(phi(tau),tau))/(-1+2*M/r(tau)) = l

>    Eq54 := isolate(Eq53p, diff(phi(tau),tau));

Eq54 := diff(phi(tau),tau) = (l*M^2*(-r(tau)+2*M)-2*J*M^2*E)/(-r(tau)^3*M^2+2*r(tau)^2*M^3-J^2*r(tau))

>    Eq54p := subs(Eq54, Eq53);

Eq54p := diff(t(tau),tau) = (-E+2*J/r(tau)*(l*M^2*(-r(tau)+2*M)-2*J*M^2*E)/(-r(tau)^3*M^2+2*r(tau)^2*M^3-J^2*r(tau)))/(-1+2*M/r(tau))

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

Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...

Solve the differential equations with initial conditions provided:

>    M := 1; J := 0.37;

M := 1

J := .37

>    l := 4.7996; E := 0.9772;

l := 4.7996

E := .9772

>    Eq77 := r(0) = 26;

Eq77 := r(0) = 26

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

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

>    Eq79 := phi(0) = 0;

Eq79 := phi(0) = 0

>    Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)});

Eq80 := D(phi)(0) = .7143004388e-2

>    ini := Eq77, Eq78, Eq79;

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

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

Eq91 := [tau = proc (tau) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc...
Eq91 := [tau = proc (tau) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc...

>    polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2500],
          scaling=constrained, color=red, legend="Kerr direct");

[Maple Plot]

>    pkerrd := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000],
          scaling=constrained, color=red, legend="Kerr direct"):

>    M := 1; J := -0.37;

M := 1

J := -.37

>    l := 4.7996; E := 0.976797;

l := 4.7996

E := .976797

>    Eq77 := r(0) = 26;

Eq77 := r(0) = 26

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

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

>    Eq79 := phi(0) = 0;

Eq79 := phi(0) = 0

>    Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)});

Eq80 := D(phi)(0) = .7053899319e-2

>    ini := Eq77, Eq78, Eq79;

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

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

Eq91 := [tau = proc (tau) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc...
Eq91 := [tau = proc (tau) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc...

>    polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2500],
          scaling=constrained, legend="Kerr retro");

[Maple Plot]

>    pkerro := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000],
          scaling=constrained, color=green, legend="Kerr retro", axesfont=[TIMES, ROMAN, 12]):

>    pkerr1 := disk([16.31*cos(3.63467), 16.31*sin(3.63467)], 0.5, color=black):

>    pkerr2 := disk([14.53*cos(3.923), 14.53*sin(3.923)], 0.5, color=black):

>    pkerr3 := disk([16.31*cos(10.9765), 16.31*sin(10.9765)], 0.5, color=black):

>    pkerr4 := disk([14.53*cos(11.7778), 14.53*sin(11.778)], 0.5, color=black):

>    display([pnewton,pschwarz,pkerrd,pkerro,pkerr1,pkerr2,psch2,pkerr3,pkerr4,psch3,pns]);

[Maple Plot]

Embedding diagram

We derive the embedding formula, see Misner, Thorne, and Wheeler, Gravitation , p. 613 for details.  A star is assumed to have uniform density:

>    M := 'M': m := rp^3/R^3*M;

m := rp^3/R^3*M

>    zin := int(sqrt(1/(1 - 2*m/rp) -1), rp=2*M..r) + C;

zin := int((1/(1-2*rp^2/R^3*M)-1)^(1/2),rp = 2*M .. r)+C

>    zout := int(sqrt(1/(1 - 2*M/rp) -1), rp=2*M..r);

zout := 2*2^(1/2)*(M/(r-2*M))^(1/2)*(r-2*M)

>    Eq201 := eval(zin, r=R) = eval(zout, r=R):

>    C := solve(Eq201, C):

>    pin := plot3d([r*cos(phi), r*sin(phi), eval(zin, {M=1, R=5.8})], phi=0..2*Pi, r=0..5.8, grid=[25,6], orientation=[88,135], shading=zhue):

>    pout := plot3d([r*cos(phi), r*sin(phi), eval(zout, {M=1, R=5.8})], phi=0..2*Pi, r=5.8..27, grid=[25,12]):

We replot the trajectory in the section "Schwarzschild metric" in this embedding diagram:

>    display([pin, pout, pschwarz3d], style=hidden);

[Maple Plot]

>   

Example 1

This example demonstrates an infalling particle with zero angular momentum encountering a Kerr black hole with maximal spin angular momentum J  ( a = M ).

>    restart:

>    a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2;

a := J/M

Delta := r(tau)^2-2*M*r(tau)+J^2/M^2

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

f := -1/2*(1-2*M/r(tau))*diff(t(tau),tau)^2+1/2*r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2+1/2*(r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)^2-2*J/r(tau)*diff(t(tau),tau)*diff(phi...
f := -1/2*(1-2*M/r(tau))*diff(t(tau),tau)^2+1/2*r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2+1/2*(r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)^2-2*J/r(tau)*diff(t(tau),tau)*diff(phi...

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

>    Epr11 := diff(f1, var2):

>    Epr12 := diff(f1, var1):

>    Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):

>    Eq14 := Epr13 = -E;

Eq14 := -(1-2*M/r(tau))*diff(t(tau),tau)-2*J/r(tau)*diff(phi(tau),tau) = -E

>    Epr21 := diff(f1, var4):

>    Epr22 := diff(f1, var3):

>    Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):

>    Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):

>    Epr25 := diff(Epr23, tau):

>    Eq26 := Epr25 - Epr24 = 0;

Eq26 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq26 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq26 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...

>    Epr31 := diff(f1, var6):

>    Epr32 := diff(f1, var5):

>    Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):

>    Eq34 := Epr33 = l;

Eq34 := (r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)-2*J/r(tau)*diff(t(tau),tau) = l

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

Eq53 := diff(t(tau),tau) = (-E+2*J/r(tau)*diff(phi(tau),tau))/(-1+2*M/r(tau))

>    Eq53p := subs(Eq53, Eq34);

Eq53p := (r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)-2*J/r(tau)*(-E+2*J/r(tau)*diff(phi(tau),tau))/(-1+2*M/r(tau)) = l

>    Eq54 := isolate(Eq53p, diff(phi(tau),tau));

Eq54 := diff(phi(tau),tau) = (-l*M^2*(-r(tau)+2*M)+2*J*M^2*E)/(r(tau)^3*M^2-2*r(tau)^2*M^3+J^2*r(tau))

>    Eq54p := subs(Eq54, Eq53);

Eq54p := diff(t(tau),tau) = (-E+2*J/r(tau)*(-l*M^2*(-r(tau)+2*M)+2*J*M^2*E)/(r(tau)^3*M^2-2*r(tau)^2*M^3+J^2*r(tau)))/(-1+2*M/r(tau))

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

Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...

>    M := 1; J := 1;

M := 1

J := 1

>    l := 0; E := 1;

l := 0

E := 1

>    Eq77 := r(0) = 5;

Eq77 := r(0) = 5

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

Eq78 := D(r)(0) = -.64

>    Eq79 := phi(0) = 0;

Eq79 := phi(0) = 0

>    Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)});

Eq80 := D(phi)(0) = 1/40

>    ini := Eq77, Eq78, Eq79;

ini := r(0) = 5, D(r)(0) = -.64, phi(0) = 0

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

Eq91 := [tau = proc (tau) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc...
Eq91 := [tau = proc (tau) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc...

>    with(plots): 

Warning, the name changecoords has been redefined

>    polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..4.49],
          scaling=constrained, numpoints=2400);

[Maple Plot]

>   

Example 2

This example demonstrates the behavior of a particle near the unstable circular orbit.  The particle orbits several times, then flies out.  An animation of the motion is included.    

>    restart:

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

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

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

>    Epr11 := diff(f1, var2):

>    Epr12 := diff(f1, var1):

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

>    Eq14 := Epr13 = -a/c;

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

>    Epr21 := diff(f1, var4):

>    Epr22 := diff(f1, var3):

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

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

>    Epr25 := diff(Epr23, s):

>    Eq26 := Epr25 - Epr24 = 0;

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

>    Epr31 := diff(f1, var6):

>    Epr32 := diff(f1, var5):

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

>    Eq34 := Epr33 = h/c;

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

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

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

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

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

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

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

>    M := 1; G := 1; c := 1;

M := 1

G := 1

c := 1

>    Eq77 := r(0) = 3.76777;

Eq77 := r(0) = 3.76777

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

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

>    Eq79 := phi(0) = 0;

Eq79 := phi(0) = 0

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

Eq80 := D(phi)(0) = .30290054

>    h := eval(lhs(Eq34), {r(s)=rhs(Eq77), diff(phi(s),s)=rhs(Eq80)});

h := 4.300003560

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

a := 1.039364786

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

>    plot([Vsq, a], x=2..25, a-0.1*a..a+0.1*a, legend=["effective potential", "energy"]);

[Maple Plot]

>    ini := Eq77, Eq78, Eq79;

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

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

Eq91 := [s = proc (s) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc := ...
Eq91 := [s = proc (s) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; end module; solnproc := ...

>    with(plots): 

Warning, the name changecoords has been redefined

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

[Maple Plot]

>    p1 := plot3d([x*cos(t), x*sin(t), sqrt(8*(x-2))], x=2..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..100], coords=cylindrical, color=black, numpoints=400):

>    display([p1, p2]);

[Maple Plot]

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

[Maple Plot]

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

>    display([p1, p3]);

[Maple Plot]

Example 3

Literature on numerical integrations for relativistic orbits can be found in the early 1970's; see M. Johnston and R. Ruffini, "Generalized Wilkins effect and selected orbits in a Kerr-Newman geometry," Phys. Rev.  D 10 , 2324-2329 (1974).  This example uses the same initial conditions for Fig. 6 of that paper (a particle initially orbits against a maximally rotating black hole).  Our calculation shows interesting details of the last orbits that was omitted in that paper perhaps due to limited computational power of that era.      

>    restart:

>    a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2;

a := J/M

Delta := r(tau)^2-2*M*r(tau)+J^2/M^2

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

f := -1/2*(1-2*M/r(tau))*diff(t(tau),tau)^2+1/2*r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2+1/2*(r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)^2-2*J/r(tau)*diff(t(tau),tau)*diff(phi...
f := -1/2*(1-2*M/r(tau))*diff(t(tau),tau)^2+1/2*r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2+1/2*(r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)^2-2*J/r(tau)*diff(t(tau),tau)*diff(phi...

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

>    Epr11 := diff(f1, var2):

>    Epr12 := diff(f1, var1):

>    Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):

>    Eq14 := Epr13 = -E;

Eq14 := -(1-2*M/r(tau))*diff(t(tau),tau)-2*J/r(tau)*diff(phi(tau),tau) = -E

>    Epr21 := diff(f1, var4):

>    Epr22 := diff(f1, var3):

>    Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):

>    Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):

>    Epr25 := diff(Epr23, tau):

>    Eq26 := Epr25 - Epr24 = 0;

Eq26 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq26 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq26 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...

>    Epr31 := diff(f1, var6):

>    Epr32 := diff(f1, var5):

>    Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
              var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):

>    Eq34 := Epr33 = l;

Eq34 := (r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)-2*J/r(tau)*diff(t(tau),tau) = l

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

Eq53 := diff(t(tau),tau) = (-E+2*J/r(tau)*diff(phi(tau),tau))/(-1+2*M/r(tau))

>    Eq53p := subs(Eq53, Eq34);

Eq53p := (r(tau)^2+J^2/M^2+2/M*J^2/r(tau))*diff(phi(tau),tau)-2*J/r(tau)*(-E+2*J/r(tau)*diff(phi(tau),tau))/(-1+2*M/r(tau)) = l

>    Eq54 := isolate(Eq53p, diff(phi(tau),tau));

Eq54 := diff(phi(tau),tau) = (l*M^2*(-r(tau)+2*M)-2*J*M^2*E)/(-r(tau)^3*M^2+2*r(tau)^2*M^3-J^2*r(tau))

>    Eq54p := subs(Eq54, Eq53);

Eq54p := diff(t(tau),tau) = (-E+2*J/r(tau)*(l*M^2*(-r(tau)+2*M)-2*J*M^2*E)/(-r(tau)^3*M^2+2*r(tau)^2*M^3-J^2*r(tau)))/(-1+2*M/r(tau))

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

Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...
Eq55 := r(tau)/(r(tau)^2-2*M*r(tau)+J^2/M^2)*diff(r(tau),tau)^2-r(tau)^2/(r(tau)^2-2*M*r(tau)+J^2/M^2)^2*diff(r(tau),tau)*(2*r(tau)*diff(r(tau),tau)-2*M*diff(r(tau),tau))+r(tau)^2/(r(tau)^2-2*M*r(tau)+...

>    M := 1; J := -1;

M := 1

J := -1

>    l := 4; E := 0.97;

l := 4

E := .97

>    A := (r^3 + a^2*r + 2*M*a^2);

A := r^3+r+2

>    B := -4*a*M*l;

B := 16

>    C := -(r - 2*M)*l^2 - r*(r^2 - 2*M*r + a^2);

C := -16*r+32-r*(r^2-2*r+1)

>    Veff := (-B + sqrt(B^2 - 4*A*C))/(2*A);

Veff := 1/2*(-16+2*(18*r^4-32*r^3+r^6-2*r^5+13*r^2+2*r)^(1/2))/(r^3+r+2)

>    plot({Veff, E}, r=1..30, 0..1.2);

[Maple Plot]

>    fsolve(Veff=E, r=20..30);

23.95403154

>    Eq77 := r(0) = %;

Eq77 := r(0) = 23.95403154

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

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

>    Eq79 := phi(0) = -5*Pi/4;

Eq79 := phi(0) = -5/4*Pi

>    Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)});

Eq80 := D(phi)(0) = .6804181347e-2

>    ini := Eq77, Eq78, Eq79;

ini := r(0) = 23.95403154, D(r)(0) = 0, phi(0) = -5/4*Pi

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

Eq91 := [tau = proc (tau) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; Data := table([(
Eq91 := [tau = proc (tau) local res, data, solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; data := module () local Data; export Get, Set; Data := table([(

>    with(plots): with(plottools): 

Warning, the name changecoords has been redefined




Warning, the name arrow has been redefined

>    pbh := disk([0,0], 1, color=black):

>    pergo := polarplot(2, phi=0..2*Pi, color=black, thickness=2):

>    porb1 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=-180..207.4], scaling=constrained, axes=boxed):

This plot resembles Fig. 6 of Johnston and Ruffini: the particle initially orbits counterclockwise, and the black hole spins clockwise.  There is an indication of reversal of direction when the particle approaches the black hole, but we can see more details in the following.  

>    display([pbh, porb1]);

[Maple Plot]

>    pergo3d := sphereplot(M + sqrt(M^2 - a^2*cos(theta)^2), phi=0..2*Pi, theta=0..Pi, scaling=constrained, style=wireframe, grid=[16,16]):

>    pkerrbh := sphere([0,0,0], 1, color=black):

This plot shows the ergosphere of a Kerr black hole: the region between the static limit described by r = M+sqrt(M^2-a^2*cos(theta)^2)  and the event horizon described by r = M .  Within the ergosphere, a particle (even a photon) must corotate with the black hole.  

>    display([pergo3d, pkerrbh]);

[Maple Plot]

>    porb3d1 := spacecurve([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), 0, tau=200..207.4], coords=cylindrical, scaling=constrained, color=red, thickness=2):

>    porb3d2 := spacecurve([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), 0, tau=207.4..207.648], coords=cylindrical, scaling=constrained, numpoints=800, color=red, thickness=2):

>    porb2 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=200..207.4], scaling=constrained, axes=boxed):

>    porb3 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=207.4..207.648], scaling=constrained, numpoints=400, axes=boxed):

The numerical solution shows that the particle reverses its direction and corotates with the black hole as it spirals toward it.  

>    display([pkerrbh, pergo3d, porb3d1, porb3d2]);

[Maple Plot]

Viewing from the top:

>    display([pbh, pergo, porb2, porb3]);

[Maple Plot]

>   

>