Indice | Previo: RelatividadGeneral.SimbolosChristoffel | Siguiente: RelatividadGeneral.InerciaYGeodesicas

2.4. Derivada total general y geodésicas

Una aplicación interesante de la derivada covariante es la deducción de una de las ecuaciones más importantes de la relatividad general: la ecuación geodésica.

2.4.1. Derivada total general

Hemos introducido ya una forma de la derivada que puede ser usada en las leyes de la física de modo que tengan una forma manifiestamente covariante bajo trasnformaciones generales. Es la derivada covariante:

$$ {A^\mu}_{;\alpha} = {A^\mu}_{,\alpha} + \Gamma^\mu_{\gamma\alpha} A^\gamma $$

Esta derivada permite definir otra cantidad muy útil, a saber la derivada total general, es decir, la derivada total de una cantidad tensorial respecto a un escalar (por ejemplo, el tiempo propio). Como hemos visto en el capítulo anterior, las derivadas totales son muy importantes en la definición de las propiedades cinemáticas y dinámicas en física, de modo que es conveniente definirlas.

Si nos movemos a lo largo de una curva en el espacio-tiempo tal que cada evento es función de una cantidad escalar o parámetro $u$, es decir $x^\mu(u)$, entonces definiremos la derivada total de un campo tensorial $\tilde A(\tilde x)$ respecto al parámetro $u$ como:

$$ \frac{\mathrm{D}}{\mathrm{D}u}A^\mu \equiv {A^\mu}_{;\alpha}\frac{\mathrm{d} x^\alpha}{\mathrm{d}u} $$

que es un análogo a la regla de la cadena en el cálculo vectorial convencional.

Nótese que esta definición es análoga a la que podemos hacer en el caso en el que tengamos transformaciones lineales de coordenadas (métrica de Minkowski), a saber $\mathrm{d}A^\mu/\mathrm{d}u={A^\mu}_{,\alpha}\mathrm{d} x^\beta/\mathrm{d}u$ con la diferencia de que en lugar de usar la derivada "," usamos la derivada covariante. De nuevo es sencillo probar que la derivada total general es en sí misma una cantidad tensorial (ver Problemas al final del capítulo).

2.4.2. Derivada total y transporte paralelo

En términos de la derivada total es posible reinterpretar el concepto de transporte paralelo.

Diremos que un vector es transportado de forma paralela a lo largo de una trayectoria si su derivada total general es proporcional es cero:

$$ \frac{\mathrm{D}}{\mathrm{D}u}A^\mu =0\;\mathrm{, transporte\;paralelo} $$

Usando la definición de derivada total general:

$$ \partial_\alpha A^\mu \frac{\mathrm{d} x^\alpha}{\mathrm{d}u} + \Gamma^\mu_{\alpha\gamma} A^\gamma \frac{\mathrm{d} x^\alpha}{\mathrm{d}u}=0 $$

El primer término es la derivada total convencional:

$$ \frac{\mathrm{d} A^\mu}{\mathrm{d}u} + \Gamma^\mu_{\alpha\gamma} A^\gamma \frac{\mathrm{d} x^\alpha}{\mathrm{d}u}=0 $$

Esta ecuación en definitiva representa la condición matemática para que un vector sea transportado de forma paralela a lo largo de una trayectoria descrita por $\tilde x(u)$.

2.4.3. Transporte paralelo y bases vectoriales

Existe una manera alternativa de obtener la fórmula de la derivada covariante y de las conexiones que es interesante. Un cuadrivector puede escribirse en términos de los vectores de una base como:

$$\mathbf{A}= A^\mu \hat{e}_\mu$$

Si derivamos respecto de un parámetro $u$ que define una trayectoria obtenemos:

$$\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}u}=\frac{\mathrm{d}A^\mu}{\mathrm{d}u}\hat{e}_\mu+A^\mu\frac{\mathrm{d}\hat{e}_\mu}{\mathrm{d}u}$$

Las derivadas de cada vector de la base respecto a las coordenadas se pueden expresar en términos de los mismos vectores de la base:

\begin{equation} \label{eq:derivada_e} \partial_\alpha\hat{e}_\mu\equiv G^\gamma_{\alpha\mu}\hat{e}_\gamma \end{equation}

donde los coeficientes $G^\gamma_{\alpha\mu}$ conectan las derivadas de cada componente de cada vector unitario con las componentes de los demás vectores unitarios de la base.

Reemplazando queda:

$$\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}u}=\frac{\mathrm{d}A^\mu}{\mathrm{d}u}\hat{e}_\mu+A^\mu G^\gamma_{\alpha\mu}\hat{e}_\gamma\frac{\mathrm{d}x^\alpha}{\mathrm{d}u}$$

que se puede escribir como: $$\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}u}=\left(\frac{\mathrm{d}A^\gamma}{\mathrm{d}u}+G^\gamma_{\alpha\mu}A^\mu \frac{\mathrm{d}x^\alpha}{\mathrm{d}u}\right)\hat{e}_\gamma$$

Para que el vector no varíe a lo largo de la trayectoria es necesario que $\mathrm{d}\mathbf{A}/{\mathrm{d}u}=0$ lo que implica que:

$$ \frac{\mathrm{d}A^\gamma}{\mathrm{d}u}+G^\gamma_{\alpha\mu}A^\mu \frac{\mathrm{d}x^\alpha}{\mathrm{d}u} = 0$$

que es justamente la condición que habíamos obtenido antes.

Bastaría ahora con probar que $G^\gamma_{\alpha\mu}$ son los símbolos de Christoffel. Para ello debemos primero reconocer que la métrica se define en términos del producto punto entre los vectores de la base:

\begin{equation} \label{eq:g_edote} g_{\mu\nu}=\hat{e}_\mu\cdot\hat{e}_\nu \end{equation}

Podemos calcular la derivada convencional de los coeficientes métricos:

$$ \begin{array}{rcl} \partial_\alpha g_{\mu\nu} & = & G^\gamma_{\alpha\mu}\hat{e}_\gamma\cdot\hat{e}_\nu+G^\gamma_{\nu\beta}\hat{e}_\mu\cdot\hat{e}_\gamma\\ & = & G^\gamma_{\alpha\mu}g_{\gamma\nu}+G^\gamma_{\nu\beta}g_{\mu\gamma} \end{array} $$

Si se repite el mismo procedimiento cambiando cíclicamente los índices de $\partial_\alpha g_{\mu\nu}$ se obtiene:

$$ G^\lambda_{\mu\nu}=\frac{1}{2}g^{\lambda\rho}(g_{\rho\mu,\nu}-g_{\mu\nu,\rho}+g_{\nu\rho,\mu}) $$

que es justamente la forma de los símbolos de Christoffel.

2.4.4. Bases y componentes vectoriales en espacio curvo

Las Ecuaciones (derivada_e) y (g_edote) ponen de relieve una importante propiedad que tienen los vectores en el espacio curvo.

Considere por ejemplo el caso en el que tenemos una métrica diagonal. En esta situación la Ec. (g_edote) se puede escribir como:

$$ |\hat{e}_\mu|=\sqrt{|g_{\mu\mu}|} $$

Es decir, cuando expresamos un vector en un espacio general los vectores de la base no son necesariamente unitarios: ¡esta si que es una sorpresa! Pero la sorpresa va más allá: en realidad este resultado no solo es válido en espacio curvo sino también en el espacio plano.

Consideremos por ejemplo el caso de la descripción del espacio plano de 3 dimensiones en coordenadas cilíndricas $(r,\theta,z)$, $g_{\mu\nu}=\mathrm{diag}(1,r^2,1)$. En este caso la magnitud de los vectores unitarios será:

$$ \begin{array}{rcl} |\hat{e}_r| & = & \sqrt{|g_{rr}|}=1\\ |\hat{e}_\theta| & = & \sqrt{|g_{\theta\theta}|}=r\\ |\hat{e}_z| & = & \sqrt{|g_{zz}|}=1 \end{array} $$

Todo lo que hemos dicho sobre transporte paralelo y derivada covariante aplica para las componentes de vectores definidos sobre esta base, por ejemplo en el caso del espacio euclidiano de 3 dimensiones y en coordenadas cilíndircas:

$$ \vec{A}=A^r\hat{e}_r+A^\theta\hat{e}_\theta+A^z\hat{e}_z $$

En el cálculo (y por lo tanto en la física, incluyendo en la relatividad especial), por otro lado hemos estado usando bases ortonormales que llamaremos en lo sucesivo, para diferenciarlas de estas bases curvas, como $\hat{u}_\mu$. En estas bases las componentes del mismo vector mencionado anteriormente se escriben como:

$$ \vec{A}=A^{\hat{r}}\hat{u}_r+A^{\hat\theta}\hat{u}_\theta+A^{\hat z}\hat{u}_z $$

donde nótese que para distinguirlas de las componentes en la base curva hemos escrito una barra sobre los indices respectivos.

¿Cómo se relacionan las componentes de un vector en la base curva con las componentes en la base ortonormal? En el caso estudiado aquí, en donde ambas bases son ortogonales y los vectores unitarios de cada una solo se diferencian por su longitud, las relaciones serán simplemente:

$$ \begin{array}{rcl} \hat{e}_r & = & \sqrt{g_{rr}}\hat{u}_r=\hat{u}_r\\ \hat{e}_\theta & = & \sqrt{g_{\theta}}\hat{u}_\theta=r\hat{u}_\theta\\ \hat{e}_z & = & \sqrt{g_{z}}\hat{u}_z=\hat{u}_z \end{array} $$

y las componentes del vector serán:

$$ \begin{array}{rcl} A^r & = & A^{\hat r}\\ A^\theta & = & r A^{\hat\theta}\\ A^z & = & A^{\hat z} \end{array} $$

Es interesante anotar que físicamente estas componentes se diferencian de manera importante. Las componentes que medimos en los laboratorios para las cantidades físicas son las componentes en la base ortonormal. Por tanto es siempre importante cuando trabajamos con métricas arbitrarias tener siempre presente el tipo de transformaciones definidas anteriormente para pasar de las componentes calculadas (normalmente las componentes en la base curva) a las componentes observadas (componentes en la base ortonormal).

Para ponerle orden a todo esto y generalizarlo al espacio-tiempo introducimos las siguientes definición:

Definición: base coordenada y base física. En cualquier punto de un espacio-tiempo descrito por una métrica $g_\mu\nu$ en un sistema de coordenadas dado siempre es posible construir dos bases vectoriales y respecto a ellas definir las componentes de los cuadrivectores de interés:

• La base física, formada por un conjunto de vectores ortonormales $\hat{u}_\mu$ que cumplen: $$\hat{u}_{\mu}\cdot\hat{u}_\nu=\eta_{\mu\nu}$$ Las componentes de un cuadrivector arbitrario en la base física se denotan como $A^{\bar\mu}$ y satisfacen: $$\tilde{A}=A^{\hat\mu}\hat{u}_\mu$$ Nuestros experimentos miden el valor de $A^{\hat\mu}$.

• La base coordenada o base curva, formada por un conjunto de vectores no necesariamente ortogonales y no necesariamente unitarios $\hat{e}_\mu$ que cumplen: $$\hat{e}_{\mu}\cdot\hat{e}_\nu=g_{\mu\nu}$$ Las componentes de un cuadrivector arbitrario en la base coordenada se denotan como $A^{\mu}$ y satisfacen: $$\tilde{A}=A^{\mu}\hat{e}_\mu$$ Nuestros cálculos en relatividad general, permiten calcular el valor de $A^{\mu}$.

Cada uno de los vectores de la base coordenada se puede expresar en términos de la base física: $$\hat{e}_\mu=e^\alpha_\mu \hat{u}_\alpha$$ y por lo tanto las componentes físicas de un vector se pueden escribir en términos de sus componentes coordenadas como: $$A^{\hat\alpha}=A^{\mu} e^\alpha_\mu$$

Algunas situaciones prácticas del uso de la base coordenada y la base física se muestran en los ejemplos a continuación.

2.4.5. Ejemplos numéricos de transporte paralelo

Vimos entonces que las componentes coordenadas de un vector transfortado de forma paralela a lo largo de una trayectoria $x^\mu(u)$ esta dado por:

$$ \frac{\mathrm{d} A^\mu}{\mathrm{d}u} + \Gamma^\mu_{\alpha\gamma} A^\gamma \frac{\mathrm{d} x^\alpha}{\mathrm{d}u}=0 $$

Esta ecuación, en general, corresponde a un conjunto de 4 (o $n$ en un espacio de $n$ dimensiones) ecuaciones diferenciales de primer orden. Cada ecuación contienen en el lado izquierdo un todal de (en general) 1+16=17 términos (o $1+n^2$ términos):

$$ \frac{\mathrm{d} A^\mu}{\mathrm{d}u} + \Gamma^\mu_{00} A^0 \frac{\mathrm{d} x^0}{\mathrm{d}u} + \Gamma^\mu_{01} A^0 \frac{\mathrm{d} x^1}{\mathrm{d}u}+\ldots+ \Gamma^\mu_{10} A^1 \frac{\mathrm{d} x^0}{\mathrm{d}u} + \Gamma^\mu_{11} A^1 \frac{\mathrm{d} x^1}{\mathrm{d}u}+\ldots=0 $$

Naturalmente las simetrías de la métrica y de los símbolos de Christoffel hacen que estas ecuaciones sean mucho más cortas.

Implementemos estas ecuaciones en una rutina que poamos integrar numéricamente:

Usando la rutina simbólica del cálculo de los símbolos de Christoffel:

def A_parallel(A,u,xfun,Gfun,fargs=(),N=4):
    """
    Calcula la derivada de las componentes de un vector A 
    respecto al parámetro u de una función xfun
    
    Parametros:
        A: Arreglo con componentes coordenadas del vector
        u: Valor del parámetro
    
    Opciones:
        xfun: función de la trayectoria (posición y derivada)
        Gfun: función que da los símbolos de Christoffel
        fargs: argumentos de la función de la trayectoria
        N: Número de dimensiones
    """
    from numpy import zeros,array
    dAdu=zeros(N)
    xmu,dxmudu=xfun(u,*fargs)
    G=array(Gfun(*xmu))
    for pi in range(N):
        for mu in range(N):
            for nu in range(N):
                dAdu[pi]+=-G[pi][mu][nu]*A[mu]*dxmudu[nu]
    return dAdu

2.4.5.1. Transporte paralelo en coordenadas cilíndrica

Para aplicar la anterior ecuación consideremos por ejemplo el transporte paralelo del siguiente vector en el espacio euclidiano de 2 dimensiones y en coordenadas cilíndricas.

N=2
from numpy import array
A0=array([1,0])

Podemos transportar el vector a través de una gran familia de trayectorias, pero debemos cuidar para hacerlo que los símbolos de Christoffel no sean singulares. Como vimos en la sección anterior $\Gamma^{\theta}_{r\theta}=1/r$ de modo que debemos evitar las trayectorias que pasan por el origen.

En la siguiente rutina construimos un conjunto de trayectorias que tienen esa propiedad.

def x_fun(u,tipo="circunferencia"):
    """
    Esta función define nuestro camino en el espacio
    """
    from numpy import zeros
    x=zeros(2)
    dxdu=zeros(2)
    
    if tipo=="circunferencia":
        x[0]=1;dxdu[0]=0;
        x[1]=u;dxdu[1]=1;
    if tipo=="cicloide":
        x[0]=2+cos(u);dxdu[0]=-sin(u);
        x[1]=u;dxdu[1]=1;
    if tipo=="espiral":
        x[0]=1+u;dxdu[0]=1;
        x[1]=u;dxdu[1]=1;
    if tipo=="elipse":
        x[0]=1/(1+0.5*cos(u));dxdu[0]=0.5*sin(u)/(1+0.5*cos(u))**2;
        x[1]=u;dxdu[1]=1;
    return x,dxdu

Para realizar el transporte paralelo necesitamos por otro lado la función que define la métrica de este espacio que es simplemente:

$$ g_{\mu\nu}=\mathrm{diag}(1,r^2) $$

El cálculo simbólico nos da:

def Gfun_cil2d():
    from export import Gamma_sym
    from sympy import symbols,diag
    s=symbols('r,theta')
    r,q=s
    gcomp=diag(1,r**2)
    Gfun=Gamma_sym(gcomp,s)
    return Gfun

Ahora podemos resolver la ecuación del transporte paralelo:

from numpy import pi,linspace
us=linspace(0,2*pi,20)

from scipy.integrate import odeint
tipo="circunferencia"
#tipo="cicloide"
#tipo="espiral"
#tipo="elipse"
Gfun=Gfun_cil2d()
As=odeint(A_parallel,A0,us,args=(x_fun,Gfun,(tipo,),N))

A = 
[[ 1.     0.   ]
 [ 0.946 -0.325]
 [ 0.789 -0.614]
 ...
 [ 0.789  0.614]
 [ 0.946  0.325]
 [ 1.     0.   ]]

¿Qué son estos números? Estas son las componentes coordenadas del vector $A$ a lo largo de la trayectoria. Recordemos que por componentes coordenadas entendemos las componentes en la base coordenada (que no es de vectores unitarios).

Si queremos representar graficamente el transporte paralelo podemos hacerlo en nuestro familiar sistema de coordenadas cartesianas (en el que funciona justamente el sistema de graficación de Python). Para ello recordemos que la base ortonormal se escribe en términos de los vectores cartesianos como:

$$ \begin{array}{rcl} \hat{u}_r & = & \cos\theta\;\hat{u}_x+\sin\theta\;\hat{u}_y\\ \hat{u}_\theta & = & -\sin\theta\;\hat{u}_x+\cos\theta\;\hat{u}_y\\ \end{array} $$

Por lo tanto la base coordenada será:

$$ \begin{array}{rcl} \hat{e}_r & = & \cos\theta\;\hat{u}_x+\sin\theta\;\hat{u}_y\\ \hat{e}_\theta & = & -r\sin\theta\;\hat{u}_x+r\cos\theta\;\hat{u}_y\\ \end{array} $$

Es justamente respecto a esta última base que se calcularon las componentes del vector transportado de forma paralela.

Un gráfico del vector transportado se muestra en la Figura (code:transporte_paralelo). En la versión electrónica del libro pueden encontrar el código usado para generar esta figura.

import warnings
warnings.filterwarnings('ignore')
%matplotlib inline

Figura 2.35. Vector transportado de forma paralela en coordenadas cilíndricas

2.4.5.2. Transporte paralelo sobre una esfera

Un caso menos trivial de transporte paralelo es el que se produce en la superficie de una esfera (superficie curva de dos dimensiones), en el que los puntos vienen dados en función de su latitud $\phi$ y longitud $\lambda$. Esta superficie tiene métrica:

$$ \mathrm{d}l^2=R^2\mathrm{d}\phi^2+R^2\cos^2\phi\;\mathrm{d}\lambda^2 $$

con coeficientes métricos:

$$ g_{ij}:\mathrm{diag}(R^2,R^2\cos^2\phi) $$

def Gfun_esf2d(R):
    from export import Gamma_sym
    from sympy import symbols,diag,cos
    s=symbols('f,l')
    f,l=s
    gcomp=diag(R**2,R**2*cos(f)**2)
    Gfun=Gamma_sym(gcomp,s)
    return Gfun

Supongamos que queremos transportar de forma paralela el vector que apunta directamente hacia el norte:

N=2
from numpy import array
A0=array([0,1])

Podemos seguir distintas trayectorias. Por ejemplo ir desde el ecuador hasta el polo siguiendo un meridiano ($\lambda=\mathrm{cte}$). O podríamos ir alrededor de la esfera sobre un paralelo ($\phi=\mathrm{cte}$). Podemos implementar estas trayectorias con esta rutina:

def x_fun_esfera(u,lat_0=0,lon_0=0,tipo="meridiano"):
    """
    Esta función define nuestro camino en el espacio
    """
    from numpy import zeros
    x=zeros(2)
    dxdu=zeros(2)
    if tipo=="meridiano":
        x[0]=u;dxdu[0]=1;
        x[1]=lon_0;dxdu[1]=0;
    if tipo=="paralelo":
        x[0]=lat_0;dxdu[0]=0;
        x[1]=u;dxdu[1]=1;
    if tipo=="paralelo_contrario":
        x[0]=lat_0;dxdu[0]=0;
        x[1]=-u;dxdu[1]=-1;
    return x,dxdu

Una integral de la ecuación de transporte paralelo se puede obtener usando:

from numpy import pi,linspace
us=linspace(0,360*pi/180,10)

from scipy.integrate import odeint
R=1
Gfun=Gfun_esf2d(R)
#tipo="meridiano"
tipo="paralelo"
#tipo="paralelo_contrario"
lat_0=45*pi/180
lon_0=0
As=odeint(A_parallel,A0,us,args=(x_fun_esfera,Gfun,(lat_0,lon_0,tipo),N))

A = 
[[ 0.     1.   ]
 [-0.335  0.881]
 [-0.59   0.551]
 ...
 [ 0.218 -0.951]
 [ 0.511 -0.691]
 [ 0.682 -0.266]]

Para representar este vector realizaremos una proyección ortográfica sobre el plano $x-y$ usando la regla:

$$ x=\theta\cos\lambda\\ y=\theta\sin\lambda $$

donde $\theta=\pi/2-\phi$ es la colatitud.

De nuevo el resultado obtenido esta expresado en la base coordenada sobre la esfera. Podemos expresar esta base en coordenadas cartesianas si hacemos primero una proyección ortonormal en el sistema de coorenadas cartesianas así:

$$ \begin{array}{rcl} \hat{u}_\phi & = & -\cos\lambda\;\hat{u}_x-\sin\lambda\;\hat{u}_y\\ \hat{u}_\lambda & = & -\sin\lambda\;\hat{u}_x+\cos\lambda\;\hat{u}_y\\ \end{array} $$

De allí los vectores coordenados en coordenadas esféricas serán:

$$ \begin{array}{rcl} \hat{e}_\phi & = & -R\cos\lambda\;\hat{u}_x-R\sin\lambda\;\hat{u}_y\\ \hat{e}_\lambda & = & -R\cos\theta\sin\lambda\;\hat{u}_x+R\cos\theta\cos\lambda\;\hat{u}_y\\ \end{array} $$

Figura 2.36. Vector transportado de forma paralela sobre la superficie de una esfera. Se usa proyección azimuthal para representar las coordenadas (malla punteada).

Figura 2.37.

2.4.6. Ecuación geodésica

De todos los campos vectoriales que pueden definirse a lo largo de una trayectoria el más importante es aquel que corresponde al vector tangente (o el vector cuadrivelocidad). El transporte paralelo de este vector permite definir un tipo de trayectoria muy especial.

Definición 4.21. Geodésica. Una trayectoria en un espacio general se define como una geodésica si mantiene a todo lo largo la misma dirección, es decir si el vector tangente $t^\alpha\equiv \mathrm{d}x^\alpha/\mathrm{d}\sigma$ en cada punto es paralelo (en el sentido de transporte paralelo) al vector tangente de cualquier otro punto. Matemáticamente:

$$\frac{\mathrm{D}}{\mathrm{D}\sigma}t^\alpha=0$$

donde $\sigma$ se conoce como el parámetro afin de la trayectoria.

Si usamos la definición dada en la sección anterior, la geodésica será la trayectoria que satisfaga la ecuación:

$$ \frac{\mathrm{d^2} x^\mu}{\mathrm{d}\sigma^2} + \Gamma^\mu_{\alpha\gamma} \frac{\mathrm{d} x^\alpha}{\mathrm{d}\sigma}\frac{\mathrm{d} x^\gamma}{\mathrm{d}\sigma}=0 $$

Figura 2.38. Definición de geodésica en el espacio-tiempo plano y sobre la superficie de una esfera.

2.4.7. Ejemplos numéricos de geodésicas

2.4.7.1. Geodésica en coordenadas cilíndricas

Nos proponemos ahora a calcular la geodésica que sigue un cuerpo en el espacio plano dada una condición inicial. La ecuación de la geodésica es:

\begin{equation} \frac{\mathrm{d^2} x^\mu}{\mathrm{d}\sigma^2} + \Gamma^\mu_{\alpha\gamma} \frac{\mathrm{d} x^\alpha}{\mathrm{d}\sigma}\frac{\mathrm{d} x^\gamma}{\mathrm{d}\sigma}=0 \end{equation}

Para programar la ecuación de la geodésica es necesario linearizarla y expresarla de forma general como:

$$ \left\{\frac{\mathrm{d}Y^\mu}{\mathrm{d}\sigma}=f^\mu(\{Y^\nu\},u)\right\} $$

En este caso podemos hacer esta asignación:

$$ \begin{array}{rcl} Y^\mu & \equiv & x^\mu\\ Y^{4+\mu} & \equiv & \mathrm{d}x^\mu/\mathrm{d}\sigma\\ \end{array} $$

Con esa identificación la ecuación de la geodésica se puede escribir como:

def ecuacion_geodesica(Y,s,Gfun,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=array(Gfun(*x))
    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

Definamos las condiciones iniciales y de integración:

N=2
from numpy import array
Y0s=array([1,0,0.1,1])
from numpy import pi,linspace
ss=linspace(0,10,100)

Y podemos integrar:

from scipy.integrate import odeint
Gfun=Gfun_cil2d()
Ys=odeint(ecuacion_geodesica,Y0s,ss,args=(Gfun,N))

Hagamos un gráfico de la trayectoria en el espacio:

import matplotlib.pyplot as plt
fig=plt.figure(figsize=(5,5))
ax=fig.gca()

from numpy import sin,cos
for i,Y in enumerate(Ys):
    
    #Coordenadas
    r=Y[0]
    teta=Y[1]
    
    #Puntos
    #Vectores unitarios
    er=array([cos(teta),sin(teta)])
    et=array([-sin(teta),cos(teta)])
    
    #Posición en la trayectoria
    rpos=r*er
    
    #Grafica puntos
    ax.plot(rpos[0],rpos[1],'k.')
    
#ax.set_xlim((0,1.5))
#ax.set_ylim((0,1.5))
ax.grid()
fig.tight_layout()

Figura 2.39.

Que coincide con lo que esperabamos: la trayectoria es una línea recta.

2.4.7.2. Geodésica sobre la superficie de una esfera

Un problema más interesante es calcular la geodésica sobre una esfera. Para ello necesitamos la métrica:

Y podemos usar los algoritmos introducidos antes para integrar la geodésica:

#Condiciones iniciales: arrancando en Medellín
N=2
from numpy import array,pi
Y0s=array([0*pi/180,-75*pi/180,0.5,0.3])
from numpy import pi,linspace
ss=linspace(0,20,100)

#Integra la ecuación de la geodésica
from scipy.integrate import odeint
R=1
Gfun=Gfun_esf2d(R)
Ys=odeint(ecuacion_geodesica,Y0s,ss,args=(Gfun,N))

#Extrae las coordenadas y las convierte a geográficas
from numpy import mod
lons=mod(Ys[:,1]*180/pi,360)
for i,lon in enumerate(lons):
    lons[i]=lon if lon<180 else lon-360
lats=Ys[:,0]*180/pi

#Grafica
import matplotlib.pyplot as plt
fig=plt.figure(figsize=(5,5))
ax=fig.gca()
ax.plot(lons,lats,'r.')

#Decoracion
ax.set_xlim((-180,180))
ax.set_ylim((-90,90))
ax.grid()
fig.tight_layout()

Figura 2.40.

Para verificar que si es una circunferencia máxima usemos el paquete Cartopy que permite representar puntos sobre un mapa realista de la Tierra.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

fig=plt.figure(figsize=(6,3))
ax=plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()
ax.plot(lons,lats,'r.')

fig.tight_layout()

Figura 2.41.

Indice | Previo: RelatividadGeneral.SimbolosChristoffel | Siguiente: RelatividadGeneral.InerciaYGeodesicas