Python3で反応拡散系のパターンを描画する


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

#用いた環境
・Python3.7
 ・MacOS Mojave 


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

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

```math
\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式になります．

```math
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()
```


#プログラム全体

```python: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()
```
このプログラムを実行すると，バクテリアと呼ばれる以下のパターンが形成されます
![Bacteria.png](https://qiita-image-store.s3.amazonaws.com/0/185351/aef2cf89-fe50-2ac9-b704-a46873da8c2f.png)



#他のパラメータでの数値シミュレーション
最後に，他のパラメータで実験してみたところ，4つのパターンが確認できました．
パラメータは以下のサイトを参考に設定しました．
http://www.aliensaint.com/uo/java/rd/

```
Du, Dv, F, k = 0.16, 0.08, 0.060, 0.062
```
Du, Dv, F, k = 0.16, 0.08, 0.060, 0.062の場合だと，Coral Patternと呼ばれるパターンになります
![Coral Pattern.png](https://qiita-image-store.s3.amazonaws.com/0/185351/b2eab264-00c6-11fb-1033-00f9e39ec74d.png)


```
Du, Dv, F, k = 0.19, 0.05, 0.060, 0.062
```
Du, Dv, F, k = 0.19, 0.05, 0.060, 0.062の場合だと，Fingerprintと呼ばれるパターンになります
![Fingerprint.png](https://qiita-image-store.s3.amazonaws.com/0/185351/9669ef70-984f-5fc9-656a-6022d315fbc2.png)

```
Du, Dv, F, k = 0.12, 0.08, 0.020, 0.050
```
Du, Dv, F, k = 0.12, 0.08, 0.020, 0.050の場合だと，Spirals Denseと呼ばれるパターンになります
![Spirals Dense.png](https://qiita-image-store.s3.amazonaws.com/0/185351/faddefa3-42e2-1e19-23ac-fb400b54be42.png)

```
Du, Dv, F, k = 0.16, 0.08, 0.035, 0.060
```
Du, Dv, F, k = 0.16, 0.08, 0.035, 0.060の場合だと，Zebrafishと呼ばれるパターンになります
![Zebrafish.png](https://qiita-image-store.s3.amazonaws.com/0/185351/b9450b82-4127-fb2e-2fee-4c963be57f25.png)


#参考文献．
・2次元Turingパターンの生成
https://qiita.com/STInverSpinel/items/b9d844c98e63dfe72adc

反応拡散系そのものの説明は次の2つの文献が詳しいと思います
・Mechanism of pigment pattern formation
https://www.fbs-osaka-kondolabo.net/pigment-pattern

・反応拡散系（Reaction-Diffusion system）
http://www.fbs.osaka-u.ac.jp/labs/skondo/ozaki/what%20is%20RD%202(outline).htm