Python（SymPy）を使って極座標のラプラシアンを求める

今回は、Python（SymPy）を勉強に活用する実践例として、Python（SymPy）で極座標のラプラシアンを求める方法を紹介したいと思います。２次元と３次元のラプラシアンは電磁気学や量子力学などではおなじみです。物理学の教科書などで扱う問題の多くは回転対称・球対称な場合なので、極座標に移って問題を議論するのはごく自然なことです。

しかし、デカルト座標で表せばあんなに簡単な形に書けるラプラシアンが、ひとたび極座標へ座標変換するやいなやその見た目を大きく変えます。極座標での表式を導出するのにもちょっと手間がかかって、あっちのサインとこっちのコサインがキャンセルしたり組み合わさったり、すったもんだの末に授業が１コマ終わってたりします。

このような、単純だけれども面倒な計算にこそコンピュータを使ってラクをしたいものです。PythonのSymPyは、そのようなことを可能にする数式処理ライブラリです。

次の記事の中ではPython（SymPy）の使い方のごく初歩的な部分を紹介しています：

pianofisica.hatenablog.com

本記事では、２次元、３次元のラプラシアンの極座標表示をPython（SymPy）を使って求めます。おまけとして、４次元の極座標で表したラプラシアンも求めてみます。５次元、６次元の場合についても結果だけ載せておきます。

（考察する対象を容易に拡張できるのもコンピュータを使うメリットです）

２次元ラプラシアン
３次元ラプラシアン
４次元ラプラシアン
５次元ラプラシアン
６次元ラプラシアン

２次元ラプラシアン

まずはウォーミングアップとして２次元のラプラシアン

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

を２次元極座標 $(r,\theta)$ で表してみましょう。 $(x,y)$ と $(r,\theta)$ の関係は

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

または

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

で与えられます。さて、微分の連鎖律から

$\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)}$

です。ここで行列 $R_2$ を

$\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)}$

によって導入しました。この行列は

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

とPythonで入力することにより

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

と計算されます。ここで $f=f(r,\theta)$ を勝手な関数として

$\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)}$

とすると

$\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}}$

によってラプラシアンを求めることができます：

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)))

に入力して上の数式をPythonで計算させることにより

$\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}}$

がわかります。

３次元ラプラシアン

次に３次元のラプラシアン

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

を３次元極座標 $(r,\theta,\varphi)$ で表してみましょう。 $(x,y,z)$ と $(r,\theta,\varphi)$ の関係は

$x=r\sin\theta\cos\varphi\\ y=r\sin\theta\sin\varphi \\ z=r\cos\theta$

または

$\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)}$

で与えられます。ここで２次元のときと同様に $f=f(r,\theta,\varphi)$ を勝手な関数として

$\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)}$

を導入します。すると

$\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}$

によってラプラシアンを求めることができます。まず変数を初期化しておきます

%reset

２次元の場合と同様に

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)))

を入力することにより

$\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}}$

がわかります。

４次元ラプラシアン

ついでに４次元のラプラシアン

$\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}}$

を４次元極座標 $(r,\theta,\varphi,\psi)$ で表してみましょう。 $(x,y,z,w)$ と $(r,\theta,\varphi,\psi)$ の関係は

$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$

または

$\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)}$

です。なんとなくパターンが見えてきますね。あとは省略しますが

%reset -f   #  -f は定義したシンボルを初期化する確認が不要である場合のオプションです

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)))

から

$\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}} }$

と求まります。

以上の結果から５次元、６次元、…の場合の結果を予想してみて、その予想が正しかったかどうか、Pythonを使って検証してみるというのも面白いかもしれません。

ちなみに、答えだけ書いておくと

５次元ラプラシアン

５次元の極座標

$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$

を使って表した５次元のラプラシアンは

$\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}}}$

です。ついでに、６次元のラプラシアンは

６次元ラプラシアン

極座標を

$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$

として

$\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}}}$

となります。一般の $n$ 次元の場合は、…もう予想できますね。

キーワード：Python、SymPy、ベクトル解析、ラプラシアン