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

Python + matplotlib.pyplot で3Dランダムウォーク

Python でランダムウォークくらいググれば誰でもできる。
そんなことを改まって書くのも無益なので、

3D、軌跡がColormapに従って変化

するものでも書いてみようかと思う。
3Dグラフは普段使わないので忘れてしまうからというのも、もちろんある。
で結果としてはこんなものだ。
さすがに一枚目のグラフは即座に描画できるけど、
二枚目のグラフは一歩ずつ描画しているのでなかなか時間がかかる。

In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

In [2]:

N = 2000
R = (np.random.rand(N)*6).astype("int")
x = np.zeros(N)
y = np.zeros(N)
z = np.zeros(N)
x[ R==0 ] = -1; x[ R==1 ] = 1
y[ R==2 ] = -1; y[ R==3 ] = 1
z[ R==4 ] = -1; z[ R==5 ] = 1
x = np.cumsum(x)
y = np.cumsum(y)
z = np.cumsum(z)

In [3]:

plt.figure()
ax = plt.subplot(1,1,1, projection='3d')
ax.plot(x, y, z,alpha=0.6)
ax.scatter(x[-1],y[-1],z[-1])
plt.show()

In [4]:

plt.figure()
ax = plt.subplot(1,1,1, projection='3d')
cm = plt.get_cmap('jet')
ax.set_prop_cycle('color',[cm(1.*i/(x.shape[-1]-1)) for i in range(x.shape[-1]-1)])
for i in range(x.shape[-1]-1):
    ax.plot([x[i+1],x[i]], [y[i+1],y[i]], [z[i+1],z[i]],alpha=0.6)
ax.scatter(x[-1],y[-1],z[-1],facecolor=cm(1))
plt.show()