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

Path: BHLectures/sage / gis_paramaters_b_lt_bc.ipynb

Views: ²⁰¹¹⁰

Kernel: SageMath 9.1.beta1

Null geodesics with $b < b_{\rm crit}$

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

To run it, you must start SageMath with sage -n jupyter.

In [1]:

%display latex

In [2]:

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

In [3]:

bc = 3*sqrt(3)
assume(b < bc)

The unique real zero of $P_b$

In [4]:

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

Check:

In [5]:

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

In [6]:

plot(u_neg, (0.1, bc), axes_labels=[r'$b/m$', r'$u_{\rm neg}$'],
     frame=True, gridlines=True)

In [7]:

u_neg(b).taylor(b, 0, 2)

In [8]:

lim(u_neg(b), b=0)

In [9]:

u0(b) = 1/4 - u_neg(b)/2
us(b) = sqrt(u_neg(b)*(3*u_neg(b) - 1)) + u_neg(b)

In [10]:

lim(u0(b), b=bc)

In [11]:

lim(us(b), b=bc)

In [12]:

g = plot(us, (0.05, bc), legend_label=r'$u_*$',
         thickness=1.5, color='red',
         axes_labels=[r'$b/m$', r'$u_0,\ u_*$'],
         frame=True, gridlines=True) \
    + plot(u0, (0.05, bc), legend_label=r'$u_0$',
           thickness=1.5, color='blue') \
    + line([(0, 1/3), (bc, 1/3)], color='black', linestyle='--')
show(g, ymin=0, ymax=2.5)

In [13]:

g.save("gis_u0_us_b.pdf", ymin=0, ymax=2.5)

In [14]:

plot(lambda b: us(b) - u0(b), (0.1, bc), axes_labels=[r'$b/m$', r'$u_* - u_0$'],
     frame=True, gridlines=True)

Factorization of $P_b$ :

In [15]:

P1 = 2*(u - u_neg(b))*((u - 1/4 + u_neg(b)/2)^2 + (6*u_neg(b) + 1)*(2*u_neg(b) - 1)/16)
(P1 - P(u, b)).simplify_full().factor()

In [16]:

P2 = 2*(u - u_neg(b))*((u - u0(b))^2 + (us(b) - u_neg(b))^2 - (u0(b) - u_neg(b))^2)
(P2 - P(u, b)).simplify_full().factor()

Elliptic integral

In [17]:

kb(b) = sqrt((us(b) - 5/2*u_neg(b) + 1/4)/(2*(us(b) - u_neg(b))))
g = plot(kb, (0.001, bc), axes_labels=[r'$b/m$', r'$k$'],
         thickness=1.5, color='red',
         frame=True, gridlines=True, axes=False)
g

In [18]:

g.save("gis_k_b_lt_bc.pdf")

In [19]:

def unf(b):
    xi = bc/b - sqrt((bc/b)^2 - 1)
    return (1 - xi^(2/3) - xi^(-2/3))/6

plot(unf, (0.1, bc), axes_labels=[r'$b/m$', r'$u_{\rm n}$'],
     frame=True, gridlines=True,ymax=0)

In [20]:

def usf(b):
    un = unf(b)
    return sqrt(un*(3*un - 1)) + un

plot(usf, (0.1, bc), axes_labels=[r'$b/m$', r'$u_*$'],
     frame=True, gridlines=True)

In [21]:

def kf(b):
    un = unf(b)
    us = usf(b)
    return sqrt((us - 5/2*un + 1/4)/(2*(us - un)))

plot(kf, (0.1, bc), axes_labels=[r'$b/m$', r'$k$'],
     frame=True, gridlines=True)

In [22]:

def phif(u, b):
    us = usf(b)
    un = unf(b)
    return arccos(abs(us - u)/(us + u - 2*un))

In [23]:

plot(lambda u: phif(u, 4), (0, 6), axes_labels=[r'$u$', r'$\phi(u)$'],
     frame=True, gridlines=True)

In [24]:

def Phi1(u, b):
    k2 = kf(b)^2
    us = usf(b)
    un = unf(b)
    aa = elliptic_kc(k2) - elliptic_f(phif(u, b), k2)
    if u > us:
        aa = - aa
    return aa / sqrt(2*(us - un))

In [25]:

def Phi(u, b):
    u = RDF(u)
    b = RDF(b)
    bc = RDF(3*sqrt(3))
    xi = bc/b - sqrt((bc/b)^2 - 1)
    un = (1 - xi^(2/3) - xi^(-2/3))/6
    us = sqrt(un*(3*un - 1)) + un
    k2 = (us - 2.5*un + 0.25)/(2*(us - un))
    phi = arccos(abs(us - u)/(us + u - 2*un))
    aa = elliptic_kc(k2) - elliptic_f(phi, k2)
    if u > us:
        aa = - aa
    return aa / sqrt(2*(us - un))

In [26]:

g = Graphics()
for b1 in [5.19, 5.15, 5, 3, 1, 1/2]:
    g += plot(lambda u: Phi(u, b1), (0.01, 1.5),
              legend_label="b={}".format(float(b1)), 
              color=hue(b1/bc), axes_labels=[r'$u$', r'$\Phi_b(u)$'], 
              frame=True, gridlines=True, axes=False)
g

In [27]:

un = var('u_n')
a12 = (6*un + 1)*(2*un - 1)/16
a12

In [28]:

b1 = 1/4 - un/2
b1

In [29]:

A = sqrt((b1 - un)^2 + a12).simplify_full()
A

In [30]:

Ab(b) = A.subs({un: u_neg(b)}).simplify_full().factor()
Ab(b)

In [31]:

plot(Ab, (0.1, bc), axes_labels=[r'$b/m$', r'$A$'],
     frame=True, gridlines=True)

In [32]:

Ab(bc)

In [33]:

k = sqrt((A + b1 - un)/(2*A))
k

In [34]:

kb(b) = k.subs({un: u_neg(b)}).simplify_full().factor()
kb

In [35]:

g = plot(kb, (0.001, bc), axes_labels=[r'$b/m$', r'$k$'],
         thickness=1.5, color='red',
         frame=True, gridlines=True, axes=False)
g

In [36]:

kb(b).taylor(b, 0, 2)

In [37]:

k0 = lim(kb(b), b=0)
k0

In [38]:

n(k0)

In [39]:

k0 = k0.canonicalize_radical()
k0

In [40]:

n(k0)

In [41]:

t = (A + un - u)/(A - un + u)
t

In [42]:

tb(u, b) = t.subs({un: u_neg(b)}).simplify_full().factor()
tb(u, b)

In [43]:

g = Graphics()
for b1 in [bc, 5, 4, 3, 2, 1, 1/2]:
    g += plot(tb(u, b1), (u, u_neg(b1), 8), 
              legend_label="b={}".format(float(b1)), 
              color=hue(b1/bc), axes_labels=[r'$u$', r'$t$'], 
              frame=True, gridlines=True)
g

In [44]:

show(g, xmin=0, xmax=0.3, ymin=0.2, ymax=0.4)

In [45]:

unx = - 1/6 - x

In [46]:

bx = bc - 81*sqrt(3)/4*x
bx

In [47]:

bool(u_neg(bx).simplify_full().taylor(x, 0, 1) == unx)

In [48]:

A.subs({un: unx}).taylor(x, 0, 1)

In [49]:

k.subs({un: unx}).taylor(x, 0, 1)

In [50]:

g = Graphics()
for b1 in [bc, 5, 4, 3, 2, 1, 1/2]:
    g += plot(arccos(tb(u, b1)), (u, u_neg(b1), 8), 
              legend_label="b={}".format(float(b1)), 
              color=hue(b1/bc), axes_labels=[r'$u$', r'$\phi_u$'], 
              frame=True, gridlines=True)
g

In [51]:

forget()
assume(u>0)
assume(u<pi/2)
integrate(1/cos(x), x, 0, u, algorithm='giac')

In [52]:

forget()
assume(u>pi/2)
assume(u<pi)
integrate(1/cos(x), x, pi/2, u, algorithm='giac')

In [ ]:

Null geodesics with b<bcritb < b_{\rm crit}b<bcrit​

The unique real zero of PbP_bPb​

Factorization of PbP_bPb​:

Elliptic integral

Null geodesics with $b < b_{\rm crit}$

The unique real zero of $P_b$

Factorization of $P_b$ :