分岐

• 存在から発展へ
  • p120
  • パラメータBの値をかえたときの挙動
• 定常状態になるまで時間がかかる
  • ほど計算させる必要がある
• 領域の大きさの影響もないとは言い切れない
  • 以下のソースは図 5.10
    • の領域内で計算するもの
    • 本では での計算
    • （）
  • 図 5.11 のような回転解がみれていない
• マスを増やして計算していると次のようなエラーがでた
  • 解読中...

R(1270,0xa020a4e0) malloc:

mmap(size=2441216) failed (error code=12)

error: can't allocate region

set a breakpoint in malloc_error_break to debug

library(rgl)

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

Lx<-1.18
dx<-0.02
Nx<-floor(Lx/dx)

Ly<-1.18
dy<-0.02
Ny<-floor(Ly/dy)

T<-1
dt<-0.001
Nt<-T/dt

# 拡散係数
Dx<-0.008
Dy<-0.004

# 補正
Dx<-Dx/dx^2*dt
Dy<-Dy/dy^2*dt


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

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


# 各点のX,Yの濃度
C<-array(c(X0),c(Nx,Ny,2))

# C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),2]<-Y0+seq(from=0,to=0.1,length=4)
C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),1]<-X0
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


# ここからがシミュレーション	
for(t in 1:Nt)
{		
	if(t %% 100 ==1){
	M<-max(C[,,1])
	m<-min(C[,,1][which(C[,,1]!=0)])
		if(M == m){
		Col<-1
		}else{
		Col<-(C[,,1]-m)/(M-m)*100*(K>0)+1
		}
	plot3d(slice.index(C[,,1],1)*Lx/Nx,slice.index(C[,,1],2)*Ly/Ny,C[,,1],alpha=ifelse(K==0,0,1),xlab="x",ylab="y",zlab="X conc.",col=rainbow(101)[Col],main=paste("time = ",(t-1)*dt,sep=""))
#	rgl.snapshot(paste(1000+(t-1)))
	}	
	
	# ブリュッセレータの計算
	tempC<-(K>0)*(A+C[,,1]^2*C[,,2]-B*C[,,1]-C[,,1])*dt
	C[,,2]<-C[,,2]+(K>0)*(B*C[,,1]-C[,,1]^2*C[,,2])*dt
	C[,,1]<-C[,,1]+tempC
	
	
	
	# 拡散の計算
	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
}

#色で表示するだけにする
library(rgl)

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

Lx<-1.18
dx<-0.02
Nx<-floor(Lx/dx)

Ly<-1.18
dy<-0.02
Ny<-floor(Ly/dy)

T<-1
dt<-0.001
Nt<-T/dt

# 拡散係数
Dx<-0.008
Dy<-0.004

# 補正
Dx<-Dx/dx^2*dt
Dy<-Dy/dy^2*dt


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

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


# 各点のX,Yの濃度
C<-array(c(X0),c(Nx,Ny,2))

# C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),2]<-Y0+seq(from=0,to=0.1,length=4)
C[floor((Nx+1)/2):floor((Nx)/2+1),floor((Ny+1)/2):floor((Ny)/2+1),1]<-X0
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


# ここからがシミュレーション	
for(t in 1:Nt)
{		
	if(t %% 100 ==1){
	M<-max(C[,,1])
	m<-min(C[,,1][which(C[,,1]!=0)])
		if(M == m){
		Col<-1
		}else{
		Col<-(C[,,1]-m)/(M-m)*100*(K>0)+1
		}
	plot3d(slice.index(K,1)*Lx/Nx,slice.index(K,2)*Ly/Ny,0,alpha=ifelse(K==0,0,1),xlab="x",ylab="y",zlab="X conc.",col=rainbow(101)[Col],main=paste("time = ",(t-1)*dt,sep=""))
	rgl.viewpoint(theta=0,phi=0)
#	rgl.snapshot(paste(1000+(t-1)))
	}	
	
	# ブリュッセレータの計算
	tempC<-(K>0)*(A+C[,,1]^2*C[,,2]-B*C[,,1]-C[,,1])*dt
	C[,,2]<-C[,,2]+(K>0)*(B*C[,,1]-C[,,1]^2*C[,,2])*dt
	C[,,1]<-C[,,1]+tempC
	
	
	
	# 拡散の計算
	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
}