Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

SageMath notebooks associated to the Black Hole Lectures (https://luth.obspm.fr/~luthier/gourgoulhon/bh16)

Project: BHLectures
Views: 20113
Kernel: SageMath 10.0.rc0

Solving the ODE for outgoing radial null geodesics in Vaidya spacetime

This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes.

The computations make use of tools developed through the SageManifolds project.

version()
'SageMath version 10.0.rc0, Release Date: 2023-04-23'
%display latex

x = var('x', domain='real')

a = var('a', latex_name=r'\alpha', domain='real')

f(x) = (1 - a*x) / (a*x^2 - x + 2)
f(x)

−αx−1αx2−x+2\displaystyle -\frac{{\alpha} x - 1}{{\alpha} x^{2} - x + 2}

Case of a high density radiation shell

assume(a > 1/8)

F(x) = integrate(f(x), x)
F(x)

arctan⁡(2 αx−18 α−1)8 α−1−12 log⁡(αx2−x+2)\displaystyle \frac{\arctan\left(\frac{2 \, {\alpha} x - 1}{\sqrt{8 \, {\alpha} - 1}}\right)}{\sqrt{8 \, {\alpha} - 1}} - \frac{1}{2} \, \log\left({\alpha} x^{2} - x + 2\right)

P(x) = a*x^2 - x + 2

x1 = 1/(2*a)

P(x1)

−14 α+2\displaystyle -\frac{1}{4 \, {\alpha}} + 2

F(x1)

−12 log⁡(−14 α+2)\displaystyle -\frac{1}{2} \, \log\left(-\frac{1}{4 \, {\alpha}} + 2\right)

exp(F(x1)).simplify_full()

28 α−1α\displaystyle \frac{2}{\sqrt{\frac{8 \, {\alpha} - 1}{{\alpha}}}}

P(1/a)

2\displaystyle 2

thb = -4*exp(-2/sqrt(16/x - 1)*atan(1/sqrt(16/x - 1)))
thb

−4 e(−2 arctan⁡(116x−1)16x−1)\displaystyle -4 \, e^{\left(-\frac{2 \, \arctan\left(\frac{1}{\sqrt{\frac{16}{x} - 1}}\right)}{\sqrt{\frac{16}{x} - 1}}\right)}

plot(thb, (x, 0.001, 15.99))
Image in a Jupyter notebook
taylor(thb, x, 0, 3)

−13840 x3−196 x2+12 x−4\displaystyle -\frac{1}{3840} \, x^{3} - \frac{1}{96} \, x^{2} + \frac{1}{2} \, x - 4

x

x\displaystyle x

f(a, x) = sqrt(2)/sqrt(a*x^2 - x + 2)*exp(1/sqrt(8*a - 1)*\
          (atan((2*a*x - 1)/sqrt(8*a - 1)) + atan(1/sqrt(8*a - 1))))
f

(α,x) ↦ 2e(arctan⁡(2 αx−18 α−1)+arctan⁡(18 α−1)8 α−1)αx2−x+2\displaystyle \left( {\alpha}, x \right) \ {\mapsto} \ \frac{\sqrt{2} e^{\left(\frac{\arctan\left(\frac{2 \, {\alpha} x - 1}{\sqrt{8 \, {\alpha} - 1}}\right) + \arctan\left(\frac{1}{\sqrt{8 \, {\alpha} - 1}}\right)}{\sqrt{8 \, {\alpha} - 1}}\right)}}{\sqrt{{\alpha} x^{2} - x + 2}}

diff(f(a, x), x)

−2(2 αx−1)e(arctan⁡(2 αx−18 α−1)+arctan⁡(18 α−1)8 α−1)2 (αx2−x+2)32+2 2αe(arctan⁡(2 αx−18 α−1)+arctan⁡(18 α−1)8 α−1)αx2−x+2((2 αx−1)28 α−1+1)(8 α−1)\displaystyle -\frac{\sqrt{2} {\left(2 \, {\alpha} x - 1\right)} e^{\left(\frac{\arctan\left(\frac{2 \, {\alpha} x - 1}{\sqrt{8 \, {\alpha} - 1}}\right) + \arctan\left(\frac{1}{\sqrt{8 \, {\alpha} - 1}}\right)}{\sqrt{8 \, {\alpha} - 1}}\right)}}{2 \, {\left({\alpha} x^{2} - x + 2\right)}^{\frac{3}{2}}} + \frac{2 \, \sqrt{2} {\alpha} e^{\left(\frac{\arctan\left(\frac{2 \, {\alpha} x - 1}{\sqrt{8 \, {\alpha} - 1}}\right) + \arctan\left(\frac{1}{\sqrt{8 \, {\alpha} - 1}}\right)}{\sqrt{8 \, {\alpha} - 1}}\right)}}{\sqrt{{\alpha} x^{2} - x + 2} {\left(\frac{{\left(2 \, {\alpha} x - 1\right)}^{2}}{8 \, {\alpha} - 1} + 1\right)} {\left(8 \, {\alpha} - 1\right)}}

a_sel = [1/4, 1/2, 1, 2, 4, 6]
graph = Graphics()
for a1 in a_sel:
    legend_label=r'$\alpha = {:.2f}$'.format(float(a1))
    graph += plot(f(a1, x), (x, 0, 6), color=hue(a1/8), 
                  legend_label=legend_label, thickness=1.5,
                  axes_labels=[r"$x$", r"$f_\alpha(x)$"])
show(graph, frame=True, axes=False, gridlines=True)
graph.save("vai_f_alpha_x.pdf", frame=True, axes=False, gridlines=True)
Image in a Jupyter notebook

Equivalent form:

h(a, x) = sqrt(2)/sqrt(a*x^2 - x + 2)*exp(1/sqrt(8*a - 1)*\
           (atan(x*sqrt(8*a - 1)/(4 - x)) + pi*heaviside(x-4)))
h

(α,x) ↦ 2e(πH(x−4)+arctan⁡(−8 α−1xx−4)8 α−1)αx2−x+2\displaystyle \left( {\alpha}, x \right) \ {\mapsto} \ \frac{\sqrt{2} e^{\left(\frac{\pi H\left(x - 4\right) + \arctan\left(-\frac{\sqrt{8 \, {\alpha} - 1} x}{x - 4}\right)}{\sqrt{8 \, {\alpha} - 1}}\right)}}{\sqrt{{\alpha} x^{2} - x + 2}}

a_sel = [1/4, 1/2, 1, 2, 4, 6]
graph = Graphics()
for a1 in a_sel:
    legend_label=r'$\alpha = {:.2f}$'.format(float(a1))
    graph += plot(h(a1, x), (x, 0, 6), color=hue(a1/8), 
                  legend_label=legend_label, thickness=1.5,
                  axes_labels=[r"$x$", r"$f_\alpha(x)$"])
show(graph, frame=True, axes=False, gridlines=True)
Image in a Jupyter notebook

Maximal value of fα(x)f_\alpha(x):

f(a, 1/a)

e(2 arctan⁡(18 α−1)8 α−1)\displaystyle e^{\left(\frac{2 \, \arctan\left(\frac{1}{\sqrt{8 \, {\alpha} - 1}}\right)}{\sqrt{8 \, {\alpha} - 1}}\right)}

Expansion for small xx:

taylor(f(a, x), x, 0, 2)

−14 (α−1)x2+12 x+1\displaystyle -\frac{1}{4} \, {\left({\alpha} - 1\right)} x^{2} + \frac{1}{2} \, x + 1

Case of a low density radiation shell

forget(a > 1/8)

assume(a < 1/8)

f(x) = (1 - a*x) / (a*x^2 - x + 2)
F(x) = integrate(f(x), x)
F(x)

log⁡(2 αx−−8 α+1−12 αx+−8 α+1−1)2 −8 α+1−12 log⁡(αx2−x+2)\displaystyle \frac{\log\left(\frac{2 \, {\alpha} x - \sqrt{-8 \, {\alpha} + 1} - 1}{2 \, {\alpha} x + \sqrt{-8 \, {\alpha} + 1} - 1}\right)}{2 \, \sqrt{-8 \, {\alpha} + 1}} - \frac{1}{2} \, \log\left({\alpha} x^{2} - x + 2\right)

x01 = (1 - sqrt(1 - 8*a))/(2*a)
x02 = (1 + sqrt(1 - 8*a))/(2*a)
x01, x02

(−−8 α+1−12 α,−8 α+1+12 α)\displaystyle \left(-\frac{\sqrt{-8 \, {\alpha} + 1} - 1}{2 \, {\alpha}}, \frac{\sqrt{-8 \, {\alpha} + 1} + 1}{2 \, {\alpha}}\right)

(x01 + x02).simplify_full()

1α\displaystyle \frac{1}{{\alpha}}

bool(P(x) == a*(x - x01)*(x - x02))

True\displaystyle \mathrm{True}

x1 = var('x_1', domain='real')
x2 = var('x_2', domain='real')
a1 = 1/(x1 + x2)
a1

1x1+x2\displaystyle \frac{1}{x_{1} + x_{2}}

f1(x) = (1 - a1*x) / (a1*(x - x1)*(x - x2))
f1(x)

−(x1+x2)(xx1+x2−1)(x−x1)(x−x2)\displaystyle -\frac{{\left(x_{1} + x_{2}\right)} {\left(\frac{x}{x_{1} + x_{2}} - 1\right)}}{{\left(x - x_{1}\right)} {\left(x - x_{2}\right)}}

integrate(f1(x), x).simplify_full()

x2log⁡(x−x1)−x1log⁡(x−x2)x1−x2\displaystyle \frac{x_{2} \log\left(x - x_{1}\right) - x_{1} \log\left(x - x_{2}\right)}{x_{1} - x_{2}}

(x01 - x02).simplify_full()

−−8 α+1α\displaystyle -\frac{\sqrt{-8 \, {\alpha} + 1}}{{\alpha}}

g = x2/(x1 - x2)*ln(abs(x - x1)) - x1/(x1 - x2)*ln(abs(x - x2))
g

x2log⁡(∣x−x1∣)x1−x2−x1log⁡(∣x−x2∣)x1−x2\displaystyle \frac{x_{2} \log\left({\left| x - x_{1} \right|}\right)}{x_{1} - x_{2}} - \frac{x_{1} \log\left({\left| x - x_{2} \right|}\right)}{x_{1} - x_{2}}

diff(g, x)

−x1(x−x2)(x1−x2)+x2(x−x1)(x1−x2)\displaystyle -\frac{x_{1}}{{\left(x - x_{2}\right)} {\left(x_{1} - x_{2}\right)}} + \frac{x_{2}}{{\left(x - x_{1}\right)} {\left(x_{1} - x_{2}\right)}}

X1 = x01.subs(a=1/9).simplify_full()
X2 = x02.subs(a=1/9).simplify_full()
X1, X2

(3,6)\displaystyle \left(3, 6\right)

n(X1), n(X2)

(3.00000000000000,6.00000000000000)\displaystyle \left(3.00000000000000, 6.00000000000000\right)

b1 = (X1/(X2 - X1)).expand()
b2 = (X2/(X2 - X1)).expand()
b1, b2

(1,2)\displaystyle \left(1, 2\right)

R(x) = abs(x/X2 - 1)^b1/abs(x/X1 - 1)^b2
R(x)

9 ∣16 x−1∣(x−3)2\displaystyle \frac{9 \, {\left| \frac{1}{6} \, x - 1 \right|}}{{\left(x - 3\right)}^{2}}

graph = plot(R(x), (x, 0, 2.6), thickness=1.5, color='red',
             axes_labels=[r"$x$", r"$r/r_0$"]) \
        + plot(R(x), (x, 2.9, 12), thickness=1.5, color='red')
graph += line([(X1, 0), (X1, 4)], color='blue', linestyle='dashed') \
         + text(r'$x=x_1$', (2.7, 1.5), fontsize=16, color='blue', rotation='vertical')
graph += line([(X2, 0), (X2, 4)], color='green', linestyle='dashed') \
         + text(r'$x=x_2$', (5.7, 1.5), fontsize=16, color='green', rotation='vertical')
graph += line([(X1 + X2, 0), (X1 + X2, 4)], color='black', linestyle='dotted') \
         + text(r'$x=\alpha^{-1}$', (8.7, 1.5), fontsize=16, color='black', rotation='vertical')
show(graph, ymax=3, frame=True, axes=False, gridlines=True)
graph.save('vai_r_x_naksing.pdf', ymax=3, frame=True, axes=False, gridlines=True)
Image in a Jupyter notebook
f2(x) = 2 /(x*(a1*(x - x1)*(x - x2)))
f2(x)

2 (x1+x2)(x−x1)(x−x2)x\displaystyle \frac{2 \, {\left(x_{1} + x_{2}\right)}}{{\left(x - x_{1}\right)} {\left(x - x_{2}\right)} x}

s = integrate(f2(x), x).simplify_full()
s

2 (x22(log⁡(x−x1)−log⁡(x))−x12log⁡(x−x2)+x12log⁡(x)+(x1log⁡(x−x1)−x1log⁡(x−x2))x2)x12x2−x1x22\displaystyle \frac{2 \, {\left(x_{2}^{2} {\left(\log\left(x - x_{1}\right) - \log\left(x\right)\right)} - x_{1}^{2} \log\left(x - x_{2}\right) + x_{1}^{2} \log\left(x\right) + {\left(x_{1} \log\left(x - x_{1}\right) - x_{1} \log\left(x - x_{2}\right)\right)} x_{2}\right)}}{x_{1}^{2} x_{2} - x_{1} x_{2}^{2}}

s.coefficient(log(x - x1)).simplify_full()

2 (x1+x2)x12−x1x2\displaystyle \frac{2 \, {\left(x_{1} + x_{2}\right)}}{x_{1}^{2} - x_{1} x_{2}}

s.coefficient(log(x - x2)).simplify_full()

−2 (x1+x2)x1x2−x22\displaystyle -\frac{2 \, {\left(x_{1} + x_{2}\right)}}{x_{1} x_{2} - x_{2}^{2}}

s.coefficient(log(x)).simplify_full()

2 (x1+x2)x1x2\displaystyle \frac{2 \, {\left(x_{1} + x_{2}\right)}}{x_{1} x_{2}}

Marginal case: v0=16mv_0 = 16 m

f(x) = (8 - x)/(x - 4)^2
f(x)

−x−8(x−4)2\displaystyle -\frac{x - 8}{{\left(x - 4\right)}^{2}}

integrate(f(x), x)

−4x−4−log⁡(x−4)\displaystyle -\frac{4}{x - 4} - \log\left(x - 4\right)