Cambio de coordenadas curvilineas con Python SymPy

Coordenadas curvilineas con Python SymPy¶

Autor: Eduardo Martín Calleja

Este post nace de mi interés de continuar estudiando las propiedades del Fondo Cósmico de Microondas. El método habitual para el análisis de sus anisotropías requiere trabajar en coordenadas esféricas descomponiendo la señal en un desarrollo mediante una serie de armónicos esféricos. Por este motivo inicio ahora unos posts dedicados al tema de los armónicos esféricos. Estos constituyen la parte angular de las soluciones de la ecuación de Laplace:

$$\nabla^2 f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} + \frac{\partial^2 f}{\partial z^2} = 0$$

Cuando esta se expresa en coordenadas esféricas. De modo que el primer paso, el cual es el objeto de este post, será expresar el operador laplaciano en coordenadas esféricas. Por supuesto el resultado se puede encontrar con facilidad en internet y libros de texto, pero me ha parecido que puede ser un ejercicio interesante para realizar utilizando la librería de matemática simbólica SymPy de Python.

Este artículo está escrito íntegramente utilizando el Notebook de Ipython.

Importaciones y referencias¶

In [1]:

from __future__ import division
import sympy as sym
sym.init_printing()

# Generar un cuadro con versiones de las librerías utilizadas en este notebook
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information sympy

Out[1]:

Software	Version
Python	2.7.9 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython	3.1.0
OS	Linux 3.13.0 49 generic x86_64 with debian jessie sid
sympy	0.7.6
Mon Apr 20 23:55:11 2015 CEST

Este artículo está fuertemente inspirado en el siguiente ejemplo: Ejemplo de cambio de variables en coordenadas curvilineas con Sympy

También será de utilidad el Artículo sobre el tensor métrico en la Wikipedia

Y, como se verá más adelante, utilizaremos también la definición del operador de Laplace-Beltrami

Algo de teoría¶

Vamos a abordar la expresión del lagrangiano, primero en un disco y despues sobre una superficie esférica. Para ello necesitaremos expresar el lagrangiano en diferentes tipos de coordenadas. Comenzaremos con algo de teoría tomando como primer ejemplo la conversión de coordenadas cartesianas a coordenadas polares en el plano.

El tensor métrico

Comenzaremos por el caso del plano $\mathbb{R}^2$. Un sistema de coordenadas, posiblemente curvilineas, $(u, v)$ sobre el plano, será una aplicación $\varphi(\mathbf p) = (u,v)$ que a cada punto $\mathbf p$ del plano le asocia un par de números reales $(u, v)$, por ejemplo, sus coordenadas polares.

Para este sistema de coordenadas, podemos asociar a cada punto $\mathbf p$ del plano, unos llamados vectores tangentes coordenados:

$$\vec r_u = \frac{\partial \varphi^{-1}}{\partial u}(\varphi(\mathbf p)) $$$$\vec r_v = \frac{\partial \varphi^{-1}}{\partial v}(\varphi(\mathbf p)) $$

Ejemplo: coordenadas polares

Por ejemplo, el sistema de coordenadas polares asocia a cada punto $\mathbf p = (x, y)$ del plano el par de números reales $(\rho, \theta)$ definidos de la forma habitual:

$$\varphi(x, y) = (\rho, \theta)$$

La aplicación coordenada inversa vendrá dada por:

$$ \varphi^{-1}(\rho, \theta) = (x, y) = (\rho \cos \theta, \rho \sin \theta)$$

Que normalmente se expresa simplemente con las ecuaciones:

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

Los vectores tangentes coordenados serán pues:

$$ \vec r_\rho = \frac{\partial \varphi^{-1}}{\partial \rho} = \left( \frac{\partial x}{\partial \rho} , \frac{\partial y}{\partial \rho} \right) = (\cos \theta, \sin \theta)$$$$ \vec r_\theta = \frac{\partial \varphi^{-1}}{\partial \theta} = \left( \frac{\partial x}{\partial \theta} , \frac{\partial y}{\partial \theta} \right) = (-r \sin \theta, r \cos \theta)$$

Definición del tensor métrico en el plano

Volvamos al caso general de unas corrdenadas cualesquiera $(u, v)$. A partir de los vectores coordenados tangentes $\vec r_u, \vec r_v$ y utilizando el producto escalar de $\mathbb{R}^2$ definiremos las cantidades:

$E = \vec r_u \cdot \vec r_u \quad , \quad F = \vec r_u \cdot \vec r_v \quad , \quad G = \vec r_v \cdot \vec r_v $

Y a partir de ellas definimos el tensor métrico como la matriz:

$$ g = \left( g_{ij} \right) = \begin{pmatrix} E & F \\ F & G \end{pmatrix} $$

Por ejemplo, en el caso de las coordenadas polares, calculando los productos escalares con los vectores tangentes $ \vec r_\rho , \vec r_\theta $ obtenidos antes, resulta el tensor métrico:

$$ g_{\: \text{polares}}= \begin{pmatrix} 1 & 0 \\ 0 & \rho^2 \end{pmatrix} $$

Cambio de coordenadas

El tensor métrico es efectivamente de un tensor por la forma en que se transforma ante un cambio de coordenadas. En efecto, supongamos que definimos otro sistema de coordenadas $(u', v')$.

$$\varphi: (x, y) \longrightarrow (u, v)$$$$\psi: (x, y) \longrightarrow (u', v')$$

Podemos hacer depender las coordenadas antiguas de las nuevas mediante la transformación:

$$\varphi \circ \psi^{-1}: (u', v') \longrightarrow (u, v)$$

Cuya matriz Jacobiana será

$$ J = \begin{pmatrix} \frac{\partial u}{\partial u'} & \frac{\partial u}{\partial v'} \\ \frac{\partial v}{\partial u'} & \frac{\partial v}{\partial v'} \end{pmatrix} $$

Y aplicando las reglas de derivación puede comprobarse que el tensor métrico expresado en las nuevas coordenadas será:

$$ \begin{pmatrix} E' & F' \\ F' & G' \end{pmatrix} = J^T \begin{pmatrix} E & F \\ F & G \end{pmatrix} J $$

Obtención del tensor métrico en coordenadas polares

A continuación, y como aplicación de lo anterior, vamos a comprobar que, mediante un cambio de coordenadas a partir del tensor métrico en coordenadas cartesianas (la matriz unidad), también llegamos a obtener el tensor métrico en coordenadas polares:

In [2]:

from sympy import symbols, sin, cos, Matrix, eye, simplify, sqrt, Function, expand, sympify

In [3]:

# En SymPy siempre hay que definir los símbolos que vamos a utilizar
theta = symbols('theta', real=True)
rho = symbols('rho', real=True, negative=False)

In [4]:

#Ecuaciones del cambio de coordenadas
x = rho * cos(theta)
y = rho * sin(theta)

In [5]:

C = [x, y]    # Coordenadas cartesianas
P = [rho, theta] # coordenadas polares
C

Out[5]:

$$\left [ \rho \cos{\left (\theta \right )}, \quad \rho \sin{\left (\theta \right )}\right ]$$

In [6]:

# Cálculo del Jacobiano de las coordenadas C respecto de las P
J = Matrix(C).jacobian(P) # Matrix se utiliza para convertir C
                        # en un objeto sympy que posea el método Jacobian
J

Out[6]:

$$\left[\begin{matrix}\cos{\left (\theta \right )} & - \rho \sin{\left (\theta \right )}\\\sin{\left (\theta \right )} & \rho \cos{\left (\theta \right )}\end{matrix}\right]$$

In [7]:

eye(2)   # este es el tensor métrico en coordenadas cartesianas

Out[7]:

$$\left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]$$

In [8]:

# Y este es el tensor métrico en coordenadas polares
g = simplify(J.T * eye(2) * J)
g

Out[8]:

$$\left[\begin{matrix}1 & 0\\0 & \rho^{2}\end{matrix}\right]$$

Elemento de longitud

A partir de la métrica se puede obtener el elemento de longitud $ds$, Su cuadrado se define mediante la expresión:

$$ ds^2 = (du, dv) \begin{pmatrix} E & F \\ F & G \end{pmatrix} \begin{pmatrix} du \\ dv \end{pmatrix} $$

Y no depende del sistema de coordenadas elegido. Por esta razón, a $ds^2$ se le llama primera forma fundamental.

In [9]:

# Ejemplo: elemento de longitud en coordenadas polares:
drho, dtheta = symbols(r'd\rho d\theta')
dP = Matrix([drho, dtheta])
ds2 = dP.T * g * dP
ds2[0]

Out[9]:

$$d\rho^{2} + d\theta^{2} \rho^{2}$$

Elemento de área

El elemento de área $da$ es otro elemento intrínseco de la superficie, que no depende del sistema de coordenadas elegido. Se obtiene como:

$$da = \sqrt{\det(g)} du dv$$

In [10]:

# Ejemplo: elemento de área en coordenadas polares:
da = sqrt(g.det())*drho*dtheta
da

Out[10]:

$$d\rho d\theta \rho$$

El operador de Laplace-Beltrami

El operador de Laplace-Beltrami es un operador diferencial que generaliza el operador de Laplace para poder ser utilizado en un sistema $(\theta_1, \theta_2, \dots \:)$ de coordenadas curvilineas. Su expresión es:

$$ \nabla^2 f = \frac{1}{\sqrt{\lvert g \rvert}} \frac{\partial}{\partial \theta_i} \left( \sqrt{\lvert g \rvert} g^{ij} \frac{\partial f}{\partial \theta_j} \right) $$

Donde:

• $g$ es el tensor métrico, $\vert g \vert$ el determinante de su matriz.
• $g^{ij}$ son las componentes de la matriz inversa de la matriz del tensor métrico.

Así pues, vamos a utilizar esa expresión para obtener el laplaciano en un sistema de coordenadas diferentes de las cartesianas.

El laplaciano en coordenadas polares¶

En primer lugar vamos a utilizar el operador de Laplace_Beltrami para expresar el operador laplaciano en coordenadas polares $(\rho, \theta)$

In [11]:

f = sym.Function('f')
f = f(*list(P))    # Expresamos f como función de las coordenadas polares
f

Out[11]:

$$f{\left (\rho,\theta \right )}$$

In [12]:

#Calculemos ahora la inversa de la métrica:
g_inv = g.inv()
g_inv

Out[12]:

$$\left[\begin{matrix}1 & 0\\0 & \frac{1}{\rho^{2}}\end{matrix}\right]$$

In [13]:

# Y el determinante de la métrica:
g_det = g.det()
g_det

Out[13]:

$$\rho^{2}$$

In [14]:

# Y finalmente podemos expresar el operador de Laplace-Beltrami siguiendo su definición:
L_B_pol = 0
for i in range(2):
    for j in range(2):
        L_B_pol += (sqrt(g_det) * g_inv[i,j] * f.diff(P[j])).diff(P[i])
L_B_pol = expand(L_B_pol / sqrt(g_det))
L_B_pol

Out[14]:

$$\frac{\partial^{2}}{\partial \rho^{2}} f{\left (\rho,\theta \right )} + \frac{1}{\rho} \frac{\partial}{\partial \rho} f{\left (\rho,\theta \right )} + \frac{1}{\rho^{2}} \frac{\partial^{2}}{\partial \theta^{2}} f{\left (\rho,\theta \right )}$$

En caso de que quisieramos utilizar esta expresión en otro notebook, se puede generar a partir de ella un string, el cual se podría copiar al portapapeles y generar la expresión de nuevo con sympify()

In [15]:

str(L_B_pol)

Out[15]:

'Derivative(f(rho, theta), rho, rho) + Derivative(f(rho, theta), rho)/rho + Derivative(f(rho, theta), theta, theta)/rho**2'

In [16]:

sympify('Derivative(f(rho, theta), rho, rho) + \
        Derivative(f(rho, theta), rho)/rho +   \
        Derivative(f(rho, theta), theta, theta)/rho**2')

Out[16]:

$$\frac{\partial^{2}}{\partial \rho^{2}} f{\left (\rho,\theta \right )} + \frac{1}{\rho} \frac{\partial}{\partial \rho} f{\left (\rho,\theta \right )} + \frac{1}{\rho^{2}} \frac{\partial^{2}}{\partial \theta^{2}} f{\left (\rho,\theta \right )}$$

El laplaciano en coordenadas esféricas¶

In [17]:

# theta es la colatitud, phi la longitud
rho, theta, phi = sym.symbols('rho theta phi', real=True, negative=False)
rho, theta, phi

Out[17]:

$$\left ( \rho, \quad \theta, \quad \phi\right )$$

In [18]:

# Ecuaciones de cambio de coordenadas
x = rho * sin(theta) * cos(phi)
y = rho * sin(theta) * sin(phi)
z = rho * cos(theta)

In [19]:

C = [x, y, z]    # Coordenadas cartesianas
E = [rho, theta, phi] # Coordenadas esféricas
C

Out[19]:

$$\left [ \rho \sin{\left (\theta \right )} \cos{\left (\phi \right )}, \quad \rho \sin{\left (\phi \right )} \sin{\left (\theta \right )}, \quad \rho \cos{\left (\theta \right )}\right ]$$

In [20]:

# Cálculo del Jacobiano de las coordenadas C respecto de las E
J = Matrix(C).jacobian(E)
J

Out[20]:

$$\left[\begin{matrix}\sin{\left (\theta \right )} \cos{\left (\phi \right )} & \rho \cos{\left (\phi \right )} \cos{\left (\theta \right )} & - \rho \sin{\left (\phi \right )} \sin{\left (\theta \right )}\\\sin{\left (\phi \right )} \sin{\left (\theta \right )} & \rho \sin{\left (\phi \right )} \cos{\left (\theta \right )} & \rho \sin{\left (\theta \right )} \cos{\left (\phi \right )}\\\cos{\left (\theta \right )} & - \rho \sin{\left (\theta \right )} & 0\end{matrix}\right]$$

In [21]:

# Y este es el tensor métrico en coordenadas esféricas
g = simplify(J.T * eye(3) * J)
g

Out[21]:

$$\left[\begin{matrix}1 & 0 & 0\\0 & \rho^{2} & 0\\0 & 0 & \rho^{2} \sin^{2}{\left (\theta \right )}\end{matrix}\right]$$

In [22]:

# Y el determinante del tensor métrico:
g_det = g.det()
g_det

Out[22]:

$$\rho^{4} \sin^{2}{\left (\theta \right )}$$

In [23]:

#Calculemos ahora la inversa de la métrica:
g_inv = g.inv()
g_inv

Out[23]:

$$\left[\begin{matrix}1 & 0 & 0\\0 & \frac{1}{\rho^{2}} & 0\\0 & 0 & \frac{1}{\rho^{2} \sin^{2}{\left (\theta \right )}}\end{matrix}\right]$$

In [24]:

f = sym.Function('f')
f = f(*list(E))    # Expresamos f como función de las coordenadas esféricas
f

Out[24]:

$$f{\left (\rho,\theta,\phi \right )}$$

In [25]:

# Y finalmente podemos obtener la expresión del operador de Laplace-Beltrami siguiendo su definición:
L_B_esf = 0
for i in range(3):
    for j in range(3):
        L_B_esf += (sqrt(g_det) * g_inv[i,j] * f.diff(E[j])).diff(E[i])
L_B_esf = expand(L_B_esf / sqrt(g_det))
L_B_esf

Out[25]:

$$\frac{\partial^{2}}{\partial \rho^{2}} f{\left (\rho,\theta,\phi \right )} + \frac{2}{\rho} \frac{\partial}{\partial \rho} f{\left (\rho,\theta,\phi \right )} + \frac{\cos{\left (\theta \right )} \frac{\partial}{\partial \theta} f{\left (\rho,\theta,\phi \right )}}{\rho^{2} \left\lvert{\sin{\left (\theta \right )}}\right\rvert} \operatorname{sign}{\left (\sin{\left (\theta \right )} \right )} + \frac{1}{\rho^{2}} \frac{\partial^{2}}{\partial \theta^{2}} f{\left (\rho,\theta,\phi \right )} + \frac{\frac{\partial^{2}}{\partial \phi^{2}} f{\left (\rho,\theta,\phi \right )}}{\rho^{2} \sin^{2}{\left (\theta \right )}}$$

Es el resultado correcto, siempre que observemos que el signo y el valor absoluto en sen($\theta$) se cancelan (sympy al parecer no lo hace), quedando simplemente $\sin(\theta)$ en el denominador.

No obstante, se puede llegar a esta expresión mediante fuerza bruta, eliminando el valor absoluto, y sustituyendo el signo por 1:

In [26]:

L_B_esf = L_B_esf.replace(sym.Abs, sym.Id).replace(sym.sign(sin(theta)), 1)
L_B_esf

Out[26]:

$$\frac{\partial^{2}}{\partial \rho^{2}} f{\left (\rho,\theta,\phi \right )} + \frac{2}{\rho} \frac{\partial}{\partial \rho} f{\left (\rho,\theta,\phi \right )} + \frac{1}{\rho^{2}} \frac{\partial^{2}}{\partial \theta^{2}} f{\left (\rho,\theta,\phi \right )} + \frac{\cos{\left (\theta \right )} \frac{\partial}{\partial \theta} f{\left (\rho,\theta,\phi \right )}}{\rho^{2} \sin{\left (\theta \right )}} + \frac{\frac{\partial^{2}}{\partial \phi^{2}} f{\left (\rho,\theta,\phi \right )}}{\rho^{2} \sin^{2}{\left (\theta \right )}}$$

El laplaciano esférico¶

Es el operador laplaciano (mejor, el operador de Laplace-Beltrami) restringido a la superficie esférica de radio $\rho = 1$. En este caso estamos ante un operador diferencial sobre una superficie cuyas coordenadas locales son $(\theta, \phi)$. Apliquemos el mismo método que antes para obtener su expresión en función de estas coordenadas:

In [27]:

# Ecuaciones de cambio de coordenadas
x = sin(theta) * cos(phi)
y = sin(theta) * sin(phi)
z = cos(theta)

C = [x, y, z]    # Coordenadas cartesianas
E2 = [theta, phi] # Coordenadas esféricas

J = Matrix(C).jacobian(E2)
g = simplify(J.T * eye(3) * J)
g_det = g.det()
g_inv = g.inv()

In [28]:

f = sym.Function('f')
f = f(*list(E2))    # Expresamos f como función de las coordenadas esféricas
f

Out[28]:

$$f{\left (\theta,\phi \right )}$$

In [29]:

# Y finalmente obtenemos la expresión del operador de Laplace-Beltrami sobre la esfera de radio 1
L_B_s_esf = 0
for i in range(2):
    for j in range(2):
        L_B_s_esf += (sqrt(g_det) * g_inv[i,j] * f.diff(E2[j])).diff(E2[i])
L_B_s_esf = expand(L_B_s_esf / sqrt(g_det))
L_B_s_esf = L_B_s_esf.replace(sym.Abs, sym.Id).replace(sym.sign(sin(theta)), 1)
L_B_s_esf

Out[29]:

$$\frac{\partial^{2}}{\partial \theta^{2}} f{\left (\theta,\phi \right )} + \frac{\frac{\partial}{\partial \theta} f{\left (\theta,\phi \right )}}{\sin{\left (\theta \right )}} \cos{\left (\theta \right )} + \frac{\frac{\partial^{2}}{\partial \phi^{2}} f{\left (\theta,\phi \right )}}{\sin^{2}{\left (\theta \right )}}$$