ブリュッセレータの分岐

• 存在から発展へ
  • きのうの記事 のつづき
  • p119
• 2次元でブリュッセレータに拡散があるモデル
  • こちら の反応を こちら の場合のように2次元で考える
• おそらく1マスの大きさによって拡散係数は変わってくうと思うので、本に出てくる図と同じになる条件を探索中
  • 拡散係数などを修正する
  • 図 5.7 は対称性の保たれたモデル（こちらへ）
  • 図 5.8 は対称性のくずれたモデル（こちらへ）

library(rgl)

A<-2
B<-4.6
X0<-A
Y0<-B/A

Lx<-0.2
dx<-0.01
Nx<-floor(Lx/dx)+1

Ly<-0.2
dy<-0.01
Ny<-floor(Ly/dy)+1

T<-30
dt<-0.001
Nt<-T/dt

# 拡散係数
Dx<-0.00325
Dy<-0.0162


# 円を作る
K<-matrix(0,Nx,Ny)
for(nx in 1:Nx){
for(ny in 1:Ny){
if((nx-Nx/2)^2/(Nx/2)^2+(ny-Ny/2)^2/(Ny/2)^2 < 0.9){K[nx,ny]<-1}
}}

# 内部は1、境界は2
K[-Nx,][which(K[-Nx,]-K[-1,]!=0)]<-K[,-Ny][which(K[,-Ny]-K[,-1]!=0)]<-2
K[1,]<-(K[1,]>0)*2;K[Nx,]<-(K[Nx,]>0)*2;K[,1]<-(K[,1]>0)*2;K[,Ny]<-(K[,Ny]>0)*2


# 各点のX,Yの濃度
C<-array(c(X0),c(Nx,Ny,2))
# C[1:floor(Nx/2),,2]<-Y0+0.001  #初期値にゆらぎを与える
 C[,,2]<-rep(Y0,Nx*Ny)+runif(Nx*Ny)*0.1

C[,,1]<-C[,,1]*(K>0)
C[,,2]<-C[,,2]*(K>0)

C[,,1][which(K==2)]<-X0
C[,,2][which(K==2)]<-Y0

# 拡散係数の計算
DX<-array(0,c(Nx-1,Ny,2))
DY<-array(0,c(Nx,Ny-1,2))
DX[,,1]<-(K[-1,]*K[-Nx,]>0)*Dx
DX[,,2]<-(K[-1,]*K[-Nx,]>0)*Dy
DY[,,1]<-(K[,-1]*K[,-Ny]>0)*Dx
DY[,,2]<-(K[,-1]*K[,-Ny]>0)*Dy

# ブリュッセレータを計算する関数
f1<-function(X){
	y<- A+X[1]^2*X[2]-B*X[1]-X[1]
	return(y)
	}
f2<-function(X){
	y<- B*X[1]-X[1]^2*X[2]
	return(y)
	}

# ここからがシミュレーション	
for(t in 1:Nt)
{	
	# ブリュッセレータの計算
	for(nx in 1:Nx){
	for(ny in 1:Ny)
	if(K[nx,ny] >= 1){
	X<-C[nx,ny,]
	dX<-c(f1(X),f2(X))*dt
	C[nx,ny,]<-X+dX
	}
	}
	
	# 拡散の計算
	dC_x<-C[-1,,]-C[-Nx,,]
	dC_y<-C[,-1,]-C[,-Ny,]

	Ix<- -DX*dC_x
	Iy<- -DY*dC_y

	C[-1,,]<-C[-1,,]+Ix
	C[-Nx,,]<-C[-Nx,,]-Ix
	C[,-1,]<- C[,-1,]+Iy
	C[,-Ny,]<-C[,-Ny,]-Iy
	
	# 固定端
#	C[,,1][which(K==2)]<-X0
#	C[,,2][which(K==2)]<-Y0
		
	if(t %% 100 ==0){
#	persp(C[,,1],theta=30,phi=45,col=4,xlab="x",ylab="y",zlab="X conc.",zlim=c(0,4),main=paste("time = ",t*dt,sep=""))
#	image(C[,,1])
	plot3d(slice.index(C[,,1],1),slice.index(C[,,1],2),C[,,1],alpha=ifelse(K==0,0,1),xlab="x",ylab="y",zlab="X conc.",col=rainbow(101)[C[,,1]/max(C[,,1])*100+1])
	}
	
}