CoCalc -- ges_null_periastron.ipynb

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

Path: BHLectures/sage / ges_null_periastron.ipynb

Views: ²⁰¹¹⁰

Kernel: SageMath 9.1.beta3

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

Click here to download the notebook file (ipynb format). To run it, you must start SageMath with sage -n jupyter.

In [2]:

version()

'SageMath version 9.1.beta3, Release Date: 2020-02-02'

In [3]:

%display latex

The cubic polynomial

In [4]:

r, b, m = var('r b m', domain='real')

In [5]:

P(r, b, m) = r^3 - b^2*r + 2*m*b^2
P(r, b, m)

Roots via Viète's trigonometric method

In [6]:

r0 = 2/sqrt(3)*b*cos(pi/3 - arccos(3*sqrt(3)*m/b)/3)
r0

In [7]:

P(r0, b, m).simplify_full().trig_reduce()

In [8]:

r1 = 2/sqrt(3)*b*cos(pi - arccos(3*sqrt(3)*m/b)/3)
r1

In [9]:

P(r1, b, m).simplify_full().trig_reduce()

In [10]:

r2 = 2/sqrt(3)*b*cos(5*pi/3 - arccos(3*sqrt(3)*m/b)/3)
r2

In [11]:

P(r2, b, m).simplify_full().trig_reduce()

In [12]:

rm(b) = r1.subs({m: 1})
rp(b) = r0.subs({m: 1})
ra(b) = r2.subs({m: 1})

In [13]:

g = plot(rp(b), (b, 3*sqrt(3), 8), color='red', legend_label=r'$r_{\rm per}$',
         axes_labels=[r'$b/m$', r'$r$']) \
    + plot(ra(b), (b, 3*sqrt(3), 8), color='green', legend_label=r'$r_{\rm apo}$') \
    + plot(rm(b), (b, 3*sqrt(3), 8), color='blue', legend_label=r'$r_{\rm neg}$')
g

In [14]:

g = plot(rp(b), (b, 3*sqrt(3), 10), color='red', thickness=1.5,
         legend_label=r'$r_{\rm per}$', axes_labels=[r'$b/m$', r'$r/m$'], 
         frame=True, gridlines=True, aspect_ratio=1) \
    + plot(ra(b), (b, 3*sqrt(3), 10), color='green', thickness=1.5,
           legend_label=r'$r_{\rm apo}$')
g

In [15]:

show(g, xmin=5, ymin=2)

In [16]:

g.save('ges_null_per_apo.pdf', xmin=5, ymin=2)

Expansions close to $b_{\rm crit}$

Expansion of the form $b = b_{\rm crit} + x$

In [17]:

bcrit = 3*sqrt(3)

In [18]:

rpx = rp(bcrit + x).taylor(x, 0, 2)
rpx

In [19]:

rax = ra(bcrit + x).taylor(x, 0, 2)
rax

In [20]:

upx = (1/rpx).taylor(x, 0, 2)
upx

In [21]:

uax = (1/rax).taylor(x, 0, 2)
uax

Expansion of the form $u_{\rm per} = \frac{1}{3} (1 - x)$

In [22]:

up = (1 - x)/3

In [23]:

bx = 1/(up*sqrt(1 - 2*up)).simplify_full()
bx

In [24]:

upx = (1/rp(bx)).simplify_full().taylor(x, 0, 3)
upx

In [25]:

P(1/upx, bx, 1).simplify_full()

In [26]:

uax = (1/ra(bx)).simplify_full().taylor(x, 0, 3)
uax

In [27]:

P(1/uax, bx, 1).simplify_full().taylor(x, 0, 3)

In [28]:

umx = (1/rm(bx)).simplify_full().taylor(x, 0, 3)
umx

In [29]:

P(1/umx, bx, 1).simplify_full().taylor(x, 0, 3)

Another check:

In [30]:

umx + upx + uax

The cubic polynomial $P_b(u)$

In [31]:

P(u, b) = 2*u^3 - u^2 + 1/b^2
P

In [32]:

b_crit = 3*sqrt(3)
b_sel = [3, 4, b_crit, 7, 20]
b_sel

In [33]:

graph = Graphics()
for b in b_sel:
    if b == b_crit:
        legend_label = r'$b = 3\sqrt{3}\, m$'
    else:
        legend_label=r'$b = {:.0f} \, m$'.format(float(b))
    graph += plot(P(u, b), (u, -0.5, 0.8), color = hue((b - b_crit)/6), 
                  legend_label=legend_label, thickness=1.5,
                  axes_labels=[r"$u$", r"$P_b(u)$"])
show(graph, ymin=-0.2, ymax=0.2)

In [34]:

graph.save('ges_polynomial_b_u.pdf', ymin=-0.2, ymax=0.2)

Roots of $P_b(u)$

The depressed polynomial:

In [35]:

b = var('b')
pd = (P(x + 1/6, b)/2).expand()
pd

In [36]:

p = -1/12
q = 1/(2*b^2) - 1/108
p, q

Roots via Viète formula:

In [37]:

3*q/(2*p)*sqrt(-3/p)

In [38]:

psi = arcsin(b_crit/b)
un(b) = 1/3*cos(2*psi/3 + 2*pi/3) + 1/6
un(b)

In [39]:

P(un(b), b).simplify_full().trig_reduce().trig_expand()

In [40]:

up(b) = 1/3*cos(2*psi/3 + 4*pi/3) + 1/6
up(b)

In [41]:

P(up(b), b).simplify_full().trig_reduce().trig_expand()

In [42]:

ua(b) = 1/3*cos(2*psi/3) + 1/6
ua(b)

In [43]:

P(ua(b), b).simplify_full().trig_reduce().trig_expand()

In [44]:

g = plot(ua, (b_crit, 16), color='green', thickness=1.5,
         legend_label=r'$u_{\rm a}$', axes_labels=[r'$b/m$', r'$u$'], 
         frame=True, gridlines=True) \
    + plot(up, (b_crit, 16), color='red', thickness=1.5,
           legend_label=r'$u_{\rm p}$') \
    + plot(un, (b_crit, 16), color='maroon', thickness=1.5,
           legend_label=r'$u_{\rm n}$')
g

Adding $u_{\rm n}$ for $b< b_{\rm c}$ :

In [45]:

xi = b_crit/b - sqrt((b_crit/b)^2 - 1)
u_neg(b) = (1 - xi^(2/3) - xi^(-2/3))/6
u_neg(b)

Check:

In [46]:

P(u_neg(b), b).simplify_full().factor()

In [47]:

g += plot(u_neg, (0.1, b_crit), color='maroon', thickness=1.5) \
     + line([(b_crit, -3.5), (b_crit, 0.6)], linestyle='dashed',
            color='black')

In [48]:

show(g, xmax=15, ymin=-0.5, ymax=0.5)

In [49]:

g.save('gis_u_per_apo_neg.pdf', xmax=15, ymin=-0.5, ymax=0.5)

Elliptic modulus $k$

In [50]:

k(b) = sqrt((up(b) - un(b))/(ua(b) - un(b))).simplify_full()
k(b)

In [51]:

g = plot(k(b), (b, 3*sqrt(3), 20), color='red', thickness=1.5,
         axes_labels=[r'$b/m$', r'$k$'], 
         frame=True, gridlines=True)
g

In [52]:

g.save("gis_elliptic_mod.pdf")

In [53]:

k1(b) = sqrt(2) / sqrt(sqrt(3)*cot(2/3*arcsin(3*sqrt(3)/b)) + 1)
k1(b)

In [54]:

(k(b)^2 - k1(b)^2).simplify_full()

Variable $t$

In [59]:

tb(u, b) = sqrt((up(b) - u)/(ua(b) - u))/k(b)
tb(u, b)

In [60]:

bc = b_crit
g = Graphics()
for b1 in [bc, 5.3, 5.5, 6, 8, 10, 100]:
    g += plot(tb(u, b1), (u, 0, up(b1)), 
              legend_label="b={}".format(float(b1)), 
              color=hue(b1/bc), axes_labels=[r'$u$', r'$t$'], 
              frame=True, gridlines=True)
g

Function $\phi(u, b)$

In [61]:

g = Graphics()
phib(u, b) = arcsin(tb(u, b))
for b1 in [bc, 5.3, 5.5, 6, 8, 10, 100]:
    g += plot(phib(u, b1), (u, 0, up(b1)), 
              legend_label="b={}".format(float(b1)), 
              color=hue(b1/bc), axes_labels=[r'$u$', r'$\phi(u, b)$'], 
              frame=True, gridlines=True)
g

Series expansion for $b$ close to $b_{\rm crit}$

In [54]:

upx = 1/3 - x

In [55]:

bx = 1/(upx*sqrt(1 - 2*upx)).simplify_full()
bx

In [56]:

up(bx).taylor(x, 0, 2)

In [57]:

bx.taylor(x, 0, 2)

In [58]:

uax = ua(bx).taylor(x, 0, 2)
uax

In [59]:

unx = un(bx).taylor(x, 0, 2)
unx

In [60]:

kx = k(bx).taylor(x, 0, 2)
kx

In [61]:

sin_ax = (sqrt(upx/uax)/kx).taylor(x, 0, 2)
sin_ax

In [62]:

ax = (arcsin(sqrt(upx/uax)/kx)).taylor(x, 0, 2)
ax

In [63]:

tan_ax = tan(arcsin(sqrt(upx/uax)/kx)).taylor(x, 0, 2)
tan_ax

In [64]:

(tan_ax * sqrt(1 - k(bx)^2)).taylor(x, 0, 2)

In [65]:

(1/_).taylor(x, 0, 2)

In [66]:

sin(arctan(1/sqrt(2))).simplify_full()

In [67]:

arctan(1/sqrt(2)).simplify_full()

In [68]:

sqrt(1 - k(bx)^2).taylor(x, 0, 2)

In [69]:

aax = arcsin(sqrt(upx/uax)/kx)
(tan(aax) + sec(aax)).taylor(x, 0, 2)

In [70]:

(3*sqrt(3)/bx).taylor(x, 0, 2)

In [71]:

(2/3)*arcsin((3*sqrt(3)/bx)).taylor(x, 0, 2)

In [ ]:

Periastron/apoastron of null geodesics in Schwarzschild spacetime and related quantities

The cubic polynomial

Roots via Viète's trigonometric method

Expansions close to bcritb_{\rm crit}bcrit​

Expansion of the form b=bcrit+xb = b_{\rm crit} + xb=bcrit​+x

Expansion of the form uper=13(1−x)u_{\rm per} = \frac{1}{3} (1 - x)uper​=31​(1−x)

The cubic polynomial Pb(u)P_b(u)Pb​(u)

Roots of Pb(u)P_b(u)Pb​(u)

Adding unu_{\rm n}un​ for b<bcb< b_{\rm c}b<bc​:

Elliptic modulus kkk

Variable ttt

Function ϕ(u,b)\phi(u, b)ϕ(u,b)

Series expansion for bbb close to bcritb_{\rm crit}bcrit​