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
.
version()
'SageMath version 9.1.beta3, Release Date: 2020-02-02'
%display latex
r, b, m = var('r b m', domain='real')
P(r, b, m) = r^3 - b^2*r + 2*m*b^2
P(r, b, m)
r0 = 2/sqrt(3)*b*cos(pi/3 - arccos(3*sqrt(3)*m/b)/3)
r0
P(r0, b, m).simplify_full().trig_reduce()
r1 = 2/sqrt(3)*b*cos(pi - arccos(3*sqrt(3)*m/b)/3)
r1
P(r1, b, m).simplify_full().trig_reduce()
r2 = 2/sqrt(3)*b*cos(5*pi/3 - arccos(3*sqrt(3)*m/b)/3)
r2
P(r2, b, m).simplify_full().trig_reduce()
rm(b) = r1.subs({m: 1})
rp(b) = r0.subs({m: 1})
ra(b) = r2.subs({m: 1})
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
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
show(g, xmin=5, ymin=2)
g.save('ges_null_per_apo.pdf', xmin=5, ymin=2)
bcrit = 3*sqrt(3)
rpx = rp(bcrit + x).taylor(x, 0, 2)
rpx
rax = ra(bcrit + x).taylor(x, 0, 2)
rax
upx = (1/rpx).taylor(x, 0, 2)
upx
uax = (1/rax).taylor(x, 0, 2)
uax
up = (1 - x)/3
bx = 1/(up*sqrt(1 - 2*up)).simplify_full()
bx
upx = (1/rp(bx)).simplify_full().taylor(x, 0, 3)
upx
P(1/upx, bx, 1).simplify_full()
uax = (1/ra(bx)).simplify_full().taylor(x, 0, 3)
uax
P(1/uax, bx, 1).simplify_full().taylor(x, 0, 3)
umx = (1/rm(bx)).simplify_full().taylor(x, 0, 3)
umx
P(1/umx, bx, 1).simplify_full().taylor(x, 0, 3)
Another check:
umx + upx + uax
P(u, b) = 2*u^3 - u^2 + 1/b^2
P
b_crit = 3*sqrt(3)
b_sel = [3, 4, b_crit, 7, 20]
b_sel
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)
graph.save('ges_polynomial_b_u.pdf', ymin=-0.2, ymax=0.2)
The depressed polynomial:
b = var('b')
pd = (P(x + 1/6, b)/2).expand()
pd
p = -1/12
q = 1/(2*b^2) - 1/108
p, q
Roots via Viète formula:
3*q/(2*p)*sqrt(-3/p)
psi = arcsin(b_crit/b)
un(b) = 1/3*cos(2*psi/3 + 2*pi/3) + 1/6
un(b)
P(un(b), b).simplify_full().trig_reduce().trig_expand()
up(b) = 1/3*cos(2*psi/3 + 4*pi/3) + 1/6
up(b)
P(up(b), b).simplify_full().trig_reduce().trig_expand()
ua(b) = 1/3*cos(2*psi/3) + 1/6
ua(b)
P(ua(b), b).simplify_full().trig_reduce().trig_expand()
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
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:
P(u_neg(b), b).simplify_full().factor()
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')
show(g, xmax=15, ymin=-0.5, ymax=0.5)
g.save('gis_u_per_apo_neg.pdf', xmax=15, ymin=-0.5, ymax=0.5)
k(b) = sqrt((up(b) - un(b))/(ua(b) - un(b))).simplify_full()
k(b)
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
g.save("gis_elliptic_mod.pdf")
k1(b) = sqrt(2) / sqrt(sqrt(3)*cot(2/3*arcsin(3*sqrt(3)/b)) + 1)
k1(b)
(k(b)^2 - k1(b)^2).simplify_full()
tb(u, b) = sqrt((up(b) - u)/(ua(b) - u))/k(b)
tb(u, b)
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
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
upx = 1/3 - x
bx = 1/(upx*sqrt(1 - 2*upx)).simplify_full()
bx
up(bx).taylor(x, 0, 2)
bx.taylor(x, 0, 2)
uax = ua(bx).taylor(x, 0, 2)
uax
unx = un(bx).taylor(x, 0, 2)
unx
kx = k(bx).taylor(x, 0, 2)
kx
sin_ax = (sqrt(upx/uax)/kx).taylor(x, 0, 2)
sin_ax
ax = (arcsin(sqrt(upx/uax)/kx)).taylor(x, 0, 2)
ax
tan_ax = tan(arcsin(sqrt(upx/uax)/kx)).taylor(x, 0, 2)
tan_ax
(tan_ax * sqrt(1 - k(bx)^2)).taylor(x, 0, 2)
(1/_).taylor(x, 0, 2)
sin(arctan(1/sqrt(2))).simplify_full()
arctan(1/sqrt(2)).simplify_full()
sqrt(1 - k(bx)^2).taylor(x, 0, 2)
aax = arcsin(sqrt(upx/uax)/kx)
(tan(aax) + sec(aax)).taylor(x, 0, 2)
(3*sqrt(3)/bx).taylor(x, 0, 2)
(2/3)*arcsin((3*sqrt(3)/bx)).taylor(x, 0, 2)