平衡点を求める
- 状態の変化を表した式から平衡状態を求めたい
- 1次元のブリュッセレータで試してみる
- パッケージrootSolveの中のmultiroot関数を使う
- multirootの中身でもRの逆行列（solve）をつかっているようだ
- 黒がX、赤がYの濃度
- 解が安定かどうかはたぶんまた別のはなし
- 解の一意性も別のはなし
- 初期条件をいろいろ変えてみる必要がありそう
- 計算の時間の刻み幅（dt）の値で結果が変わってしまうので拡散係数の修正が必要
- こちらのような数値計算で方程式の解をだせるようにしたい
library(rootSolve) # ブリュッセレータのパラメータ A<-2 B<-4.17 X0<-c(A,B/A) # 平衡となるXとYの値 # シミュレーションのパラメータ L<-1 dr<-0.01 Nr<-floor(L/dr)+1 dt<-0.00001 # 拡散係数 Dx<-0.0016 Dy<-0.008 # 各点のX,Yの濃度、初期条件 C<-matrix(X0,nr=2,nc=Nr) Nint<-2 C<-C+c(sin(Nint*pi*1:Nr/(L*Nr)),sin(Nint*pi*1:Nr/(L*Nr))) # 初期条件 C[,1]<-C[,Nr]<-X0 #固定端 par(mfrow=c(1,2)) plot(seq(from=0,to=1,length=Nr),C[1,],type="l",xlab="r",ylab="X conc.",ylim=range(C),main="initial concentration") par(new=T) plot(seq(from=0,to=1,length=Nr),C[2,],type="l",xlab="r",ylab="Y conc.",ylim=range(C),main="initial concentration",col=2) # ブリュッセレータ f<-function(X){ Y<-c(0,0) Y[1] <- A+X[1]^2*X[2]-B*X[1]-X[1] Y[2] <- B*X[1]-X[1]^2*X[2] return(Y) } # ブリュッセレータと拡散 h<-function(C){ C<-matrix(C,nr=2,nc=Nr) C0<-C # ブリュッセレータの計算 for(l in 1:Nr) { X<-C[,l] dX<-f(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]<-X0 return(C-C0) # 変化量をかえす } # パッケージrootSolveの中のmultirootを使う root<-multiroot(h,C,maxiter=100)$root # print(matrix(h(root),nr=2,nc=Nr)) CC<-matrix(root,nr=2,nc=Nr) plot(seq(from=0,to=1,length=Nr),CC[1,],type="l",xlab="r",ylab="X conc.",ylim=range(CC),main=paste(" equilibirium")) par(new=T) plot(seq(from=0,to=1,length=Nr),CC[2,],type="l",xlab="r",ylab="Y conc.",ylim=range(CC),main=paste(" equilibirium"),col=2) # plot(seq(from=0,to=1,length=Nr),h(CC)[1,],type="l",ylim=c(-0.1,0.1),xlab="r",ylab="X conc.",main=paste("flux"))
- 下はシミュレーション用のコード
A<-2 B<-4.17 X0<-c(A,B/A) # 平衡となるXとYの値 L<-1 dr<-0.01 Nr<-floor(L/dr)+1 T<-3 dt<-0.00001 Nt<-T/dt # 拡散係数 Dx<-0.0016 Dy<-0.008 # 各点のX,Yの濃度 C<-matrix(X0,nrow=2,ncol=Nr) Nint<-2 C<-C+c(sin(Nint*pi*1:Nr/(L*Nr)),sin(Nint*pi*1:Nr/(L*Nr))) range<-range(C) plot(seq(from=0,to=1,length=Nr),C[1,],type="l",ylim=range,xlab="r",ylab="X conc.",main="time = 0") par(new=T) plot(seq(from=0,to=1,length=Nr),C[2,],type="l",ylim=range,xlab="r",ylab="Y conc.",main="time = 0",col=2) # ブリュッセレータ f<-function(X){ Y<-c(0,0) Y[1] <- A+X[1]^2*X[2]-B*X[1]-X[1] Y[2] <- B*X[1]-X[1]^2*X[2] return(Y) } for(t in 1:Nt) { # ブリュッセレータの計算 for(l in 1:Nr) { X<-C[,l] dX<-f(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]<-X0 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="")) par(new=T) 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=""),col=2) } }
library(fields) # 連立方程式 f<-function(x){ z<-rep(0,dim) z[1]<-(x[1]-3)^2+x[2]^2-3 z[2]<-sin(x[1])+exp(x[2]-1)-1 return(z) } # 偏微分係数を近似計算するときに差分を与える関数 g<-function(i){ v<-rep(0,dim) v[i]<-1 return(v) } d<-0.001 # 偏微分係数の近似のための差分 dim<-2 A<-matrix(0,nr=dim,nc=dim) #ヤコビアン vectors<-c() errors<-c() # 初期値の集合 by<-0.05 x0<-seq(from=0,to=5,by=by) y0<-seq(from=-2.5,to=2.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){ x<-XY<-points[l,] # ヤコビアンの計算 for(i in 1:dim){ for(j in 1:dim){ A[i,j]<-(f(x+d*g(j))[i]-f(x-d*g(j))[i])/(2*d) }} dXY<- try(-solve(A)%*%f(x)) if(is(dXY)[1] != "try-error"){# 連立方程式が解けた時 vectors[l,]<-dXY }else{ # 連立方程式が解けなかった時 errors<-cbind(errors,XY) vectors[l,]<-c(0,0) # 0 を入れておく } } # ベクトルでプロット #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()
- こちらに追加
