pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

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

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

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

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

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

pianofisica.hatenablog.com

本記事では、2次元、3次元のラプラシアン極座標表示をPython(SymPy)を使って求めます。おまけとして、4次元の極座標で表したラプラシアンも求めてみます。5次元、6次元の場合についても結果だけ載せておきます。

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



2次元ラプラシアン

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

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

を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}}

がわかります。



3次元ラプラシアン

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

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

を3次元極座標  (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)}

で与えられます。ここで2次元のときと同様に  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

2次元の場合と同様に

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

がわかります。

 


4次元ラプラシアン

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

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

を4次元極座標  (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}}  }

と求まります。


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



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

5次元ラプラシアン

5次元の極座標

 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

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

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

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

6次元ラプラシアン

極座標

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、ベクトル解析、ラプラシアン