\(U^n_{i+j}\)に関する時間発展方程式 \[U^{n+1}_{ij}=\max(U^{n}_{ij},U^{n}_{i-1j},U^{n}_{i+1j},U^{n}_{ij-1},U^{n}_{ij+1},U^{n-1}_{ij})-U^{n-1}_{ij}\;\tag{1}\] を考える。 \(i,j\)は2次元の空間の格子点の座標を表し,\(n\)(整数)は時刻を表す。差分格子は\(-\infty<i<\infty,-\infty<j<\infty\)の範囲とする。 この式において、次の時刻\(n+1\)の値を得るためには\(n-1,n\)での\(U\)の値が必要である。よって、初期値としては連続する2つの時刻での\(U\)の空間分布を与える。

ここでは、\(-30\leqq i\leqq30,-30\leqq j\leqq 30,1\leqq n\leqq 100\)の範囲で上の方程式を考えて、得られたパターンを見ていく。

例

初期値を\(U_{0,0}^0=U_{0,0}^1=1\)として時間発展を調べると次のようになる。(ターゲットパターン)

t_end <- 100
m <- 30

mat <- array(0,dim=c(t_end,2*m+1,2*m+1)) #60*60の値が0の格子点を用意する。

ind <- function(x){x+m+1}
f <- function(t,x,y){
  if(x< -m || x>m || y < -m || y>m){return(0)}else{return(mat[t,ind(x),ind(y)])}
}

mat[1,ind(0),ind(0)] <- 1 #初期値を定める。
mat[2,ind(0),ind(0)] <- 1 
for(t in  3:t_end){
  for(x in -m:m){
    for(y in -m:m){
      mat[t,ind(x),ind(y)] <- max(f(t-1,x,y),f(t-1,x-1,y),f(t-1,x+1,y),f(t-1,x,y-1),f(t-1,x,y+1),f(t-2,x,y))-f(t-2,x,y) #時間発展方程式(1)
    }
  }
}

時間発展の様子は次のようになる。 \(n=1\)のとき

image(mat[1,,])

\(n=2\)のとき

image(mat[2,,])

\(n=3\)のとき

image(mat[3,,])

\(n=4\)のとき

image(mat[4,,])

\(n=5\)のとき

image(mat[5,,])

\(n=100\)のとき

image(mat[t_end,,])

スパイラルパターン

BZパターン

2次元パターンの現れる現象の例としてBZ反応があげられる。この反応を表現する拡散方程式は次のようになる。 \[ \begin{cases} \tau u_t=D_u\Delta u+f(u,v)-v\\ u_t=D_v\Delta v+u-v \end{cases} \] ただし、\(f(u,v)=u(1-u)-\dfrac{pv(u-q)}{u+q}\)で\(D_u,D_v\)は\(u,v\)の拡散係数である。 ここで、スパイラルパターンを、(1)の初期値として、\(U_{1i}^0=U_{0i}^0=U_{1i}^1=U_{2i}^1=1\;(-5\leqq i\leqq 5)\)とすれば再現できる。

t_end1 <- 50
m <- 30

mat1 <- array(0,dim=c(t_end1,2*m+1,2*m+1))

ind <- function(x){x+m+1}
f1 <- function(t,x,y){
  if(x< -m || x>m || y < -m || y>m){return(0)}else{return(mat1[t,ind(x),ind(y)])}
}

for(i in -5:5){
  mat1[1,ind(0),ind(i)] <- 1
  mat1[1,ind(1),ind(i)] <- 1
  mat1[2,ind(1),ind(i)] <- 1
  mat1[2,ind(2),ind(i)] <- 1
}

for(t in  3:t_end1){
  for(x in -m:m){
    for(y in -m:m){
      mat1[t,ind(x),ind(y)] <- max(f1(t-1,x,y),f1(t-1,x-1,y),f1(t-1,x+1,y),f1(t-1,x,y-1),f1(t-1,x,y+1),f1(t-2,x,y))-f1(t-2,x,y)
    }
  }
}
for(i in 1:t_end1){
  image(mat1[i,,])
}

スパイラルパターンは(1)より、\(U\)に関する時間2階差分方程式である。2時刻分の初期値を用意すれば時間発展を追うことができる。ここで、周囲の値に影響されずに時間発展する局所的なパターンをコアパターンという。

生物の増殖

(1)において、\(V_{ij}^n=U_{ij}^{n-1}\)とすれば、(1)は\(U,V\)の連立方程式に書き換えられる。 \[ \begin{cases} U^{n+1}_{ij}=\max(U^{n}_{ij},U^{n}_{i-1j},U^{n}_{i+1j},U^{n}_{ij-1},U^{n}_{ij+1},V^{n}_{ij})-V^{n}_{ij}\\ V_{i,j}^{n+1}=U_{ij}^n \tag{2} \end{cases} \] \(U_{ij}^n,V_{ij}^n\)を時刻\(n\)における\((i,j)\)格子点の生物種A,Bの個体数だとすれば、(2)は増え続けるAをBがたべるという生存競争の単純なモデルとみなすことができる。 生物種Aの増殖をみるには\(V^n_{ij}=0\)とすればよい。 \[U^{n+1}_{ij}=\max(U^{n}_{ij},U^{n}_{i-1j},U^{n}_{i+1j},U^{n}_{ij-1},U^{n}_{ij+1})\] この式において、\(U^0_{0,0}\)だけを1として残りを0とすれば、Aの増殖を再現できる。

t_end2 <- 50
m <- 30

mat2 <- array(0,dim=c(t_end2,2*m+1,2*m+1))

ind <- function(x){x+m+1}
f2 <- function(t,x,y){
  if(x< -m || x>m || y < -m || y>m){return(0)}else{return(mat2[t,ind(x),ind(y)])}
}

mat2[1,ind(0),ind(0)] <- 1
for(t in  2:t_end2){
  for(x in -m:m){
    for(y in -m:m){
      mat2[t,ind(x),ind(y)] <- max(f2(t-1,x,y),f2(t-1,x-1,y),f2(t-1,x+1,y),f2(t-1,x,y-1),f2(t-1,x,y+1))
    }
  }
}
for(i in 1:t_end2){
  image(mat2[i,,])
}

　参考文献　

「あるパターン方程式のダイナミクス　Dynamics of Pattern Formation Equations（志田篤彦, 高橋大輔）」　（早稲田大学理工学部数理科学科）