{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Fast Kerr geodesics" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%display latex" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the manifold catalog (https://trac.sagemath.org/ticket/25869):" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "K. = manifolds.Kerr(2, -1, coordinates=\"BL\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the metric to normalize the initial 4-velocity." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n", "\\frac{4 \\, r}{r^{2} + \\cos\\left(\\mathit{th}\\right)^{2}} - 1 & 0 & 0 & -\\frac{4 \\, r \\sin\\left(\\mathit{th}\\right)^{2}}{r^{2} + \\cos\\left(\\mathit{th}\\right)^{2}} \\\\\n", "0 & \\frac{r^{2} + \\cos\\left(\\mathit{th}\\right)^{2}}{r^{2} - 4 \\, r + 1} & 0 & 0 \\\\\n", "0 & 0 & r^{2} + \\cos\\left(\\mathit{th}\\right)^{2} & 0 \\\\\n", "-\\frac{4 \\, r \\sin\\left(\\mathit{th}\\right)^{2}}{r^{2} + \\cos\\left(\\mathit{th}\\right)^{2}} & 0 & 0 & {\\left(r^{2} + \\frac{4 \\, r \\sin\\left(\\mathit{th}\\right)^{2}}{r^{2} + \\cos\\left(\\mathit{th}\\right)^{2}} + 1\\right)} \\sin\\left(\\mathit{th}\\right)^{2}\n", "\\end{array}\\right)$ " ], "text/plain": [ "[ 4*r/(r^2 + cos(th)^2) - 1 0 0 -4*r*sin(th)^2/(r^2 + cos(th)^2)]\n", "[ 0 (r^2 + cos(th)^2)/(r^2 - 4*r + 1) 0 0]\n", "[ 0 0 r^2 + cos(th)^2 0]\n", "[ -4*r*sin(th)^2/(r^2 + cos(th)^2) 0 0 (r^2 + 4*r*sin(th)^2/(r^2 + cos(th)^2) + 1)*sin(th)^2]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "g = K.metric()\n", "g[:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\tau$ is the proper time." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "tau = var('tau')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We choose a starting point $p$ and a normalized initial 4-speed $v$ in the tangent space $T_p$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "p = K((0, 14.98, pi/2, 0))\n", "Tp = K.tangent_space(p)\n", "v = Tp((2, 0, 0.005, 0.05))\n", "v = v / sqrt(-g.at(p)(v, v))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Declaration of the geodesic:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "c = K.integrated_geodesic(g, (tau, 0, 3000), v)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fast integration:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "sol = c.solve(step = 1, method=\"ode_int\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Interpolation:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "interp = c.interpolate()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An embedding in $R^3$ is needed to plot the result:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "E. = manifolds.Euclidean(3)\n", "phi = K.diff_map(E, [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting the result:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "Graphics3d Object" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "P = c.plot_integrated(mapping=phi, color=\"red\", thickness=2, plot_points=3000)\n", "P = P + sage.plot.plot3d.shapes.Sphere(4, color='grey')\n", "P.show(aspect_ratio=[1, 1, 1], viewer='threejs', online=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A few speed tests:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 248 ms, sys: 8.38 ms, total: 257 ms\n", "Wall time: 244 ms\n" ] } ], "source": [ "%time sol = c.solve(step=1, method=\"ode_int\")" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 2min 28s, sys: 597 ms, total: 2min 29s\n", "Wall time: 2min 47s\n" ] } ], "source": [ "%time sol = c.solve(step=1, method=\"rk4_maxima\")" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10000 loops, best of 3: 87.5 µs per loop\n" ] } ], "source": [ "%timeit K.default_chart().valid_coordinates(1, 1, 1, 1)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The slowest run took 13.03 times longer than the fastest. This could mean that an intermediate result is being cached.\n", "1000000 loops, best of 3: 933 ns per loop\n" ] } ], "source": [ "%timeit K.default_chart().valid_coordinates_numerical(1, 1, 1, 1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 8.3.rc1", "language": "", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.15" } }, "nbformat": 4, "nbformat_minor": 2 }