---
title: "モデル方程式"
author: "山本　裕樹"
date: "2018年1月15日"
output: html_document
---

$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$として時間発展を調べると次のようになる。(ターゲットパターン)

```{r}
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$のとき

```{r}
image(mat[1,,])
```

$n=2$のとき

```{r}
image(mat[2,,])
```

$n=3$のとき

```{r}
image(mat[3,,])
```

$n=4$のとき

```{r}
image(mat[4,,])
```

$n=5$のとき

```{r}
image(mat[5,,])
```

$n=100$のとき

```{r}
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)$とすれば再現できる。
```{r}
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)
    }
  }
}

```

```{r hist1, interval=0.5, fig.height=8, fig.show="animate", fig.width=8, aniopts='controls,autoplay'}
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の増殖を再現できる。

```{r}
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))
    }
  }
}

```

```{r hist2, interval=0.5, fig.height=8, fig.show="animate", fig.width=8, aniopts='controls,autoplay'}
for(i in 1:t_end2){
  image(mat2[i,,])
}
```

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