# coding: utf-8 # # Timelike geodesic in Schwarzschild spacetime # # This Jupyter/SageMath worksheet presents the numerical computation of a timelike geodesic in Schwarzschild spacetime, given an initial point and tangent vector. It uses the integrated_geodesic functionality introduced by Karim Van Aelst in **SageMath 8.1**, in the framework of the [SageManifolds](http://sagemanifolds.obspm.fr) project. # # Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.1/SM_simple_goed_Schwarz.ipynb) to download the worksheet file (ipynb format). To run it, you must start SageMath within the Jupyter notebook, via the command sage -n jupyter # A version of SageMath at least equal to 8.1 is required to run this worksheet: # In[1]: version() # In[2]: get_ipython().run_line_magic('display', 'latex # LaTeX rendering turned on') # We define first the spacetime manifold $M$ and the standard Schwarzschild-Droste coordinates on it: # In[3]: M = Manifold(4, 'M') X. = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:\phi') X # For graphical purposes, we introduce $\mathbb{R}^3$ and some coordinate map $M\to \mathbb{R}^3$: # In[4]: R3 = Manifold(3, 'R^3', latex_name=r'\mathbb{R}^3') X3. = R3.chart() to_R3 = M.diff_map(R3, {(X, X3): [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)]}) to_R3.display() # Then, we define the Schwarzschild metric: # In[5]: g = M.lorentzian_metric('g') m = var('m'); assume(m >= 0) g[0,0], g[1,1] = -(1-2*m/r), 1/(1-2*m/r) g[2,2], g[3,3] = r^2, (r*sin(th))^2 g.display() # We pick an initial point and an initial tangent vector: # In[6]: p0 = M.point((0, 8*m, pi/2, 1e-12), name='p_0') v0 = M.tangent_space(p0)((1.297513, 0, 0, 0.0640625/m), name='v_0') v0.display() # We declare a geodesic with such initial conditions, denoting by $s$ the affine parameter (proper time), with $(s_{\rm min}, s_{\rm max})=(0, 1500\,m)$: # In[7]: s = var('s') geod = M.integrated_geodesic(g, (s, 0, 1500), v0); geod # We ask for the numerical integration of the geodesic, providing some numerical value for the parameter $m$, and then plot it in terms of the Cartesian chart X3 of $\mathbb{R}^3$: # In[8]: sol = geod.solve(parameters_values={m: 1}) # numerical integration interp = geod.interpolate() # interpolation of the solution for the plot graph = geod.plot_integrated(chart=X3, mapping=to_R3, plot_points=500, thickness=2, label_axes=False) # the geodesic graph += p0.plot(chart=X3, mapping=to_R3, size=4, parameters={m: 1}) # the starting point graph += sphere(size=2, color='grey') # the event horizon show(graph, viewer='threejs', online=True) # A Tachyon view of the geodesic: # In[9]: show(graph, viewer='tachyon', aspect_ratio=1, figsize=10) # Some details about the system solved to get the geodesic: # In[10]: geod.system(verbose=True) # We recognize in the above list the Christoffel symbols of the metric $g$: # In[11]: g.christoffel_symbols_display() # In[12]: g.riemann().display_comp()