Python3에서 반응 확산 시스템의 패턴 그리기

18499 단어 수치 계산 Python3

이 기사의 목적

반응 확산 시스템에서의 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식의 우변 제1항은 확산항이라고 하고, 제2, 3항은 반응항이라고 부른다.

수치 계산으로 풀리는 식이 다음의 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()