Curvilinear coordinates with Python SymPy¶
Author: Eduardo Martín Calleja
This post comes from my interest to continue studying the properties of the cosmic microwave background. The usual method for the analysis of their anisotropies requires working in spherical coordinates decomposing the signal into a series of spherical harmonics. Therefore I'm now beginning some posts devoted to the subject of spherical harmonics. These come as the angular component of the solution to Laplace's equation:
$$\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$$when this is expressed in spherical coordinates. So the first step, which is the subject of this post, is to write the Laplacian operator in spherical coordinates. Of course the result can be found easily on the internet and textbooks, but I thought it might be interesting to do it using the SymPy symbolic math library for Python as an exercise.
This article has been written entirely using the Notebook of IPython.
Imports and references¶
from __future__ import division
import sympy as sym
sym.init_printing()
# This IPython magic generates a table with version information
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information sympy
This article is heavily inspired by the following example: Example of changing variables in curvilinear coordinates with SymPy
Also useful the metric tensor article on Wikipedia
And, as will be seen below, we'll also need the definition of the Laplace-Beltrami operator
Some theory¶
Let's address the expression of the Lagrangian in coordinates different of the cartesian ones. We'll begin with some theory taking as a first example the conversion of cartesian coordinates to polar coordinates in the plane.
The metric tensor
Let's begin with the case of the plane $\mathbb{R}^2$.A coordinate system, possibly curvilinear, $(u, v)$ on the plane, is an application $\varphi(\mathbf p) = (u,v)$ which associates to each point $\mathbf p$ of the plane a pair of real numbers $(u, v)$, for example, its polar coordinates.
Given a coordinate system, we can associate to each point $\mathbf p$ in the plane, two vectors called coordinate tangent vectors
$$\vec r_u = \frac{\partial \varphi^{-1}}{\partial u}(\varphi(\mathbf p)) $$$$\vec r_v = \frac{\partial \varphi^{-1}}{\partial v}(\varphi(\mathbf p)) $$Example: polar coordinates
For example, the polar coordinate system associate to each point $\mathbf p = (x, y)$ in the plane two real numbers $(\rho, \theta)$ defined in the usual manner.
$$\varphi(x, y) = (\rho, \theta)$$The inverse coordinates application will be given by:
$$ \varphi^{-1}(\rho, \theta) = (x, y) = (\rho \cos \theta, \rho \sin \theta)$$That is simply expressed by the equations:
$$x = \rho \cos \theta $$$$y = \rho \sin \theta $$The coordinate tangent vectors will therefore be:
$$ \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)$$Definition of the metric tensor in the plane
Returning to the general case of any $(u, v)$ coordinates, taking the coordinate tangent vectors and using the dot product of $\mathbb{R}^2$, we can define the following quantities:
$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 $
And from them we can define the metric tensor as the matrix:
$$ g = \left( g_{ij} \right) = \begin{pmatrix} E & F \\ F & G \end{pmatrix} $$For example, in the case of polar coordinates, calculating the scalar products with the tangent vectors $ \vec r_\rho , \vec r_\theta $ obtained above, the metric tensor is:
$$ g_{\: \text{polar}}= \begin{pmatrix} 1 & 0 \\ 0 & \rho^2 \end{pmatrix} $$Change of coordinates
The metric tensor is indeed a tensor by the way it is transformed in a change of coordinates. Indeed, suppose that we have another coordinate system $(u', v')$.
$$\varphi: (x, y) \longrightarrow (u, v)$$$$\psi: (x, y) \longrightarrow (u', v')$$We can make the old coordinates dependent on the new ones by the transformation:
$$\varphi \circ \psi^{-1}: (u', v') \longrightarrow (u, v)$$Whose Jacobian matrix is
$$ 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} $$And applying the rules of derivation it can be seen that the metric tensor expressed in the new coordinates is:
$$ \begin{pmatrix} E' & F' \\ F' & G' \end{pmatrix} = J^T \begin{pmatrix} E & F \\ F & G \end{pmatrix} J $$Getting the metric tensor in polar coordinates
Now, as an application of the above, we will check that, by a change of coordinates from the metric tensor in cartesian coordinates (the unit matrix), we can get the metric tensor in polar coordinates:
from sympy import symbols, sin, cos, Matrix, eye, simplify, sqrt, Function, expand, sympify
# In SymPy we must always define the symbols that we will use
theta = symbols('theta', real=True)
rho = symbols('rho', real=True, negative=False)
#coordinate change equations
x = rho * cos(theta)
y = rho * sin(theta)
C = [x, y] # cartesian coordinates
P = [rho, theta] # polar coordinates
# Calculate the Jacobian of C coordinates with respect to P
J = Matrix(C).jacobian(P) # Matrix() is used to convert C
# into a SymPy object with the Jacobian method
J
eye(2) # this is the metric tensor in cartesian coordinates
# And this is the metric tensor in polar coordinates
g = simplify(J.T * eye(2) * J)
g
Line element
From the metric it can be obtained the element of length $ds$, it's square is defined by the expression:
$$ ds^2 = (du, dv) \begin{pmatrix} E & F \\ F & G \end{pmatrix} \begin{pmatrix} du \\ dv \end{pmatrix} $$This doesn't depend on the chosen coordinate system. For this reason, $ds^2$ is called first fundamental form.
# Example: line element in polar coordinates:
drho, dtheta = symbols(r'd\rho d\theta')
dP = Matrix([drho, dtheta])
ds2 = dP.T * g * dP
ds2[0]
Differential element of area
The element of area gives another intrinsic element of the surface, which does not depend on the chosen coordinate system. We can express it as:
$$da = \sqrt{\det(g)} \; du \; dv$$#Example: element of area in polar coordinates:
da = sqrt(g.det())*drho*dtheta
da
The Laplace-Beltrami operator
The Laplace-Beltrami operator is a differential operator that generalizes the Laplace operator to be used with a system $(\theta_1, \theta_2, \dots \:)$ of curvilinear coordinates. Its expression is:
$$ \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) $$Where:
- $g$ is the metric tensor, $\vert g \vert$ the determinant of its matrix.
- $g^{ij}$ are the components of the inverse of the matrix of the metric tensor.
So, we will use this expression to obtain the Laplacian on a system of coordinates different of the cartesian ones.
The Laplacian in polar coordinates¶
First we will use the Laplace_Beltrami operator to express the Laplacian operator in polar coordinates $(\rho, \theta)$
f = sym.Function('f')
f = f(*list(P)) # Express f as a function of the polar coordinates
#Computing the inverse of the metric:
g_inv = g.inv()
g_inv
# And the determinant of the metric:
g_det = g.det()
g_det
# And finally we can express the Laplace-Beltrami operator following its definition:
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
If you would like to use this expression in another notebook, you can generate a string from it, which could be copied to the clipboard and generate the expression again with sympify()
str(L_B_pol)
sympify('Derivative(f(rho, theta), rho, rho) + \
Derivative(f(rho, theta), rho)/rho + \
Derivative(f(rho, theta), theta, theta)/rho**2')
The Laplacian in spherical coordinates¶
# Theta is the co-latitude, phi the longitude
rho, theta, phi = sym.symbols('rho theta phi', real=True, negative=False)
rho, theta, phi
# Change of coordinates equations:
x = rho * sin(theta) * cos(phi)
y = rho * sin(theta) * sin(phi)
z = rho * cos(theta)
C = [x, y, z] # cartesian coordinates
S = [rho, theta, phi] # spherical coordinates
C
# Calculate the Jacobian of C coordinates with respect to S
J = Matrix(C).jacobian(S)
J
# This is the metric tensor in spherical coordinates
g = simplify(J.T * eye(3) * J)
g
# And the determinant of the metric tensor:
g_det = g.det()
g_det
#The inverse of the metric:
g_inv = g.inv()
g_inv
f = sym.Function('f')
f = f(*list(S)) # express f as a function of the spherical coordinates
f
# And finally we get the expression of the Laplace-Beltrami operator:
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(S[j])).diff(S[i])
L_B_esf = expand(L_B_esf / sqrt(g_det))
L_B_esf
This is the correct result, provided that we notice that the sign and absolute value in sin($\theta$) cancel out (SymPy apparently does not), just letting $\sin(\theta)$ in the denominator.
However, you can simplify this expression by brute force, eliminating the absolute value sign, and replacing the sign by 1:
L_B_esf = L_B_esf.replace(sym.Abs, sym.Id).replace(sym.sign(sin(theta)), 1)
L_B_esf
The spherical Laplacian¶
It is the Laplacian operator (better, the Laplace-Beltrami operator) restricted to the spherical surface of radius $\rho = 1$. In this case we have a differential operator on a surface whose local coordinates are the angular ones: $(\theta, \phi)$. Let us apply the same method as above to obtain its expression in terms of these coordinates:
# Equations of change of coordinates
x = sin(theta) * cos(phi)
y = sin(theta) * sin(phi)
z = cos(theta)
C = [x, y, z] # cartesian coordinates
S2 = [theta, phi] # Spherical coordinates
J = Matrix(C).jacobian(S2)
g = simplify(J.T * eye(3) * J)
g_det = g.det()
g_inv = g.inv()
f = sym.Function('f')
f = f(*list(S2)) # express f as a function of the spherical coordinates
f
# Finally we obtain the expression of the Laplace-Beltrami operator on the unit sphere
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(S2[j])).diff(S2[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
