本記事の目的

・反応拡散系におけるGray-Scottモデルという自己触媒反応の数理モデルを2次元の正方形の空間にてPython3によるシミュレートを行う．

用いた環境

・Python3.7
・MacOS Mojave

扱う系

今回考える2次元の正方形の空間とは具体的にいうと，2次元の平面の領域において$dx$という有限の長さに分割(離散化)する．また時間も離散化され，0以上の整数値のみをとるものとする．

陽解法

今回用いるGray-Scottモデルは以下の数式で表され，陽解法を用いて数値計算を行います．

\begin{array}{l}{
\frac{\partial u}{\partial t}=D_{u} \Delta v -u^{2}v + F(1-v)} \\ 
{\frac{\partial v}{\partial t}=D_{v} \Delta u +u^{2}v -(F+K) u}
\end{array}

ある時刻$t$において，$u$と$v$は濃度で$u = u(x,y,t)$，$v = v(x,y,t)$と表すことにする．$D_{u}$，$D_{v}$は拡散係数であり，係数$F, K$は0以上の実数とする．2式の右辺第一項は拡散項と呼ばれ，第二，三項は反応項と呼ばれるものです．

数値計算で解く式が以下の2式になります．

u(x,y,t+dt) = u(x,y,t) + dt \Biggl(\frac{D_{u}}{dx^{2}}\biggr(u(x+dx,y,t)+u(x−dx,y,t)+u(x,y+dy,t)+u(x,y−dy,t)−4u(x,y,t) \biggr) −u(x,y,t)^{2}v(x,y,t)+F(1−v(x,y,t)) \Biggr) \\
v(x,y,t+dt) = v(x,y,t) + dt \Biggr(\frac{D_{v}}{dx^{2}}\biggr( v(x+dx,y,t)+v(x−dx,y,t)+v(x,y+dy,t)+v(x,y−dy,t)−4v(x,y,t) \biggr) +u(x,y,t)^{2}v(x,y,t)−(F+K)u(x,y,t)) \Biggr)

ただし，$dx, dy$は空間方向における有限値の差分間隔とします．

用いたライブラリ

import numpy as np
import matplotlib.pyplot as plt
import math

各パラメータの設定

ここの右辺に4種のパラメータを代入します．

Du, Dv, F, k = 0.14, 0.06, 0.035, 0.065

数値計算と描画のプログラム

以下に，数値計算とグレースケールでの描画をするプログラムです．Vの値を0から1までのグレースケールで出力するように設定しました．また，1万回計算するようにプログラムを組んでいるので，実行時間が少々かかるかと思います．

n  = 200
Z = np.zeros((n+2,n+2), [('U', np.double), ('V', np.double)] )
U,V = Z['U'], Z['V']
u,v = U[1:-1,1:-1], V[1:-1,1:-1]#それぞれ-1から１までを範囲とする

r = 20
u[...] = 1.0

U[n//2-r:n//2+r, n//2-r:n//2+r] = 0.50
V[n//2-r:n//2+r, n//2-r:n//2+r] = 0.25

#初期濃度はランダム関数で与える
u += 0.05*np.random.random((n,n))
v += 0.05*np.random.random((n,n))

plt.ion()

size = np.array(Z.shape)
dpi = 72.0
figsize= size[1]/float(dpi),size[0]/float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi, facecolor="white")
fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=False)
im = plt.imshow(V, interpolation='bicubic', cmap=plt.cm.gray_r)
plt.xticks([])
plt.yticks([])


for i in range(10000):#陽解法
    Lu = ( U[0:-2,1:-1] + U[1:-1,0:-2] - 4*U[1:-1,1:-1] + U[1:-1,2:] + U[2:  ,1:-1] )
    Lv = ( V[0:-2,1:-1] + V[1:-1,0:-2] - 4*V[1:-1,1:-1] + V[1:-1,2:] + V[2:  ,1:-1] )

    u += (Du*Lu - u*v*v +  F  *(1-u) )#解く式1
    v += (Dv*Lv + u*v*v - (F+k)* v   )#解く式2

im.set_data(V)#Vを出力
im.set_clim(vmin=V.min(), vmax=V.max())
plt.draw()

plt.ioff()   
plt.show()

プログラム全体

main.py

# -*- coding: utf-8 -*-
#python3に対応化@2019/3/24

import numpy as np
import matplotlib.pyplot as plt
import math

#パラメータはここを参考にした。
# from http://www.aliensaint.com/uo/java/rd/
# -----------------------------------------------------

Du, Dv, F, k = 0.16, 0.08, 0.035, 0.065

n  = 200
Z = np.zeros((n+2,n+2), [('U', np.double), ('V', np.double)] )
U,V = Z['U'], Z['V']
u,v = U[1:-1,1:-1], V[1:-1,1:-1]#それぞれ-1から１までを範囲とする

r = 20
u[...] = 1.0

U[n//2-r:n//2+r, n//2-r:n//2+r] = 0.50
V[n//2-r:n//2+r, n//2-r:n//2+r] = 0.25

#初期濃度はランダム関数で与える
u += 0.05*np.random.random((n,n))
v += 0.05*np.random.random((n,n))

plt.ion()

size = np.array(Z.shape)
dpi = 72.0
figsize= size[1]/float(dpi),size[0]/float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi, facecolor="white")
fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=False)
im = plt.imshow(V, interpolation='bicubic', cmap=plt.cm.gray_r)
plt.xticks([])
plt.yticks([])


for i in range(10000):#陽解法
    Lu = ( U[0:-2,1:-1] + U[1:-1,0:-2] - 4*U[1:-1,1:-1] + U[1:-1,2:] + U[2:  ,1:-1] )
    Lv = ( V[0:-2,1:-1] + V[1:-1,0:-2] - 4*V[1:-1,1:-1] + V[1:-1,2:] + V[2:  ,1:-1] )

    u += (Du*Lu - u*v*v +  F  *(1-u) )##解く式1
    v += (Dv*Lv + u*v*v - (F+k)* v   )##解く式2
im.set_data(V)#Vを出力
im.set_clim(vmin=V.min(), vmax=V.max())
plt.draw()

plt.ioff()   
plt.show()