保存量が浮動する

  • 昨日のこの記事で、定常起動のドリフトのことを書いた
  • 定常起動は「保存量が保存された状態」だから、そのドリフトは「保存量のドリフト」
  • 保存量は正円座標系では半径に相当するからそれをドリフトさせてみる
library(rootSolve)
t<-seq(from=0,to=1,by=0.01)*2*pi
# 正円のパラメタ使用の座標
x1<-cos(t)
y1<-sin(t)
# 保存量にドリフトを入れる
Diam<-rep(0,length(t))
Diam[1]<-1
for(i in 2:length(Diam)){
	Diam[i]<-Diam[i-1]*(1+0.01*rnorm(1))
}
x1<-x1*Diam
y1<-y1*Diam
# x^2 + y^2 = 1
# x-log(x)-1 + y-log(y)-1 = 1

# 座標変換後の値格納行列
x2<-matrix(0,length(x1),2)
y2<-matrix(0,length(x1),2)
# 解探索用関数
funx <- function (x) x-log(x)-1-R^2/(2*sqrt(2))
funy <- funx
#funy <- function (x) x+1/x^2 -3/2^(2/3)- (2*R)^2
#x1<-seq(from=0,to=3,by=0.01)
#x2<-matrix(0,length(x1),2)
# uniroot()を使って2つの対応値を計算
for(i in 1:length(x1)){
	R<-x1[i]
	All <- uniroot.all(funx,c(0,100),maxiter=10000)
	if(length(All)>1)x2[i,1:length(All)]<-All
	if(length(All)==1)x2[i,]<-rep(All,length(x2[1,]))
	#curve(fun(x),0,2*max(All),main="uniroot.all")
	#points(All,y=rep(0,length(All)),pch=16,cex=2)
}

for(i in 1:length(x1)){
	R<-y1[i]
	All <- uniroot.all(funy,c(0,100),maxiter=10000)
	if(length(All)>1)y2[i,1:length(All)]<-All
	if(length(All)==1)y2[i,]<-rep(All,length(y2[i,]))
	#curve(fun(x),0,2*max(All),main="uniroot.all")
	#points(All,y=rep(0,length(All)),pch=16,cex=2)
}
minx2<-apply(x2,1,min)
maxx2<-apply(x2,1,max)
miny2<-apply(y2,1,min)
maxy2<-apply(y2,1,max)
# (x0,y0)が周期の始点座標
lotkaVolterra2 <- function (a, dt, x0,y0, n) { 
               ### 1パラメータバージョン
               ### x,y,t をそれぞれ c/d,b/a,a 倍する、d/a を新たなaとする
 x <- rep(0,n)
 y <- rep(0,n)
 x[1] <- x0    # 不動点からちょっとずらす
 y[1] <- y0     # 不動点
 for (ii in 2:n) {
   x[ii] = x[ii-1] + x[ii-1]*(1-y[ii-1])*dt   # 離散化
   y[ii] = y[ii-1] + a*y[ii-1]*(x[ii-1]-1)*dt # 離散化
 }
 #x11(8,4)
 #matplot(cbind(x,y),type="l")   # 時系列表示
 #x11(5,5)
 #plot(x,y,pch=".")
 list(x=x,y=y)
}

xlim<-ylim<-c(min(x2,y2),max(x2,y2))
xy<-lotkaVolterra2(a=1,0.001,minx2[10],miny2[10],10000)

plot(xy$x,xy$y,xlim=xlim,ylim=ylim,cex=1,col="light blue")
par(new=TRUE)
plot(c(minx2,minx2,maxx2,maxx2),c(miny2,maxy2,miny2,maxy2),cex=0.1,xlim=xlim,ylim=ylim,type="l")