- 昨日のこの記事で、定常起動のドリフトのことを書いた
- 定常起動は「保存量が保存された状態」だから、そのドリフトは「保存量のドリフト」
- 保存量は正円座標系では半径に相当するからそれをドリフトさせてみる
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
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
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,]))
}
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,]))
}
minx2<-apply(x2,1,min)
maxx2<-apply(x2,1,max)
miny2<-apply(y2,1,min)
maxy2<-apply(y2,1,max)
lotkaVolterra2 <- function (a, dt, x0,y0, n) {
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
}
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")