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

BHLectures/sage / Vaidya_solve_ode_out.ipynb

²¹⁰⁶⁰ views

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.

In [1]:

version()

Out[1]:

'SageMath version 10.0.rc0, Release Date: 2023-04-23'

In [2]:

%display latex

In [3]:

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

In [4]:

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

In [5]:

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

Out[5]:

$\displaystyle -\frac{{\alpha} x - 1}{{\alpha} x^{2} - x + 2}$

Case of a high density radiation shell

In [6]:

assume(a > 1/8)

In [7]:

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

Out[7]:

$\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)$

In [8]:

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

In [9]:

x1 = 1/(2*a)

In [10]:

P(x1)

Out[10]:

$\displaystyle -\frac{1}{4 \, {\alpha}} + 2$

In [11]:

F(x1)

Out[11]:

$\displaystyle -\frac{1}{2} \, \log\left(-\frac{1}{4 \, {\alpha}} + 2\right)$

In [12]:

exp(F(x1)).simplify_full()

Out[12]:

$\displaystyle \frac{2}{\sqrt{\frac{8 \, {\alpha} - 1}{{\alpha}}}}$

In [13]:

P(1/a)

Out[13]:

$\displaystyle 2$

In [14]:

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

Out[14]:

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

In [15]:

plot(thb, (x, 0.001, 15.99))

Out[15]:

In [16]:

taylor(thb, x, 0, 3)

Out[16]:

$\displaystyle -\frac{1}{3840} \, x^{3} - \frac{1}{96} \, x^{2} + \frac{1}{2} \, x - 4$

In [17]:

Out[17]:

$\displaystyle x$

In [18]:

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

Out[18]:

$\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}}$

In [19]:

diff(f(a, x), x)

Out[19]:

$\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)}}$

In [20]:

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)

Out[20]:

Equivalent form:

In [21]:

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

Out[21]:

$\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}}$

In [22]:

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)

Out[22]:

Maximal value of $f_\alpha(x)$ :

In [23]:

f(a, 1/a)

Out[23]:

$\displaystyle e^{\left(\frac{2 \, \arctan\left(\frac{1}{\sqrt{8 \, {\alpha} - 1}}\right)}{\sqrt{8 \, {\alpha} - 1}}\right)}$

Expansion for small $x$ :

In [24]:

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

Out[24]:

$\displaystyle -\frac{1}{4} \, {\left({\alpha} - 1\right)} x^{2} + \frac{1}{2} \, x + 1$

Case of a low density radiation shell

In [25]:

forget(a > 1/8)

In [26]:

assume(a < 1/8)

In [27]:

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

Out[27]:

$\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)$

In [28]:

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

Out[28]:

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

In [29]:

(x01 + x02).simplify_full()

Out[29]:

$\displaystyle \frac{1}{{\alpha}}$

In [30]:

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

Out[30]:

$\displaystyle \mathrm{True}$

In [31]:

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

Out[31]:

$\displaystyle \frac{1}{x_{1} + x_{2}}$

In [32]:

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

Out[32]:

$\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)}}$

In [33]:

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

Out[33]:

$\displaystyle \frac{x_{2} \log\left(x - x_{1}\right) - x_{1} \log\left(x - x_{2}\right)}{x_{1} - x_{2}}$

In [34]:

(x01 - x02).simplify_full()

Out[34]:

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

In [35]:

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

Out[35]:

$\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}}$

In [36]:

diff(g, x)

Out[36]:

$\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)}}$

In [37]:

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

Out[37]:

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

In [38]:

n(X1), n(X2)

Out[38]:

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

In [39]:

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

Out[39]:

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

In [40]:

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

Out[40]:

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

In [41]:

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)

Out[41]:

In [42]:

f2(x) = 2 /(x*(a1*(x - x1)*(x - x2)))
f2(x)

Out[42]:

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

In [43]:

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

Out[43]:

$\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}}$

In [44]:

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

Out[44]:

$\displaystyle \frac{2 \, {\left(x_{1} + x_{2}\right)}}{x_{1}^{2} - x_{1} x_{2}}$

In [45]:

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

Out[45]:

$\displaystyle -\frac{2 \, {\left(x_{1} + x_{2}\right)}}{x_{1} x_{2} - x_{2}^{2}}$

In [46]:

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

Out[46]:

$\displaystyle \frac{2 \, {\left(x_{1} + x_{2}\right)}}{x_{1} x_{2}}$

Marginal case: $v_0 = 16 m$

In [47]:

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

Out[47]:

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

In [48]:

integrate(f(x), x)

Out[48]:

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

Solving the ODE for outgoing radial null geodesics in Vaidya spacetime

Case of a high density radiation shell

Case of a low density radiation shell

Marginal case: $v_0 = 16 m$

Product

Resources

Company

Solving the ODE for outgoing radial null geodesics in Vaidya spacetime

Case of a high density radiation shell

Case of a low density radiation shell

Marginal case: v0=16mv_0 = 16 mv0​=16m

Marginal case: $v_0 = 16 m$