やること
pythonでラプラシアンを含む微分方程式のシミュレーションをしたかったのですが,つい先日にラプラシアンの差分をfor文以外で処理する方法を思いついた(車輪の再発明をした)ので拡散方程式のアニメを書きました.
思いついてから改めて調べてみたら数年前に自分のやったことの上位互換をされてる方がいらっしゃいました.その記事はかなり参考にさせていただきました: NumPy・SciPyを用いた数値計算の高速化 : 応用その1
解いている微分方程式は次式で表される1次元拡散方程式です:
$$
\frac{\partial }{\partial t}\psi(x,t)= D\frac{\partial^2}{\partial x^2}\psi(x,t)
$$
初期条件はガウシアンで与え,時間発展はとりあえずオイラー法にしました.
書いたコード
%matplotlib nbagg
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import ArtistAnimation
# calculate Laplacian with periodic boundary condition
def Lap(L, M):
# expand vector
VL = [L[0]]
VR = [L[-1]]
Le = np.hstack([VR ,L, VL])
Lc = np.dot(M, Le) # central diff
return Lc[1:N+1]
# generate 1-D lattice and initialize
N = 64
x = np.linspace(-10, 10, N)
psi = np.exp(-x**2/7)
# time and diffusion const
dt = 0.05
T = np.arange(0.0, 37, dt)
D = 3
# make (N+2)x(N+2) matrix to calc Laplacian
I = np.eye(N+2,N+2)
S = np.vstack((I[1:],np.array([0]*(N+2))))
M = S + S.T -2*I
# make animation
fig = plt.figure()
ims = []
for i in range(len(T)):
psi += dt*(D*Lap(psi, M))
line = plt.plot(x, psi, 'b')
ims.append(line)
ani = ArtistAnimation(fig, ims, interval=5, blit=True, repeat=True)
plt.show()
実行例
工夫した点
ラプラシアンの計算を行う際に元の$N$次元配列$\psi$の右端と左端の値をそれぞれ左端と右端にnp.hstack
でくっつけて$N+2$次元配列にしました.これにより,予め用意した$(N+2)\times(N+2)$行列$M$と内積を取れば自動的に周期境界条件を含む拡散項が計算できました.境界条件用に追加した部分は不要なので関数の最後に捨てています.
総括
課題点:とりあえず見える形にしたかったので,数値的な安定性や誤差は1mmも考慮していません.NumPyに高速かつ安全に数値積分してくれる何かがあるとは思うので調べて何とかしようと思います.これらの課題に関して,今後何らかの進捗が生まれたら追記します.
今後:今回の1次元バージョンを適当に拡張すれば2次元バージョンも作れそうなので,反応拡散系のシミュレーションをやってみたいなあと考えています.(よっぽどエグい非線形項を含む方程式でなければ何とかなるっしょ...)
追記:
(11/23) 前進差分で妥協すれば,3行で非線形項を入れられました.一応1次元のNavier-Stokes方程式もできそうです.
def nl(L, N):
Lf = np.hstack((L,L[0]))
return L * np.diff(Lf)
前進差分だとレイノルズ数のみならず初期条件や時間の刻み幅も慎重に選ばないとすぐに計算が発散しました.さぁ,どうしましょうか