Scribd
Upload
55 views8 pages

Kerr - Solve Ivp

This document discusses numerically solving geodesic equations in Kerr spacetime using the solve_ivp function from SciPy. It defines functions for the Kerr metric components, geodesic equations of motion, and initial conditions. Initial conditions are provided and solve_ivp is used to integrate the geodesic equations from t=0 to 120, producing the trajectory of a photon around a black hole. The solution is plotted.

Uploaded by

yulieth andrea ramirez romero
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
55 views8 pages

Kerr - Solve Ivp

This document discusses numerically solving geodesic equations in Kerr spacetime using the solve_ivp function from SciPy. It defines functions for the Kerr metric components, geodesic equations of motion, and initial conditions. Initial conditions are provided and solve_ivp is used to integrate the geodesic equations from t=0 to 120, producing the trajectory of a photon around a black hole. The solution is plotted.

Uploaded by

yulieth andrea ramirez romero
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
You are on page 1/ 8

Kerr.solve_ivp 16/11/23, 5:01 a. m.

In [140… import numpy as np


import time
from timeit import timeit
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from numpy import sqrt, sin, cos, arccos, arctan, pi

In [141… def g(x):


'''
Componentes del tensor métrico

'''
Sigma = x[1]**2 + a**2*np.cos(x[2])**2
Delta = x[1]**2 -2*M*x[1] + a**2
gtt= -(1-(2*M*x[1])/Sigma)
gtphi= -(2*a*M*x[1]*np.sin(x[2])**2)/Sigma
grr= Sigma/Delta
gthth= Sigma
gphph= (x[1]**2 + a**2 + (2*M*x[1]*(a*np.sin(x[2]))**2)/Sigma)*np.sin(x[

return [gtt, gtphi, grr, gthth, gphph]

http://localhost:8888/nbconvert/html/Desktop/tesis/Kerr/Kerr.solve_ivp%20.ipynb?download=false Página 1 de 8
Kerr.solve_ivp 16/11/23, 5:01 a. m.

In [142… def geodesics(tau, q):

'''
Ecuaciones geodéicas

'''
Sigma = q[1]**2 + a**2*np.cos(q[2])**2
Delta = q[1]**2 -2*M*q[1] + a**2

A=(q[1]**2+a**2)**2-Delta*a**2*np.sin(q[2])**2

drA=4*q[1]*(a**2+q[1]**2)-2*(q[1]-M)*a**2*np.sin(q[2])**2
drS=2*q[1]
drD=2*(q[1]-M)
d0A=-2*Delta*a**2*np.sin(q[2])*np.cos(q[2])
d0S=-2*a**2*np.cos(q[2])*np.sin(q[2])
d0D=0
drg00=1/(Sigma*Delta)*(-drA+A/Sigma*drS+A/Delta*drD)
drg03=(2*a*M)/(Sigma*Delta)*(-1+q[1]/Sigma*drS+q[1]/Delta*drD)
drg11=1/Sigma*(drD-Delta/Sigma*drS)
drg22=-1/Sigma**2*drS
drg33=1/(Delta*Sigma*np.sin(q[2])**2)*((a**2*np.sin(q[2])**2)/Delta*drD-
d0g00=1/(Sigma*Delta)*(-d0A+(A/Sigma)*d0S)
d0g03=(2*a*M*q[1])/(Sigma**2*Delta)*d0S
d0g11=-1/Sigma**2*d0A
d0g22=-1/Sigma**2*d0S
d0g33=-1/(Delta*Sigma*np.sin(q[2])**2)*((2*Delta*np.cos(q[2]))/np.sin(q[
dttau =-A/(Sigma*Delta)*q[4]-(2*a*M*q[1])/(Sigma*Delta)*q[7]
drdtau = Delta/Sigma*q[5]
dthdtau = 1/Sigma*q[6]
dphidtau = (2*a*M*q[1])/(Sigma*Delta)*q[4]+(Delta-a**2*np.sin(q[2])**2)/
dk_tdtau = 0.
dk_rdtau = -1/2*drg00*q[4]**2-drg03*q[4]*q[7]-1/2*drg11*q[5]**2-1/2*drg22
dk_thdtau = -1/2*d0g00*q[4]**2-d0g03*q[4]*q[7]-1/2*d0g11*q[5]**2-1/2*d0g22
dk_phidtau = 0.
return [dttau, drdtau, dthdtau, dphidtau, dk_tdtau, dk_rdtau, dk_thdtau,

http://localhost:8888/nbconvert/html/Desktop/tesis/Kerr/Kerr.solve_ivp%20.ipynb?download=false Página 2 de 8
Kerr.solve_ivp 16/11/23, 5:01 a. m.

In [143… def initCond(x, k, metric, M):


'''
Condiciones iniciales y componente kt para bajar indices
'''

kr, kth, kphi= k


g=metric(x)

# Metric components
g_tt, g_tphi, g_rr, g_thth, g_phph = metric(x)

kt=(-g_tphi*kphi+np.sqrt((g_tphi*kphi)**2-g_tt*(g_rr*(kr)**2+g_thth*(kth

# Lower k-indices
k_t = g_tt*kt+g_tphi*kphi
k_r = g_rr*kr
k_th = g_thth*kth
k_phi = g_phph*kphi+g_tphi*kt

return [x[0], x[1], x[2], x[3], k_t, k_r, k_th, k_phi]

In [144… M = 1
a=0.5

# Initial Conditions
t0 = 0.
r0 = 7
theta0 = np.pi/3.5
phi0 = 0.

kr0 = 0.5
kth0 = 0.
kphi0 = 0.5

x = [t0, r0, theta0, phi0]


k = [kr0, kth0,kphi0]

ic = initCond(x, k, g, a)
lmbda_span=(0.0, 120.0)

sol = solve_ivp(geodesics, lmbda_span, ic)


sol.t.shape

(13,)
Out[144]:

In [151… type(sol)

scipy.integrate._ivp.ivp.OdeResult
Out[151]:

http://localhost:8888/nbconvert/html/Desktop/tesis/Kerr/Kerr.solve_ivp%20.ipynb?download=false Página 3 de 8
Kerr.solve_ivp 16/11/23, 5:01 a. m.

In [152… print(sol)

message: 'The solver successfully reached the end of the integration inter
val.'
nfev: 74
njev: 0
nlu: 0
sol: None
status: 0
success: True
t: array([0.00000000e+00, 3.12677971e-04, 3.43945768e-03, 3.47072548
e-02,
3.47385226e-01, 1.45491771e+00, 3.07992129e+00, 5.59918516e+00,
1.05745621e+01, 2.22299465e+01, 4.88446070e+01, 1.16551653e+02,
1.20000000e+02])
t_events: None
y: array([[ 0.00000000e+00, -1.05742795e-03, -1.16311214e-02,
-1.17308812e-01, -1.16756628e+00, -4.76943311e+00,
-9.70865663e+00, -1.68054088e+01, -2.99081284e+01,
-5.91063313e+01, -1.23893281e+02, -2.86383349e+02,
-2.94628235e+02],
[ 7.00000000e+00, 7.00015637e+00, 7.00172324e+00,
7.01771013e+00, 7.20879582e+00, 8.29647276e+00,
1.07474441e+01, 1.56008327e+01, 2.65038011e+01,
5.35419483e+01, 1.16393181e+02, 2.76971737e+02,
2.85155182e+02],
[ 8.97597901e-01, 8.97597907e-01, 8.97598657e-01,
8.97674485e-01, 9.04861092e-01, 9.93035676e-01,
1.15509153e+00, 1.32170858e+00, 1.46600145e+00,
1.56662020e+00, 1.61895743e+00, 1.64463642e+00,
1.64516831e+00],
[ 0.00000000e+00, 1.62339903e-04, 1.78533663e-03,
1.79739857e-02, 1.74736510e-01, 6.22050809e-01,
9.97421005e-01, 1.26252201e+00, 1.45663701e+00,
1.58440137e+00, 1.65032947e+00, 1.68281619e+00,
1.68349090e+00],
[ 2.37395220e+00, 2.37395220e+00, 2.37395220e+00,
2.37395220e+00, 2.37395220e+00, 2.37395220e+00,
2.37395220e+00, 2.37395220e+00, 2.37395220e+00,
2.37395220e+00, 2.37395220e+00, 2.37395220e+00,
2.37395220e+00],
[ 6.96413970e-01, 6.96666172e-01, 6.99187171e-01,
7.24294275e-01, 9.64163387e-01, 1.62249255e+00,
2.11888289e+00, 2.37026919e+00, 2.44996620e+00,
2.43777741e+00, 2.40967457e+00, 2.39037323e+00,
2.38993429e+00],
[ 0.00000000e+00, 1.96187646e-03, 2.15757909e-02,
2.17206526e-01, 2.10139912e+00, 7.05368344e+00,
1.01937685e+01, 1.15972619e+01, 1.21281427e+01,
1.22347857e+01, 1.22127477e+01, 1.21825019e+01,
1.21817374e+01],
[ 1.53603362e+01, 1.53603362e+01, 1.53603362e+01,
1.53603362e+01, 1.53603362e+01, 1.53603362e+01,

http://localhost:8888/nbconvert/html/Desktop/tesis/Kerr/Kerr.solve_ivp%20.ipynb?download=false Página 4 de 8
Kerr.solve_ivp 16/11/23, 5:01 a. m.

1.53603362e+01, 1.53603362e+01, 1.53603362e+01,


1.53603362e+01, 1.53603362e+01, 1.53603362e+01,
1.53603362e+01]])
y_events: None

In [145… x=sol.y[1]*np.cos(sol.y[3])
y=sol.y[1]*np.sin(sol.y[3])
BH = plt.Circle((0, 0), 2*M, color='k')

fig, ax = plt.subplots(figsize=(5,5))
ax.plot(x,y, color='crimson', label='Fotón')
ax.add_patch(BH)
ax.set_xlim(-20,20)
ax.set_ylim(-20,20)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.grid(alpha=0.2)
ax.axvline(0, c='k', alpha=0.1)
ax.axhline(0, c='k', alpha=0.1)
plt.legend()
plt.show()

http://localhost:8888/nbconvert/html/Desktop/tesis/Kerr/Kerr.solve_ivp%20.ipynb?download=false Página 5 de 8
Kerr.solve_ivp 16/11/23, 5:01 a. m.

In [146… def plot3D(sol):


x=sol.y[1]*np.sin(sol.y[2])*np.cos(sol.y[3])
y=sol.y[1]*np.sin(sol.y[2])*np.sin(sol.y[3])
z=sol.y[1]*np.cos(sol.y[2])

u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
xs = 2*M*np.cos(u)*np.sin(v)
ys = 2*M*np.sin(u)*np.sin(v)
zs = 2*M*np.cos(v)

ax = plt.figure().add_subplot(projection='3d')
ax.plot_surface(xs, ys, zs, color='k')
ax.plot(x, y, z, color='crimson')
ax.set_xlim(-20,20)
ax.set_ylim(-20,20)
ax.set_zlim(-15,15)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$z$')
plt.show()

In [147… plot3D(sol)

In [148… sol.y

http://localhost:8888/nbconvert/html/Desktop/tesis/Kerr/Kerr.solve_ivp%20.ipynb?download=false Página 6 de 8
Kerr.solve_ivp 16/11/23, 5:01 a. m.

array([[ 0.00000000e+00, -1.05742795e-03, -1.16311214e-02,


Out[148]:
-1.17308812e-01, -1.16756628e+00, -4.76943311e+00,
-9.70865663e+00, -1.68054088e+01, -2.99081284e+01,
-5.91063313e+01, -1.23893281e+02, -2.86383349e+02,
-2.94628235e+02],
[ 7.00000000e+00, 7.00015637e+00, 7.00172324e+00,
7.01771013e+00, 7.20879582e+00, 8.29647276e+00,
1.07474441e+01, 1.56008327e+01, 2.65038011e+01,
5.35419483e+01, 1.16393181e+02, 2.76971737e+02,
2.85155182e+02],
[ 8.97597901e-01, 8.97597907e-01, 8.97598657e-01,
8.97674485e-01, 9.04861092e-01, 9.93035676e-01,
1.15509153e+00, 1.32170858e+00, 1.46600145e+00,
1.56662020e+00, 1.61895743e+00, 1.64463642e+00,
1.64516831e+00],
[ 0.00000000e+00, 1.62339903e-04, 1.78533663e-03,
1.79739857e-02, 1.74736510e-01, 6.22050809e-01,
9.97421005e-01, 1.26252201e+00, 1.45663701e+00,
1.58440137e+00, 1.65032947e+00, 1.68281619e+00,
1.68349090e+00],
[ 2.37395220e+00, 2.37395220e+00, 2.37395220e+00,
2.37395220e+00, 2.37395220e+00, 2.37395220e+00,
2.37395220e+00, 2.37395220e+00, 2.37395220e+00,
2.37395220e+00, 2.37395220e+00, 2.37395220e+00,
2.37395220e+00],
[ 6.96413970e-01, 6.96666172e-01, 6.99187171e-01,
7.24294275e-01, 9.64163387e-01, 1.62249255e+00,
2.11888289e+00, 2.37026919e+00, 2.44996620e+00,
2.43777741e+00, 2.40967457e+00, 2.39037323e+00,
2.38993429e+00],
[ 0.00000000e+00, 1.96187646e-03, 2.15757909e-02,
2.17206526e-01, 2.10139912e+00, 7.05368344e+00,
1.01937685e+01, 1.15972619e+01, 1.21281427e+01,
1.22347857e+01, 1.22127477e+01, 1.21825019e+01,
1.21817374e+01],
[ 1.53603362e+01, 1.53603362e+01, 1.53603362e+01,
1.53603362e+01, 1.53603362e+01, 1.53603362e+01,
1.53603362e+01, 1.53603362e+01, 1.53603362e+01,
1.53603362e+01, 1.53603362e+01, 1.53603362e+01,
1.53603362e+01]])

In [149… sol1=sol.y[:,1]
sol1

array([-1.05742795e-03, 7.00015637e+00, 8.97597907e-01, 1.62339903e-04,


Out[149]:
2.37395220e+00, 6.96666172e-01, 1.96187646e-03, 1.53603362e+01])

In [150… sol2=sol.y[:,2]
sol2

array([-1.16311214e-02, 7.00172324e+00, 8.97598657e-01, 1.78533663e-03,


Out[150]:
2.37395220e+00, 6.99187171e-01, 2.15757909e-02, 1.53603362e+01])

http://localhost:8888/nbconvert/html/Desktop/tesis/Kerr/Kerr.solve_ivp%20.ipynb?download=false Página 7 de 8
Kerr.solve_ivp 16/11/23, 5:01 a. m.

In [139… plt.plot(sol1,sol2)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Lab DLS')
plt.show()

In [ ]:

In [ ]:

http://localhost:8888/nbconvert/html/Desktop/tesis/Kerr/Kerr.solve_ivp%20.ipynb?download=false Página 8 de 8

You might also like