チューリングパターン 反応拡散系

• こちらのpythonをRに直す

a <- 2.8 * 10^(-4)
b <- 5 * 10^(-3)
tau <- 0.1
k <- -0.005

size <- 100
dx <- 2/size
T <- 10.0
dt <- 0.9 * dx^2/2
n <- floor(T/dt)

U <- matrix(runif(size^2),size,size)
V <- matrix(runif(size^2),size,size)

ctr <- 2:(size-1)
plus <- 3:size
minus <- 1:(size-2)

U.hx <- list()
V.hx <- list()
U.hx[[1]] <- U
V.hx[[1]] <- V
for(i in 1:n){
  deltaU <- (U[minus,ctr] + U[plus,ctr] + U[ctr,minus] + U[ctr,plus] - 4 * U[ctr,ctr])/(dx^2)
  deltaV <- (V[minus,ctr] + V[plus,ctr] + V[ctr,minus] + V[ctr,plus] - 4 * V[ctr,ctr])/(dx^2)
  
  Uctr <- U[ctr,ctr]
  Vctr <- V[ctr,ctr]
  
  new.U <- U
  new.V <- V
  new.U[ctr,ctr] <- Uctr + dt * (a*deltaU + Uctr - Uctr^3 -Vctr + k)
  new.V[ctr,ctr] <- Vctr + dt * (b*deltaV + Uctr - Vctr)/tau
  
  new.U[1,] <- new.U[2,]
  new.U[size,] <- new.U[size-1,]
  new.U[,1] <- new.U[,2]
  new.U[,size] <- new.U[,size-1]
  new.V[1,] <- new.V[2,]
  new.V[size,] <- new.V[size-1,]
  new.V[,1] <- new.V[,2]
  new.V[,size] <- new.V[,size-1]
  U <- new.U
  V <- new.V
  
  #U.hx[[i+1]] <- U
  #V.hx[[i+1]] <- V
  
  if(i%%5000 == 1){
		image(U,col=topo.colors(100))
	}
}