Pythonで2次元チューリングパターンを陽解法で計算して描画してみました

考える系や前提条件などの説明は下記の記事を参照ください。

拡散方程式

p(x,y,t+dt)=p(x,y,t)+\Biggl(f(p(x,y,t),q(x,y,t))+\frac{d_p}{dx^2}\Bigl(p(x+dx,y,t)+p(x-dx,y,t)+p(x,y+dy,t)+p(x,y-dy,t)-4p(x,y,t)\Bigr)\Biggr)dt\\
q(x,y,t+dt)=q(x,y,t)+\Biggl(g(p(x,y,t),q(x,y,t))+\frac{d_q}{dx^2}\Bigl(q(x+dx,y,t)+q(x-dx,y,t)+q(x,y+dy,t)+q(x,y-dy,t)-4q(x,y,t)\Bigl)\Biggr)dt

ここで、$d_p,d_q$はそれぞれ活性因子、抑制因子の拡散係数。境界条件は周期境界条件を用います。

定数と関数の宣言

import numpy as np
import matplotlib.pyplot as plt

dx=0.04
dt=0.02
N=int(1/dx)

Dp=2*10**-4
Dq=0.01

def f(a,b):
    return 0.6*a-1.0*b-a*a*a

def g(a,b):
    return 1.5*a-2.0*b

def laplacian(v):
    res=np.roll(v,1,axis=0)\
       +np.roll(v,-1,axis=0)\
       +np.roll(v,1,axis=1)\
       +np.roll(v,-1,axis=1)\
       -4*v
    
    return res

計算セル

初期濃度にはrand()を使ってある程度の揺らぎを与えておきます。

%%time
p=np.random.rand(N,N,1)*0.01
q=np.random.rand(N,N,1)*0.01

tmax=2000
for i in range(tmax):
    next_p=p[:,:,-1]+dt*(f(p[:,:,-1],q[:,:,-1])+Dp*laplacian(p[:,:,-1])/(dx*dx))
    next_q=q[:,:,-1]+dt*(g(p[:,:,-1],q[:,:,-1])+Dq*laplacian(q[:,:,-1])/(dx*dx))
    p=np.dstack([p,next_p])
    q=np.dstack([q,next_q])

CPU times: user 2.47 s, sys: 2.89 s, total: 5.36 s
Wall time: 5.39 s

Gif作成用

%matplotlib nbaggはJupyter notebook上で動画を表示するために必要です。
容量削減のために10枚ごとに描画しています。

%matplotlib nbagg
from matplotlib.animation import PillowWriter,FuncAnimation

def animate(i):
    clock=dt*i*10
    step=i*10
    plt.cla()
    plt.imshow(p[:,:,i*10],interpolation='nearest',vmin=-1,vmax=1,cmap='jet')
    
    plt.title("t={0:.1f} ,step={1:3d}".format(clock,step))

fig=plt.figure()
plt.colorbar(plt.imshow(p[:,:,-1],interpolation='nearest',vmin=-1,vmax=1,cmap='jet'))

anim=FuncAnimation(fig,animate,repeat=False,frames=200,interval=10)
#anim.save("2D_Turing_Pattern.gif", writer='pillow',fps=10) #Gifを保存する場合はコメントアウトを外す
#plt.show() #Gifを表示する場合はコメントアウトを外す