A continuación encontrarán algunos códigos útiles usados en el curso.
El cálculo numérico de los símbolos de Christoffel puede ser útil para futuros cálculos.
La rutina para calcular los símbolos de Christoffel será:
def Gamma(xmu,gfun,gargs=(),N=4,dxmax=1e-6,alpha=-1,mu=0,nu=0):
"""
Calcula todos los símbolos de Christoffel
gfun: función métrica
xmu: evento
"""
from scipy.misc import derivative
from numpy import where,arange
from numpy import zeros
#Indices
index=arange(N)
#Gamma
G=zeros((N,N,N))
for pi in range(N):
for nu in range(N):
#Inversa
gpipi=1/gfun(xmu,pi,*gargs) #g^pipi
#Coeficientes diagonales
xd=xmu[pi] #Punto en el que estoy derivando
dx=max(dxmax,dxmax*abs(xd))
gnunu_pi=derivative(lambda x:gfun(where(index==pi,x,xmu),nu,*gargs),xd,dx)
G[pi,nu,nu]=-0.5*gpipi*gnunu_pi
#Coeficientes mixtos
if nu==pi:continue
xd=xmu[nu] #Punto en el que estoy derivando
dx=max(dxmax,dxmax*abs(xd))
gpipi_nu=derivative(lambda x:gfun(where(index==nu,x,xmu),pi,*gargs),xd,dx)
G[pi,pi,nu]=0.5*gpipi*gpipi_nu
G[pi,nu,pi]=G[pi,pi,nu]
return G if alpha<0 else G[alpha,mu,nu]
Para poner a prueba nuestra rutina, podemos usar algunas de las métricas que habíamos visto antes.
Por ejemplo la métrica en coordenadas cilíndricas:
from numpy import array
def g_cilindricas_4d(xmu,mu):
"""
Coeficiente métrico g_mumu calculados en el evento xmu
para espacio-tiempo plano con coordenadas cilíndricas.
g_munu=diag(1,-1,-r^2,-1)
"""
from numpy import sin
t,r,teta,z=xmu
if mu==0:
g=1
elif mu==1:
g=-1
elif mu==2:
g=-r**2
elif mu==3:
g=-1
return g
Escogemos ahora un punto en el cuál calcular los símbolos de Christoffel:
xmu=array([0,2,0,1])
E invocamos nuestra rutina:
G=Gamma(xmu,g_cilindricas_4d,N=4)
Donde vemos que lo que habíamos determinado en nuestro primer ejercicio es correcto. Solo dos símbolos de Christoffel son diferentes de cero:
print(f"G^1_{22} = {G[1,2,2]}")
print(f"G^2_{21} = {G[2,2,1]}")
y su valor coincide con el valor teóricamente esperado de $\Gamma^1_{22}=-r$ y $\Gamma^2_{21}=1/r$.
Alguna métricas:
def g_cilindricas_2d(xmu,mu):
"""
Coeficiente métrico g_mumu calculados en el evento xmu
para espacio-tiempo plano con coordenadas cilíndricas.
g_munu=diag(1,r^2)
"""
r,teta=xmu
if mu==0:
g=1
elif mu==1:
g=r**2
return g
def g_esfera_2d(xmu,mu,R=1):
"""
Coeficiente métrico g_mumu calculados sobre la
superficie de una esfera de radio R.
g_munu=diag(1,r^2)
"""
from numpy import cos
fi,lamb=xmu
if mu==0:
g=R**2
elif mu==1:
g=R**2*cos(fi)**2
return g
from numpy import array
def g_newtoniana_4d(xmu,mu,R=1):
"""
Coeficiente métrico g_mumu calculados en el evento xmu
para espacio-tiempo plano con coordenadas cilíndricas.
g_munu=diag(A,-1,-r^2,-r^2 sin^2 teta)
"""
from numpy import sin
t,r,teta,fi=xmu
A=(1-R/r)
if mu==0:
g=A
elif mu==1:
g=-1
elif mu==2:
g=-r**2
elif mu==3:
g=-r**2*sin(teta)**2
return g
def ecuacion_geodesica_num(Y,s,gfun,gargs,N=4):
"""
Opciones:
gfun: función que da la métrica
N: Número de dimensiones
"""
from export import Gamma
from numpy import zeros
dYds=zeros(2*N)
x=Y[:N]
dxds=Y[N:]
dYds[:N]=dxds
G=Gamma(x,gfun,gargs,N)
for pi in range(N):
for mu in range(N):
for nu in range(N):
dYds[N+pi]+=-G[pi,mu,nu]*dxds[mu]*dxds[nu]
return dYds
El código a continuación fue adaptado de Hartle y sirve para el cálculo de las cantidades geométricas y tensoriales útiles en relatividad general.
(*
Instrucciones:
1) Abra una cuenta en https://lab.wolframcloud.com.
2) Cree un nuevo Notebook.
3) Copie el contenido de este archivo dentro del Notebook.
4) Ejecute el código de ejemplo
*)
(*Variables básicas*)
Clear[n,
Inversa,Matriz,Tabla,
Gama,Geodesica,Riemann,
GeodesicaTabla ,ListaDos,ListaTres,ListaCuatro];
n = 4;
(*Genericos*)
Inversa[g_]:=Simplify[Inverse[g]]
(*Muestra*)
Matriz[g_]:= g//MatrixForm
Tabla[Lista_]:=TableForm[Partition[DeleteCases[Flatten[Lista], Null], 2],TableSpacing -> {2, 2}]
GeodesicaTabla[geodesica_]:=TableForm[
Table[{"d/d\[Tau]" ToString[u[i-1]], "=",geodesica[[i]]},
{i, 1, n}], TableSpacing -> {2}]
ListaDos[Tensor_,simbolo_,coords_]:=Table[
If[UnsameQ[Tensor[[j, l]], 0],
{ToString[simbolo[coords[[j]], coords[[l]]]]<>" = ",Tensor[[j, l]]}] ,
{j, 1, n}, {l, 1, j}];
ListaTres[Tensor_,simbolo_,coords_] := Table[
If[UnsameQ[Tensor[[i, j, k]], 0],
{ToString[simbolo[coords[[i]], coords[[j]],coords[[ k]]]]<>" = ", Tensor[[i, j, k]]}] ,
{i, 1, n}, {j, 1, n}, {k, 1, j}];
ListaCuatro[Tensor_,simbolo_,coords_]:=Table[
If[UnsameQ[Tensor[[i, j, k, l]], 0],
{ToString[simbolo[coords[[i]], coords[[j]],coords[[ k]],coords[[ l]]]]<>" = ",Tensor[[i, j, k, l]]}] ,
{i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1,n}];
(*Relatividad General*)
Gama[g_,gi_,coord_] := Simplify[
Table[(1/2)*Sum[(gi[[i, s]])*
(D[g[[s, j]], coord[[k]] ] +
D[g[[s, k]], coord[[j]] ] -
D[g[[j, k]], coord[[s]] ]), {s, 1, n}],
{i, 1, n}, {j, 1, n}, {k, 1, n}] ];
Riemann[G_,coord_] := -Simplify[Table[D[G[[i, j, k]],coord[[l]] ] - D[G[[i, j, l]], coord[[k]] ] +
Sum[G[[s, j, k]]G[[i, l, s]] - G[[s, j, l]]G[[i, k, s]], {s, 1, n}],
{i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}] ]
Ricci[Rm_]:=Simplify[Table[
Sum[Rm[[i, j, i, l]], {i, 1, n}], {j, 1, n}, {l, 1, n}] ]
Scalar[gi_,Ri_]:= Simplify[ Sum[gi[[i, j]] Ri[[i, j]], {i, 1, n}, {j, 1, n}] ]
Einstein[Ri_,R_,g_] :=Simplify[Ri - (1/2) R*g]
Geodesica[G_] := Simplify[Table[-Sum[G[[i, j, k]] u[j-1] u[k-1], {j, 1, n}, {k, 1, n}], {i, 1, n}]]
(*
--------------------------------------------------------------------------
Ejemplo agujero de gusano:
coords = {T,r, \[Theta], \[Phi]};
gij = {{1,0,0,0},-{0, 1, 0, 0}, -{0,0, r^2 + b^2, 0}, -{0,0,
0, (r^2 + b^2) Sin[\[Theta]]^2}};
nij = Inversa[gij];
Gijk=Gama[gij,nij,coords];
Rijkl=Riemann[Gijk,coords];
Rij=Ricci[Rijkl];
Rs=Scalar[nij,Rij];
Gij=Einstein[Rij,Rs,gij];
dxidt=Geodesica[Gijk];
Matriz[gij]
Matriz[nij]
Tabla[ListaTres[Gijk,\[CapitalGamma],coords]]
Tabla[ListaCuatro[Rijkl,Riemann,coords]]
Tabla[ListaDos[Rij,Ricci,coords]]
Rs
Tabla[ListaDos[Gij,G,coords]]
GeodesicaTabla[dxidt]
--------------------------------------------------------------------------
Ejemplo métrica simétrica:
coords = {T,r, \[Theta], \[Phi]};
gij = {{Exp[2A[r]],0,0,0},-{0, Exp[2B[r]], 0, 0}, -{0,0, r^2 , 0}, -{0,0,
0, r^2 Sin[\[Theta]]^2}};
nij = Inversa[gij];
Gijk=Gama[gij,nij,coords];
Rijkl=Riemann[Gijk,coords];
Rij=Ricci[Rijkl];
Rs=Scalar[nij,Rij];
Gij=Einstein[Rij,Rs,gij];
dxidt=Geodesica[Gijk];
Matriz[gij]
Matriz[nij]
Tabla[ListaTres[Gijk,\[CapitalGamma],coords]]
Tabla[ListaCuatro[Rijkl,Riemann,coords]]
Tabla[ListaDos[Rij,Ricci,coords]]
Rs
Tabla[ListaDos[Gij,G,coords]]
GeodesicaTabla[dxidt]
*)