ブリュッセレーター

R シミュレーション グラフィクス

KABIRA2011-07-21

  • こちら のつづき
  • リミットサイクル
    • リミットサイクルになる条件
    • B > 1 + A^2

library(rgl)

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

dt<-0.1
T<-30
Nt<-T/dt
M<-matrix(0,3,Nt+1)
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<-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(fields)


A<-0.5
B<-2
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=5,by=by)
y0<-seq(from=0,to=5,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



# どっち？
# vmat or vmat * pmat ?? 

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="")