version()
'SageMath version 9.0.beta4, Release Date: 2019-11-02'
%display latex
M = Manifold(4, 'M', structure='Lorentzian')
# X.<t, r, th, ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):periodic:\phi')
X.<t, r, th, ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:\phi')
X
g = M.metric()
m = 1
g[0,0] = -(1-2*m/r)
g[1,1] = 1/(1-2*m/r)
g[2,2] = r^2
g[3,3] = (r*sin(th))^2
g.display()
Xcart.<t, x, y, z> = M.chart()
X_to_Xcart = X.transition_map(Xcart, [t, r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])
X_to_Xcart.display()
M.top_charts()
M._top_charts = [X]
M.top_charts()
M.identity_map().coord_functions(X, Xcart)
We consider a null geodesic in the equatorial plane $\theta=\pi/2$ arising from far away of the black hole with an impact parameter very close to the critical value $b_c = 3\sqrt{3} m$:
bc = 3*sqrt(3)
n(bc)
b = 5.196155
n(b - bc)
We introduce a function to compute the initial tangent vector from the value of $b$, the conserved energy $E$ and the coordinates $(r_0,\varphi_0)$ of the initial point:
def initial_vector(r0, b, ph0=0, E=1, inward=False):
t0, th0 = 0, pi/2
L = b*E
vt0 = 1/(1-2*m/r0)
vr0 = sqrt(E^2 - L^2/r0^2*(1-2*m/r0))
if inward:
vr0 = - vr0
vth0 = 0
vph0 = L / r0^2
p0 = M((t0, r0, th0, ph0), name='p_0')
return M.tangent_space(p0)((vt0, vr0, vth0, vph0), name='v_0')
r0 = 15 # distance of the source
import time
graph = circle((0, 0), 2*m, edgecolor='black', fill=True, facecolor='grey', alpha=0.5)
ph0 = asin(b/r0)
v0 = initial_vector(r0, b, ph0=ph0, inward=True)
s = var('s')
geod0 = M.integrated_geodesic(g, (s, 0, 50), v0)
time0 = time.time()
#sol = geod0.solve(step=0.1, method="odeint")
sol = geod0.solve(step=0.1, method="odeint", rtol=1.e-12, atol=1.e-12)
print("CPU time in solve(): {} s".format(time.time() - time0))
interp = geod0.interpolate()
graph += geod0.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='green')
show(graph, xmin=-10, xmax=10, ymin=-10, ymax=10, axes_labels=[r'$x/m$', r'$y/m$'])
CPU time in solve(): 0.026593446731567383 s
geod1 = M.integrated_geodesic(g, (s, 0, 50), v0)
time0 = time.time()
sol = geod1.solve(step=0.1, method="rkf45")
print("CPU time in solve(): {} s".format(time.time() - time0))
interp = geod1.interpolate()
graph += geod1.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='red')
show(graph, xmin=-10, xmax=10, ymin=-10, ymax=10, axes_labels=[r'$x/m$', r'$y/m$'])
CPU time in solve(): 4.370987176895142 s
# graph.save('test_geod.pdf', xmin=-10, xmax=10, ymin=-10, ymax=10, axes_labels=[r'$x/m$', r'$y/m$'])
graphics_array([plot(lambda s: X(geod0(s))[0] - X(geod1(s))[0], (0, 50),
axes_labels=['$s$', r'$\Delta t$']),
plot(lambda s: X(geod0(s))[1] - X(geod1(s))[1], (0, 50),
axes_labels=['$s$', r'$\Delta r$']),
plot(lambda s: X(geod0(s))[3] - X(geod1(s))[3], (0, 50),
axes_labels=['$s$', r'$\Delta \varphi$'])])
Comparison with step=0.02
instead of 0.1
:
geod = M.integrated_geodesic(g, (s, 0, 50), v0)
time0 = time.time()
sol = geod.solve(step=0.02, method="rkf45")
print("CPU time in solve(): {} s".format(time.time() - time0))
interp = geod.interpolate()
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='red')
show(graph, xmin=-10, xmax=10, ymin=-10, ymax=10, axes_labels=[r'$x/m$', r'$y/m$'])
CPU time in solve(): 21.649514198303223 s
Comparison of the two rkf45 methods, one with step=0.02
and the otjer with step=0.1
:
graphics_array([plot(lambda s: X(geod(s))[0] - X(geod1(s))[0], (0, 50),
axes_labels=['$s$', r'$\Delta t$']),
plot(lambda s: X(geod(s))[1] - X(geod1(s))[1], (0, 50),
axes_labels=['$s$', r'$\Delta r$']),
plot(lambda s: X(geod(s))[3] - X(geod1(s))[3], (0, 50),
axes_labels=['$s$', r'$\Delta \varphi$'])])
geod = M.integrated_geodesic(g, (s, 0, 50), v0)
time0 = time.time()
sol = geod.solve(step=0.1, method="rk4_maxima")
print("CPU time in solve(): {} s".format(time.time() - time0))
interp = geod.interpolate()
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='purple')
show(graph, xmin=-10, xmax=10, ymin=-10, ymax=10, axes_labels=[r'$x/m$', r'$y/m$'])
CPU time in solve(): 13.302617311477661 s
graphics_array([plot(lambda s: X(geod0(s))[0] - X(geod(s))[0], (0, 50),
axes_labels=['$s$', r'$\Delta t$']),
plot(lambda s: X(geod0(s))[1] - X(geod(s))[1], (0, 50),
axes_labels=['$s$', r'$\Delta r$']),
plot(lambda s: X(geod0(s))[3] - X(geod(s))[3], (0, 50),
axes_labels=['$s$', r'$\Delta \varphi$'])])
geod = M.integrated_geodesic(g, (s, 0, 50), v0)
time0 = time.time()
sol = geod.solve(step=0.1, method="gear2")
print("CPU time in solve(): {} s".format(time.time() - time0))
interp = geod.interpolate()
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='gold')
show(graph, xmin=-10, xmax=10, ymin=-10, ymax=10, axes_labels=[r'$x/m$', r'$y/m$'])
CPU time in solve(): 22.341139554977417 s
graphics_array([plot(lambda s: X(geod0(s))[0] - X(geod(s))[0], (0, 50),
axes_labels=['$s$', r'$\Delta t$']),
plot(lambda s: X(geod0(s))[1] - X(geod(s))[1], (0, 50),
axes_labels=['$s$', r'$\Delta r$']),
plot(lambda s: X(geod0(s))[3] - X(geod(s))[3], (0, 50),
axes_labels=['$s$', r'$\Delta \varphi$'])])
geod = M.integrated_geodesic(g, (s, 0, 50), v0)
time0 = time.time()
sol = geod.solve(step=0.1, method="rk4")
print("CPU time in solve(): {} s".format(time.time() - time0))
interp = geod.interpolate()
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='cyan')
show(graph, xmin=-10, xmax=10, ymin=-10, ymax=10, axes_labels=[r'$x/m$', r'$y/m$'])
CPU time in solve(): 7.399372816085815 s
graphics_array([plot(lambda s: X(geod0(s))[0] - X(geod(s))[0], (0, 50),
axes_labels=['$s$', r'$\Delta t$']),
plot(lambda s: X(geod0(s))[1] - X(geod(s))[1], (0, 50),
axes_labels=['$s$', r'$\Delta r$']),
plot(lambda s: X(geod0(s))[3] - X(geod(s))[3], (0, 50),
axes_labels=['$s$', r'$\Delta \varphi$'])])
geod = M.integrated_geodesic(g, (s, 0, 50), v0)
time0 = time.time()
sol = geod.solve(step=0.1, method="dopri5")
print("CPU time in solve(): {} s".format(time.time() - time0))
interp = geod.interpolate()
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='orange')
show(graph, xmin=-10, xmax=10, ymin=-10, ymax=10, axes_labels=[r'$x/m$', r'$y/m$'])
CPU time in solve(): 0.07653164863586426 s
graphics_array([plot(lambda s: X(geod0(s))[0] - X(geod(s))[0], (0, 50),
axes_labels=['$s$', r'$\Delta t$']),
plot(lambda s: X(geod0(s))[1] - X(geod(s))[1], (0, 50),
axes_labels=['$s$', r'$\Delta r$']),
plot(lambda s: X(geod0(s))[3] - X(geod(s))[3], (0, 50),
axes_labels=['$s$', r'$\Delta \varphi$'])])