Python (SymPy) Programming for finding Laplacian in Polar Coordinates

This article introduces a way to find the Laplace operator (Laplacian) in polar coordinates. The Laplace operator of 2-dimensional and 3-dimensional spaces often appears in problems in electromagnetism, quantum mechanics and so on, for those with rotational symmetry.

Though the Laplacian in the Cartesian coordinates takes such a simple form, once one moves to other coordinates, e.g. polar coordinates, it converts into such a complex form. The process of the conversion from Cartesian to polar coordinates asks your patience with dealing with many terms made of trigonometric functions. Simplifying those expressions is such time-consuming and less productive. There is a way to avoid it: USE your computer!

SymPy is such a Python library that makes it possible. The following article introduces basics for some usage of SymPy:

pianofisica.hatenablog.com

This article demonstrates how to obtain expressions for the Laplacian in polar coordinates in 2-dimensional and 3-dimensional spaces using SymPy. In addition to them, 4-dim. case is examined and resultants for 5- and 6-dim. are provided without showing the source code. (It is a great advantage of using computers that one can carry out extensions in such an easy fashion.)

Laplacian in 2-dimensional space
Laplacian in 3-dimensional space
Laplacian in 4-dimensional space
Laplacian in 5-dimensional space
Laplacian in 6-dimensional space

Laplacian in 2-dimensional space

Let it start with the 2-dimensional Laplacian:

$\qquad\displaystyle{\triangle^{(2)}=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}}$

The relation between 2-dim. polar coordinates $(r,\theta)$ and the Cartesian coordintes $(x,y)$ is given by

$\qquad\displaystyle{x=r\cos\theta\\ y=r\sin\theta}$

equivalently

$\qquad\displaystyle{r=\sqrt{x^2+y^2}\\ \theta=\arctan\left(\frac{y}{x}\right)}$

The Leibniz' chain rule gives

$\qquad\displaystyle{\left(\begin{array}{c} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \end{array} \right)=\left(\begin{array}{c} \frac{\partial r}{\partial x}\frac{\partial }{\partial r}+\frac{\partial \theta}{\partial x}\frac{\partial }{\partial \theta} \\ \frac{\partial r}{\partial y}\frac{\partial }{\partial r}+\frac{\partial \theta}{\partial y}\frac{\partial }{\partial \theta} \end{array}\right)=R_2\left(\begin{array}{c}\frac{\partial }{\partial r}\\ \frac{\partial }{\partial \theta} \end{array}\right)}$

where the matrix $R_2$ is introduced by

$\qquad\displaystyle{R_2=\left(\begin{array}{cc} \frac{\partial r}{\partial x} & \frac{\partial \theta}{\partial x} \\ \frac{\partial r}{\partial y} & \frac{\partial \theta}{\partial y} \end{array}\right)}$

The expression for the matrix $R_2$ is obtained by running the following code in Python:

import sympy as sp

sp.var('x, y')
r = sp.sqrt(x**2+y**2)
theta = sp.atan(y/x)

r2 = sp.simplify(sp.Matrix([
[sp.diff(r,x),sp.diff(theta,x)],
[sp.diff(r,y),sp.diff(theta,y)]
]))

del(r, theta)
sp.var('r, theta', positive = True)

R2 = sp.trigsimp(r2.subs(x, r*sp.cos(theta)).subs(y, r*sp.sin(theta)))

R2

resulting in

$\qquad\displaystyle{R_2=\left(\begin{array}{cc} \cos\theta & -\frac{\sin \theta}{r} \\ \sin\theta & \frac{\cos\theta}{r} \end{array}\right)}$

Introducing an arbitrary test function $f=f(r,\theta)$ , and a matrix $T_2$ via

$\qquad\displaystyle{T_2=\left(\begin{array}{c} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \end{array} \right)=R_2\left(\begin{array}{c} \frac{\partial f}{\partial r} \\ \frac{\partial f}{\partial \theta} \end{array} \right)}$

then one can find the Laplacian as follows:

$\qquad\displaystyle{\begin{aligned}\triangle^{(2)} f&=\left(\begin{array}{cc} \frac{\partial}{\partial x} & \frac{\partial}{\partial y} \end{array} \right)\left(\begin{array}{c} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \end{array} \right)\\&=\sum_{i=1}^2\left((R_2)_{i1}\frac{\partial}{\partial r}+(R_2)_{i2}\frac{\partial}{\partial \theta}\right)(T_2)_i\end{aligned}}$

Running the following code in Python

from sympy.core.function import Function

f = Function('f')(r,theta)

T2 = R2*sp.Matrix([[sp.diff(f,r)],[sp.diff(f,theta)]])

i = sp.symbols('i',integer=True)

sp.trigsimp(sp.summation(
R2[i,0]*sp.diff(T2[i,0],r)+R2[i,1]*sp.diff(T2[i,0],theta),
(i,0,1)))

one may obtain

$\qquad\displaystyle{\triangle^{(2)}=\frac{\partial^2}{\partial r^2}+\frac{1}{r}\frac{\partial}{\partial r} +\frac{1}{r^2}\frac{\partial^2}{\partial \theta^2}}$

Laplacian in 3-dimensional space

Let it move to 3-dim. case:

$\qquad\displaystyle{\triangle^{(3)}=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}}$

The relation between the polar coordinates $(r,\theta,\varphi)$ and the Cartesian coordinates $(x,y,z)$ is

$\qquad\displaystyle{x=r\sin\theta\cos\varphi\\ y=r\sin\theta\sin\varphi \\ z=r\cos\theta}$

equivalently

$\qquad\displaystyle{r=\sqrt{x^2+y^2+z^2}\\ \theta=\arctan\left(\frac{\sqrt{x^2+y^2}}{z}\right)\\ \varphi=\arctan\left(\frac{y}{x}\right)}$

Let $f=f(r,\theta,\varphi)$ be an arbitrary test function and $R_3$ be

$\qquad\displaystyle{R_3=\left(\begin{array}{ccc} \frac{\partial r}{\partial x} & \frac{\partial \theta}{\partial x} & \frac{\partial \varphi}{\partial x} \\ \frac{\partial r}{\partial y} & \frac{\partial \theta}{\partial y} & \frac{\partial \varphi}{\partial y} \\ \frac{\partial r}{\partial z} & \frac{\partial \theta}{\partial z} & \frac{\partial \varphi}{\partial z} \end{array}\right) \\ T_3=R_3\left(\begin{array}{c} \frac{\partial f}{\partial r} \\ \frac{\partial f}{\partial \theta} \\ \frac{\partial f}{\partial \varphi} \end{array} \right)}$

then one can find the Laplacian via

$\qquad\displaystyle{\triangle^{(3)} f=\sum_{i=1}^3\left((R_3)_{i1}\frac{\partial}{\partial r}+(R_3)_{i2}\frac{\partial}{\partial \theta}+(R_3)_{i3}\frac{\partial}{\partial \varphi}\right)(T_3)_i}$

First, let all the symbols initialized in the code:

%reset

then run the following code

import sympy as sp

sp.var('x, y, z')
r = sp.sqrt(x**2+y**2+z**2)
theta = sp.atan(sp.sqrt(x**2+y**2)/z)
phi = sp.atan(y/x)

r3 = sp.simplify(sp.Matrix([
[sp.diff(r,x),sp.diff(theta,x),sp.diff(phi,x)],
[sp.diff(r,y),sp.diff(theta,y),sp.diff(phi,y)],
[sp.diff(r,z),sp.diff(theta,z),sp.diff(phi,z)]
]))

del(r, theta, phi)
sp.var('r, theta, phi', nonnegative = True)

R3 = sp.expand(sp.trigsimp(
r3.subs(x, r*sp.sin(theta)*sp.cos(phi)).
subs(y, r*sp.sin(theta)*sp.sin(phi)).subs(z, r*sp.cos(theta))
), trig = True).subs(sp.Abs(sp.sin(theta)), sp.sin(theta))

T3 = R3*sp.Matrix([
[sp.diff(r,x),sp.diff(theta,x),sp.diff(phi,x)],
[sp.diff(r,y),sp.diff(theta,y),sp.diff(phi,y)],
[sp.diff(r,z),sp.diff(theta,z),sp.diff(phi,z)]
])

from sympy.core.function import Function

f = Function('f')(r,theta,phi)

T3 = R3*sp.Matrix([[sp.diff(f,r)],[sp.diff(f,theta)],[sp.diff(f,phi)]])

i = sp.symbols('i',integer=True)

sp.trigsimp(sp.summation(
R3[i,0]*sp.diff(T3[i,0],r)+R3[i,1]*sp.diff(T3[i,0],theta)+R3[i,2]*sp.diff(T3[i,0],phi),
(i,0,2)))

then one may find

$\qquad\displaystyle{\triangle^{(3)}=\frac{\partial^2}{\partial r^2}+\frac{2}{r}\frac{\partial}{\partial r} +\frac{1}{r^2}\frac{\partial^2}{\partial \theta^2}+\frac{\cos\theta}{r^2\sin\theta}\frac{\partial}{\partial \theta} +\frac{1}{r^2\sin^2\theta}\frac{\partial^2}{\partial \varphi^2}}$

Laplacian in 4-dimensional space

The 4-dim. Laplacian

$\qquad\displaystyle{\triangle^{(4)}=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}+\frac{\partial^2}{\partial w^2}}$

is converted on the polar coordinates $(r,\theta,\varphi,\psi)$ whose relation to $(x,y,z,w)$ is given by

$\qquad\displaystyle{x=r\sin\theta\sin\varphi\sin\psi \\ y=r\sin\theta\sin\varphi\cos\psi\\ z=r\sin\theta\cos\varphi\\ w=r\cos\theta}$

$\qquad\displaystyle{r=\sqrt{x^2+y^2+z^2+w^2}\\ \theta=\arctan\left(\frac{\sqrt{x^2+y^2+z^2}}{w}\right)\\ \varphi=\arctan\left(\frac{\sqrt{x^2+y^2}}{z}\right)\\ \psi=\arctan\left(\frac{\sqrt{x^2}}{y}\right)}$

One may be able to figure out the pattern. Without explaining the details any more, but running the following code

%reset

import sympy as sp

sp.var('x, y, z, w')
r = sp.sqrt(x**2+y**2+z**2+w**2)
theta = sp.atan(sp.sqrt(x**2+y**2+z**2)/w)
phi = sp.atan(sp.sqrt(x**2+y**2)/z)
psi =  sp.atan(sp.sqrt(x**2)/y)

r4 = sp.simplify(sp.Matrix([
[sp.diff(r,x),sp.diff(theta,x),sp.diff(phi,x),sp.diff(psi,x)],
[sp.diff(r,y),sp.diff(theta,y),sp.diff(phi,y),sp.diff(psi,y)],
[sp.diff(r,z),sp.diff(theta,z),sp.diff(phi,z),sp.diff(psi,z)],
[sp.diff(r,w),sp.diff(theta,w),sp.diff(phi,w),sp.diff(psi,w)]
]))

del(r, theta, phi, psi)
sp.var('r, theta, phi, psi', nonnegative = True)

R4 = sp.expand(sp.trigsimp(
r4.subs(x, r*sp.sin(theta)*sp.sin(phi)*sp.sin(psi)).
subs(y, r*sp.sin(theta)*sp.sin(phi)*sp.cos(psi)).
subs(z, r*sp.sin(theta)*sp.cos(phi)).
subs(w, r*sp.cos(theta))).
subs(sp.Abs(sp.sin(theta)), sp.sin(theta)).
subs(sp.Abs(sp.sin(phi)), sp.sin(phi)).
subs(sp.Abs(sp.sin(psi)), sp.sin(psi)), trig = True)

T4 = R4*sp.Matrix([
[sp.diff(r,x),sp.diff(theta,x),sp.diff(phi,x),sp.diff(psi,x)],
[sp.diff(r,y),sp.diff(theta,y),sp.diff(phi,y),sp.diff(psi,y)],
[sp.diff(r,z),sp.diff(theta,z),sp.diff(phi,z),sp.diff(psi,z)],
[sp.diff(r,w),sp.diff(theta,w),sp.diff(phi,w),sp.diff(psi,w)]
])

from sympy.core.function import Function

f = Function('f')(r,theta,phi,psi)

T4 = R4*sp.Matrix([[sp.diff(f,r)],[sp.diff(f,theta)],[sp.diff(f,phi)],[sp.diff(f,psi)]])

i = sp.symbols('i',integer=True)

sp.trigsimp(sp.summation(
R4[i,0]*sp.diff(T4[i,0],r)+R4[i,1]*sp.diff(T4[i,0],theta)
+R4[i,2]*sp.diff(T4[i,0],phi)+R4[i,3]*sp.diff(T4[i,0],psi),
(i,0,3)))

gives

$\displaystyle{\triangle^{(4)}={{\partial ^2}\over{\partial r^2}}+{{3}\over{r}}{{\partial }\over{\partial r}}+{{1}\over{r^2}}{{\partial^2}\over{\partial \theta^2}}+{{2\cos \theta}\over{r^2\sin \theta}}{{\partial }\over{\partial \theta}}\\\qquad\qquad\qquad+{{1}\over{r^2\sin ^2\theta}}{{\partial^2}\over{\partial \varphi^2}}+{{\cos \varphi}\over{r^2\sin ^2\theta\sin \varphi}}{{\partial}\over{\partial \varphi}}+{{1}\over{r^2\sin ^2\theta\sin ^2\varphi}}{{\partial^2}\over{\partial \psi^2}} }$

It will be a good exercise to make your hypothesis for 5-dim, 6-dim, etc. and check them via Python.

Here are the answers:

Laplacian in 5-dimensional space

The polar coord. in 5-dim.

$\qquad\displaystyle{x=r\sin\theta\sin\varphi\sin\psi\sin\chi\\y=r\sin\theta\sin\varphi\sin\psi\cos\chi\\z=r\sin\theta\sin\varphi\cos\psi\\w=r\sin\theta\cos\varphi\\s=r\cos\theta}$

converts the Laplacian into

$\displaystyle{\triangle^{(5)}=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}+\frac{\partial^2}{\partial w^2}+\frac{\partial^2}{\partial s^2}\\\quad\ \ \ \, ={{\partial^2}\over{\partial r^2}}+{{4}\over{r}}{{\partial}\over{\partial r}}+{{1}\over{r^2}}{{\partial^2}\over{\partial\theta^2}}+{{3\cos \theta}\over{r^2\sin \theta}}{{\partial}\over{\partial\theta}}+{{1}\over{r^2\sin ^2\theta}}{{\partial^2}\over{\partial\varphi^2}}+{{2\cos \varphi}\over{r^2\sin^2\theta\sin\varphi}}{{\partial}\over{\partial\varphi}}\\\qquad\quad+{{1}\over{r^2\sin^2\theta\sin^2\varphi}}{{\partial^2}\over{\partial\psi^2}}+{{\cos\psi}\over{r^2\sin^2\theta\sin^2\varphi\sin\psi}}{{\partial}\over{\partial\psi}}+{{1}\over{r^2\sin ^2\theta\sin ^2\varphi\sin ^2\psi}}{{\partial^2}\over{\partial\chi^2}}}$

Laplacian in 6-dimensional space

The polar coord. in 6-dim.

$\qquad\displaystyle{x=r\sin\theta\sin\varphi\sin\psi\sin\chi\sin\eta\\y=r\sin\theta\sin\varphi\sin\psi\sin\chi\cos\eta\\z=r\sin\theta\sin\varphi\sin\psi\cos\chi\\w=r\sin\theta\sin\varphi\cos\psi\\s=r\sin\theta\cos\varphi\\t=r\cos\theta}$

converts the Laplacian into

$\displaystyle{\triangle^{(6)}=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}+\frac{\partial^2}{\partial w^2}+\frac{\partial^2}{\partial s^2}+\frac{\partial^2}{\partial t^2}\\\quad\ \ \ \, ={{\partial^2}\over{\partial r^2}}+{{5}\over{r}}{{\partial}\over{\partial r}}+{{1}\over{r^2}}{{\partial^2}\over{\partial\theta^2}}+{{4\cos \theta}\over{r^2\sin \theta}}{{\partial}\over{\partial\theta}}+{{1}\over{r^2\sin ^2\theta}}{{\partial^2}\over{\partial\varphi^2}}+{{3\cos \varphi}\over{r^2\sin^2\theta\sin\varphi}}{{\partial}\over{\partial\varphi}}\\\qquad\quad+{{1}\over{r^2\sin^2\theta\sin^2\varphi}}{{\partial^2}\over{\partial\psi^2}}+{{2\cos\psi}\over{r^2\sin^2\theta\sin^2\varphi\sin\psi}}{{\partial}\over{\partial\psi}}\\\qquad\quad+{{1}\over{r^2\sin ^2\theta\sin ^2\varphi\sin ^2\psi}}{{\partial^2}\over{\partial\chi^2}}+{{\cos\chi}\over{r^2\sin ^2\theta\sin ^2\varphi\sin ^2\psi\sin\chi}}{{\partial}\over{\partial\chi}}\\\qquad\quad+{{1}\over{r^2\sin ^2\theta\sin ^2\varphi\sin ^2\psi\sin ^2\chi}}{{\partial^2}\over{\partial\eta^2}}}$

Now, can you guess $n$ -dim. case??

Keywords: Python, SymPy, Vector Calculus, Laplacian