Python Sympy で極座標のラプラシアン

恥ずかしながら、Sympyは使ったことがなかった。

座標変換、代数演算の微分積分をする必要があって、
Walfram alpha を使ったら、
思いの外便利で高性能だった。
Walfram alpha だと複数行の入力はできないので、
もっと汎用なものがPythonでできればいい。

ということで、Sympyを触ってみた健忘録だ。
Sympy は matplotlib.pyplot と同じように、
当然のように Jupyter notebook との相性が抜群だ。

今回は、極座標のラプラシアンを求めてみた。

In [1]:

import sympy as sy
sy.init_printing()

In [2]:

r = sy.Symbol('r',real=True, positive=True)
th = sy.Symbol('theta',real=True, positive=True,domain=sy.Interval(0,sy.pi))
ph = sy.Symbol('\phi',real=True, positive=True,domain=sy.Interval(0,2*sy.pi))

In [3]:

x = r * sy.sin(th) * sy.cos(ph)
y = r * sy.sin(th) * sy.sin(ph)
z = r * sy.cos(th)

In [4]:

xyz = [x, y, z]
xyz

Out[4]:

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

In [5]:

rtp = [r, th, ph]
rtp

Out[5]:

$$\left [ r, \quad \theta, \quad \phi\right ]$$

In [6]:

f = sy.Function('f')(r , th , ph)
f

Out[6]:

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

スケールファクター¶

In [7]:

hr = sy.simplify(sy.sqrt(sy.diff(x,r)**2+sy.diff(y,r)**2+sy.diff(z,r)**2))
ht = sy.simplify(sy.sqrt(sy.diff(x,th)**2+sy.diff(y,th)**2+sy.diff(z,th)**2))
hp = sy.simplify(sy.sqrt(sy.diff(x,ph)**2+sy.diff(y,ph)**2+sy.diff(z,ph)**2))

In [8]:

hr,ht,hp

Out[8]:

$$\left ( 1, \quad r, \quad r \left|{\sin{\left (\theta \right )}}\right|\right )$$

ヤコビアンとその逆行列¶

In [9]:

# ヤコビアン
J = sy.Matrix(3 , 3 , [sy.Derivative(c1 , c2) for c2 in rtp for c1 in xyz])
J

Out[9]:

$$\left[\begin{matrix}\frac{\partial}{\partial r}\left(r \sin{\left (\theta \right )} \cos{\left (\phi \right )}\right) & \frac{\partial}{\partial r}\left(r \sin{\left (\phi \right )} \sin{\left (\theta \right )}\right) & \frac{\partial}{\partial r}\left(r \cos{\left (\theta \right )}\right)\\\frac{\partial}{\partial \theta}\left(r \sin{\left (\theta \right )} \cos{\left (\phi \right )}\right) & \frac{\partial}{\partial \theta}\left(r \sin{\left (\phi \right )} \sin{\left (\theta \right )}\right) & \frac{\partial}{\partial \theta}\left(r \cos{\left (\theta \right )}\right)\\\frac{\partial}{\partial \phi}\left(r \sin{\left (\theta \right )} \cos{\left (\phi \right )}\right) & \frac{\partial}{\partial \phi}\left(r \sin{\left (\phi \right )} \sin{\left (\theta \right )}\right) & \frac{\partial}{\partial \phi}\left(r \cos{\left (\theta \right )}\right)\end{matrix}\right]$$

In [10]:

# 計算しておく
J = J.doit()
J

Out[10]:

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

In [11]:

# ヤコビアンの逆行列
InvJ = sy.simplify(J.inv())
InvJ

Out[11]:

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

In [12]:

a = sy.Matrix(3 , 1 , [sy.Derivative(f , c) for c in rtp])
a

Out[12]:

$$\left[\begin{matrix}\frac{\partial}{\partial r} f{\left (r,\theta,\phi \right )}\\\frac{\partial}{\partial \theta} f{\left (r,\theta,\phi \right )}\\\frac{\partial}{\partial \phi} f{\left (r,\theta,\phi \right )}\end{matrix}\right]$$

In [13]:

# grad(f)
gradf = InvJ*a
gradf

Out[13]:

$$\left[\begin{matrix}\sin{\left (\theta \right )} \cos{\left (\phi \right )} \frac{\partial}{\partial r} f{\left (r,\theta,\phi \right )} - \frac{\sin{\left (\phi \right )} \frac{\partial}{\partial \phi} f{\left (r,\theta,\phi \right )}}{r \sin{\left (\theta \right )}} + \frac{1}{r} \cos{\left (\phi \right )} \cos{\left (\theta \right )} \frac{\partial}{\partial \theta} f{\left (r,\theta,\phi \right )}\\\sin{\left (\phi \right )} \sin{\left (\theta \right )} \frac{\partial}{\partial r} f{\left (r,\theta,\phi \right )} + \frac{1}{r} \sin{\left (\phi \right )} \cos{\left (\theta \right )} \frac{\partial}{\partial \theta} f{\left (r,\theta,\phi \right )} + \frac{\cos{\left (\phi \right )} \frac{\partial}{\partial \phi} f{\left (r,\theta,\phi \right )}}{r \sin{\left (\theta \right )}}\\\cos{\left (\theta \right )} \frac{\partial}{\partial r} f{\left (r,\theta,\phi \right )} - \frac{1}{r} \sin{\left (\theta \right )} \frac{\partial}{\partial \theta} f{\left (r,\theta,\phi \right )}\end{matrix}\right]$$

In [14]:

dxx = (InvJ * sy.Matrix(3,1,[ sy.Derivative(gradf[0],e) for e in rtp ]).doit())[0].doit()
dxx

Out[14]:

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

In [15]:

dyy = (InvJ * sy.Matrix(3,1,[ sy.Derivative(gradf[1],e) for e in rtp ]).doit())[1].doit()
dyy

Out[15]:

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

In [16]:

dzz = (InvJ * sy.Matrix(3,1,[ sy.Derivative(gradf[2],e) for e in rtp ]).doit())[2].doit()
dzz

Out[16]:

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

In [17]:

# ラプラシアン
sy.simplify(dxx+dyy+dzz)

Out[17]:

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

In [18]:

# 一行で書くと
sy.simplify(sum([ (InvJ * sy.Matrix(3,1,[ sy.Derivative(gradf[i],e) for e in rtp ]).doit())[i] for i in range(3) ]))

Out[18]:

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