化学反応 - KABIRAの日記

KABIRA2011-07-20

存在から発展へ
- p104
解が安定なものの例
- それぞれの反応式の反応係数を $k_1, \, k_2, \, k_3$ として
- $\begin{matrix} \frac{dX}{dt} & = & k_1AX &-& k_2XY \\ \frac{dY}{dt} & = & k_2XY &-& k_3Y \end{matrix}$
- 定常解は $\begin{cases} X_0 &=& \frac{k_3}{k_2} \\ Y_0 &=& \frac{k_1}{k_2}A\end{cases}$
- 定常解の周辺に無数の周期的な軌道が存在する
リミットサイクル
- ブリュッセレーター
- $\begin{matrix} \frac{dX}{dt} &=& A + X^2Y -BX -X \\ \frac{dY}{dt} &=& BX-X^2Y \end{matrix}$
- 定常解は $\begin{cases} X_0 &=& A \\ Y_0 &=& \frac{B}{A} \end{cases}$
- 下のソースは一点に収束する
- リミットサイクルの条件にしたのはこちら

A<-2
k1<-0.1
k2<-0.2
k3<-0.2

#これは定常解
X0<-k3/k2
Y0<-k1*A/k2

#定常解からずらす
X<-X0+runif(1)
Y<-Y0+runif(1)

T<-100
dt<-0.005
Nt<-T/dt

Xn<-c()
Yn<-c()

for(t in 1:Nt){
	Xn[t]<-X
	Yn[t]<-Y	
	dX<-(k1*A*X-k2*X*Y)*dt
	dY<-(k2*X*Y-k3*Y)*dt
	X<-X+dX
	Y<-Y+dY
}

plot(Xn,Yn,xlim=c(0,2),ylim=c(0,2))
par(new=TRUE)
plot(X0,Y0,xlim=c(0,2),ylim=c(0,2),col=4)

#ブリュッセレーター
A<-B<-0.5

X0<-A
Y0<-B/A

X<-X0+runif(1)/2
Y<-Y0+runif(1)/2

Xn<-Yn<-c()

T<-10
dt<-0.01
Nt<-T/dt

for(n in 1:Nt){
Xn[n]<-X
Yn[n]<-Y

dX<-A+X^2*Y-B*X-X
dY<-B*X-X^2*Y

X<-X+dX
Y<-Y+dY

}

plot(Xn,Yn,xlim=c(X0-0.5,X0+0.5),ylim=c(Y0-0.5,Y0+0.5))

初期位置をふって時間ごとにプロットしてみる (こちら参照)

library(rgl)

dt<-0.1
T<-30
Nt<-T/dt
M<-matrix(0,3,Nt+1)
N<-c()

A<-2
k1<-0.1
k2<-0.2
k3<-0.2

#これは定常解
X0<-k3/k2
Y0<-k1*A/k2

f1<-function(X){
	y<- k1*A*X[1]-k2*X[1]*X[2]
	return(y)
	}
f2<-function(X){
	y<- k2*X[1]*X[2]-k3*X[2]
	return(y)
	}


#table<-matrix(c(0,0,1,1,2,2,3,3),nrow=4,byrow=TRUE)
table<-data.matrix(expand.grid(-5:5/10+X0,-5:5/10+Y0))

dimnames(table)<-NULL
k<-nrow(table)

for(i in 1:k){

X<-table[i,]
M[,1]<-c(X,0)

for(nt in 1:Nt){
dX<-c(f1(X),f2(X))*dt
X<-X+dX
#print(X)
M[,nt+1]<-c(X,nt)
}
N<-cbind(N,M)
}


range<-c(min(N[1:2,]),max(N[1:2,]))

# 3dプロット（rgl）
# plot3d(N[1,],N[2,],N[3,],xlab="x1",ylab="x2",zlab="time",xlim=range,ylim=range,zlim=c(1,Nt+1)),col=rainbow(Nt+1)[N[3,]+1],type="p")  # 時間ごとに色分け
plot3d(N[1,],N[2,],N[3,],xlab="x1",ylab="x2",zlab="time",xlim=range,ylim=range,zlim=c(1,Nt+1),col=rainbow(k)[sort(rep(1:k,Nt+1))],type="p")   # 点ごとに色分け

# 線でプロット
open3d()
for(j in 1:k){
lines3d(t(N[,1:(Nt+1)+(j-1)*(Nt+1)]),col=rainbow(k)[j])   #点ごとに色分け
# lines3d(t(N[,1:(Nt+1)+(j-1)*(Nt+1)]),col=rainbow(1+Nt)[N[3,]])  #時間ごとに色分け
par3d(scale=c(1/dist(range),1/dist(range),1/(Nt+1)))
}


#2dプロット
for(j in 1:(Nt+1)){	
#	plot(N[1,which(N[3,]==j-1)],N[2,which(N[3,]==j-1)],pch=16,xlab="x1",ylab="x2",xlim=range,ylim=range,col=rainbow(Nt+1)[j])  # 時間ごとに色分け
	plot(N[1,which(N[3,]==j-1)],N[2,which(N[3,]==j-1)],pch=16,xlab="x1",ylab="x2",xlim=range,ylim=range,col=rainbow(k)[1:k])    #点ごとに色分け
	par(new=TRUE)  #重ね書きの有無
	}

#ブリュッセレーターを経時的にプロット
library(rgl)

A<-B<-1
X0<-A
Y0<-B/A

dt<-0.1
T<-20
Nt<-T/dt
M<-matrix(0,3,Nt+1)
#A[,1]<-c(X,0)
N<-c()

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)
	}

table<-data.matrix(expand.grid(0:20/10,0:20/10))
dimnames(table)<-NULL
k<-nrow(table)

for(i in 1:k){

X<-table[i,]
M[,1]<-c(X,0)

for(nt in 1:Nt){
dX<-c(f1(X),f2(X))*dt
X<-X+dX
#print(X)
M[,nt+1]<-c(X,nt)
}
N<-cbind(N,M)
}


range<-c(min(N[1:2,]),max(N[1:2,]))

#3dプロット（rgl）
# plot3d(N[1,],N[2,],N[3,],xlab="x1",ylab="x2",zlab="time",xlim=range,ylim=range,zlim=c(1,Nt),col=rainbow(Nt+1)[N[3,]+1],type="p")  #時間ごとに色分け
plot3d(N[1,],N[2,],N[3,],xlab="x1",ylab="x2",zlab="time",xlim=range,ylim=range,zlim=c(1,Nt),col=rainbow(k)[sort(rep(1:k,Nt))],type="p")   #点ごとに色分け

#2dプロット
for(j in 1:(Nt+1)){	
#	plot(N[1,which(N[3,]==j-1)],N[2,which(N[3,]==j-1)],pch=16,xlab="x1",ylab="x2",xlim=range,ylim=range,col=rainbow(Nt+1)[j])  #時間ごとに色分け
	plot(N[1,which(N[3,]==j-1)],N[2,which(N[3,]==j-1)],pch=16,xlab="x1",ylab="x2",xlim=range,ylim=range,col=rainbow(k)[1:k])   #点ごとに色分け
	par(new=TRUE)  #重ね書きの有無
	}

library(fields)


A<-0.5
B<-0.5

dim<-2

# 連立方程式
f<-function(x){
z<-rep(0,dim)
z[1]<-A+x[1]^2*x[2]-B*x[1]-x[1]
z[2]<-B*x[1]-x[1]^2*x[2]
return(z)
}

vectors<-c()
errors<-c()

# 初期値の集合 
by<-0.025
x0<-seq(from=0,to=2,by=by)
y0<-seq(from=0,to=2,by=by)

points<-data.matrix(expand.grid(x0,y0))

lx<-length(x0)
ly<-length(y0)
L<-nrow(points)

pmat_x<-matrix(points[,1],lx,ly)
pmat_y<-matrix(points[,2],lx,ly)

vectors<-points

for(l in 1:L){
XY<-points[l,]
dXY<-f(XY)
vectors[l,]<-dXY

}

# ベクトルでプロット
dev.new()
#png("vector field.png")
plot(points,xlim=range(points[,1]),ylim=range(points[,2]),type="n")

arrow.plot(points[,1],points[,2],vectors[,1],vectors[,2],arrow.ex=.2, length=.1,col=tim.colors(5)[(atan2(vectors[,2],vectors[,1])+pi)/(2*pi)*4+1], lwd=2,xlim=range(points[1,]),ylim=range(points[2,])) 
# dev.off()



# 境界を求める
vmat_x<-matrix(vectors[,1],lx,ly)
vmat_y<-matrix(vectors[,2],lx,ly)

bound_x<-(vmat_x[-1,]*vmat_x[-lx,]<0)
bound_y<-(vmat_y[,-1]*vmat_y[,-ly]<0)

# xの境界
bx_x<-(pmat_x[-1,][bound_x == 1]+pmat_x[-lx,][bound_x == 1])/2
bx_y<-(pmat_y[-1,][bound_x == 1]+pmat_y[-lx,][bound_x == 1])/2

# yの境界
by_x<-(pmat_x[,-1][bound_y == 1]+pmat_x[,-ly][bound_y == 1])/2
by_y<-(pmat_y[,-1][bound_y == 1]+pmat_y[,-ly][bound_y == 1])/2


xstable<-(vmat_x[-1,]<0)*(vmat_x[-lx,]>0)
xstable_x<-(pmat_x[-1,][xstable == 1]+pmat_x[-lx,][xstable == 1])/2
xstable_y<-(pmat_y[-1,][xstable == 1]+pmat_y[-lx,][xstable == 1])/2

xunstable<-(vmat_x[-1,]>0)*(vmat_x[-lx,]<0)
xunstable_x<-(pmat_x[-1,][xunstable == 1]+pmat_x[-lx,][xunstable == 1])/2
xunstable_y<-(pmat_y[-1,][xunstable == 1]+pmat_y[-lx,][xunstable == 1])/2

ystable<-(vmat_y[,-1]<0)*(vmat_y[,-ly]>0)
ystable_x<-(pmat_x[,-1][ystable == 1]+pmat_x[,-ly][ystable == 1])/2
ystable_y<-(pmat_y[,-1][ystable == 1]+pmat_y[,-ly][ystable == 1])/2

yunstable<-(vmat_y[,-1]>0)*(vmat_y[,-ly]<0)
yunstable_x<-(pmat_x[,-1][yunstable == 1]+pmat_x[,-ly][yunstable == 1])/2
yunstable_y<-(pmat_y[,-1][yunstable == 1]+pmat_y[,-ly][yunstable == 1])/2








# 角度と境界をプロット
dev.new()
# png("angle.png")
par(mfrow=c(2,2))

#1
plot(points,xlim=range(points[,1]),ylim=range(points[,2]),col=tim.colors(501)[(atan2(vectors[,2],vectors[,1])+pi)/(2*pi)*500+1], lwd=2,pch=15,cex=1,main="angle",xlab="X",ylab="Y") 

#2
plot(points,xlim=range(points[,1]),ylim=range(points[,2]),col=tim.colors(501)[(atan2(vectors[,2],vectors[,1])+pi)/(2*pi)*500+1], lwd=2,pch=15,cex=1,main="merge",xlab="X",ylab="Y") 

par(new=T)
plot(xstable_x,xstable_y,pch=16,xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")
par(new=T)
plot(xunstable_x,xunstable_y,pch=3,xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")
par(new=T)
plot(ystable_x,ystable_y,pch=16,col="red",xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")
par(new=T)
plot(yunstable_x,yunstable_y,pch=3,col="red",xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")

#3
plot(xstable_x,xstable_y,pch=16,xlim=range(points[,1]),ylim=range(points[,2]),main="boundary of x",xlab="X",ylab="Y")
par(new=T)
plot(xunstable_x,xunstable_y,pch=3,xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="",)


#4
plot(ystable_x,ystable_y,pch=16,col="red",xlim=range(points[,1]),ylim=range(points[,2]),main="boundary of y",xlab="X",ylab="Y")
par(new=T)
plot(yunstable_x,yunstable_y,pch=3,col="red",xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")

# dev.off()

# 安定性
dev.new()
# png("Stability.png",width=960)
par(mfrow=c(1,2))

#1
plot(points,xlim=range(points[,1]),ylim=range(points[,2]),col=tim.colors(501)[ifelse((vmat_x<0),100,500)], lwd=2,pch=15,cex=1,main="x stability",xlab="X",ylab="Y")

par(new=T)
plot(xstable_x,xstable_y,pch=16,xlim=range(points[,1]),ylim=range(points[,2]),main="",xlab="",ylab="")
par(new=T)
plot(xunstable_x,xunstable_y,pch=3,xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")

#2
plot(points,xlim=range(points[,1]),ylim=range(points[,2]),col=tim.colors(501)[ifelse((vmat_y<0),100,500)], lwd=2,pch=15,cex=1,main="y stability",xlab="X",ylab="Y")

par(new=T)
plot(ystable_x,ystable_y,pch=16,col="red",xlim=range(points[,1]),ylim=range(points[,2]),main="",xlab="",ylab="")
par(new=T)
plot(yunstable_x,yunstable_y,pch=3,col="red",xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")
# dev.off()


dev.new()
plot(points,xlim=range(points[,1]),ylim=range(points[,2]),type="n")

arrow.plot(points[,1],points[,2],vectors[,1],vectors[,2],arrow.ex=.2, length=.1,col=tim.colors(5)[(atan2(vectors[,2],vectors[,1])+pi)/(2*pi)*4+1], lwd=2,xlim=range(points[1,]),ylim=range(points[2,])) 

par(new=T)
plot(xstable_x,xstable_y,pch=16,xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")
par(new=T)
plot(xunstable_x,xunstable_y,pch=3,xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")
par(new=T)
plot(ystable_x,ystable_y,pch=16,col="red",xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")
par(new=T)
plot(yunstable_x,yunstable_y,pch=3,col="red",xlim=range(points[,1]),ylim=range(points[,2]),xlab="",ylab="")