分岐：ブリュッセレーター

R シミュレーション 拡散

KABIRA2011-07-24

  • 存在から発展へ
  • ブリュッセレーターに拡散を考慮したもの
    • \begin{matrix} \frac{dX}{dt} &=& A + X^2Y -BX -X + D_X\frac{\par^2X}{\par r^2}\\ \frac{dY}{dt} &=& BX-X^2Y + D_Y\frac{\par^2Y}{\par r^2 } \end{matrix}
  • A = 2, \, B = 4.17, \, D_X = 0.0016, \, D_Y = 0.008 として p117 の定常状態が得られる

  • 下のソースは定常状態にはならずに波をうつ
    • ムービーを参照
    • p118ページ相当 と思われる
A<-2
B<-6
X0<-A
Y0<-B/A

L<-1
dr<-0.01
Nr<-floor(L/dr)+1

T<-3
dt<-0.0001
Nt<-T/dt

# 拡散係数
Dx<-0.05
Dy<-0.02

# 各点のX,Yの濃度
C<-matrix(c(X0,Y0),2,Nr)
Nint<-2
C<-C+c(sin(Nint*pi*1:Nr/(L*Nr)),sin(Nint*pi*1:Nr/(L*Nr)))

range<-c(min(C),max(C))
plot(seq(from=0,to=1,length=Nr),C[1,],type="l",ylim=range,xlab="r",ylab="X conc.",main="time = 0")
# plot(seq(from=0,to=1,length=Nr),C[2,],type="l",ylim=range,xlab="r",ylab="Y conc.",main="time = 0")

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(l in 1:Nr)
	{
	X<-C[,l]
	dX<-c(f1(X),f2(X))*dt
	C[,l]<-X+dX
	}
	
	# 拡散の計算
	dC_r<-C[,-1]-C[,-Nr]
	Ir<- -c(Dx,Dy)*dC_r
	C[,-1]<-C[,-1]+Ir
	C[,-Nr]<-C[,-Nr]-Ir
	
	# 固定端
	C[,1]<-C[,Nr]<-c(X0,Y0)

	if(t %% 100 == 0){
	plot(seq(from=0,to=1,length=Nr),C[1,],type="l",ylim=range,xlab="r",ylab="X conc.",main=paste("time = ",t*dt,sep=""))
#	plot(seq(from=0,to=1,length=Nr),C[2,],type="l",ylim=range,xlab="r",ylab="Y conc.",main=paste("time = ",t*dt,sep=""))
	}

}