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
.
%display latex
P(u, b) = 2*u^3 - u^2 + 1/b^2
P(u, b)
bc = 3*sqrt(3)
assume(b < bc)
xi = bc/b - sqrt((bc/b)^2 - 1)
u_neg(b) = (1 - xi^(2/3) - xi^(-2/3))/6
u_neg(b)
Check:
P(u_neg(b), b).simplify_full().factor()
plot(u_neg, (0.1, bc), axes_labels=[r'$b/m$', r'$u_{\rm neg}$'],
frame=True, gridlines=True)
u_neg(b).taylor(b, 0, 2)
lim(u_neg(b), b=0)
u0(b) = 1/4 - u_neg(b)/2
us(b) = sqrt(u_neg(b)*(3*u_neg(b) - 1)) + u_neg(b)
lim(u0(b), b=bc)
lim(us(b), b=bc)
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)
g.save("gis_u0_us_b.pdf", ymin=0, ymax=2.5)
plot(lambda b: us(b) - u0(b), (0.1, bc), axes_labels=[r'$b/m$', r'$u_* - u_0$'],
frame=True, gridlines=True)
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()
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()
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
g.save("gis_k_b_lt_bc.pdf")
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)
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)
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)
def phif(u, b):
us = usf(b)
un = unf(b)
return arccos(abs(us - u)/(us + u - 2*un))
plot(lambda u: phif(u, 4), (0, 6), axes_labels=[r'$u$', r'$\phi(u)$'],
frame=True, gridlines=True)
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))
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))
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
un = var('u_n')
a12 = (6*un + 1)*(2*un - 1)/16
a12
b1 = 1/4 - un/2
b1
A = sqrt((b1 - un)^2 + a12).simplify_full()
A
Ab(b) = A.subs({un: u_neg(b)}).simplify_full().factor()
Ab(b)
plot(Ab, (0.1, bc), axes_labels=[r'$b/m$', r'$A$'],
frame=True, gridlines=True)
Ab(bc)
k = sqrt((A + b1 - un)/(2*A))
k
kb(b) = k.subs({un: u_neg(b)}).simplify_full().factor()
kb
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
kb(b).taylor(b, 0, 2)
k0 = lim(kb(b), b=0)
k0
n(k0)
k0 = k0.canonicalize_radical()
k0
n(k0)
t = (A + un - u)/(A - un + u)
t
tb(u, b) = t.subs({un: u_neg(b)}).simplify_full().factor()
tb(u, b)
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
show(g, xmin=0, xmax=0.3, ymin=0.2, ymax=0.4)
unx = - 1/6 - x
bx = bc - 81*sqrt(3)/4*x
bx
bool(u_neg(bx).simplify_full().taylor(x, 0, 1) == unx)
A.subs({un: unx}).taylor(x, 0, 1)
k.subs({un: unx}).taylor(x, 0, 1)
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
forget()
assume(u>0)
assume(u<pi/2)
integrate(1/cos(x), x, 0, u, algorithm='giac')
forget()
assume(u>pi/2)
assume(u<pi)
integrate(1/cos(x), x, pi/2, u, algorithm='giac')