曲線を円の変形と考える

  • こちらで、2次元正円を2次元周回曲線に対応付ける関数について考えた
  • 対応関数は解けないことも多いので、RのrootSolveパッケージを使うことも書いた
  • ロトカ=ヴォルテラのx,y対称な曲線を描いてみる(x-\log(x)+1 + y -\log(y)+1 = R^2
  • うまくいっているように見える(こちらで紹介したロトカ=ヴォルテラの曲線描図プログラムを少し改変して、始点を指定できるようにすることで、2方法の曲線が重なることを示した)
  • 別の関数(x+\frac{1}{x}-2)+(y+\frac{1}{y}-2) = R^2でも描ける。この関数ではx=1で値が0でありそれは最小値。
  • さらに、ロトカ=ヴォルテラと同様に0 < x,y \le \inftyが定義域
library(rootSolve)
t<-seq(from=0,to=1,by=0.01)*2*pi
# 正円のパラメタ使用の座標
x1<-cos(t)
y1<-sin(t)

# 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))
plot(c(minx2,minx2,maxx2,maxx2),c(miny2,maxy2,miny2,maxy2),cex=1,col="light blue",xlim=xlim,ylim=ylim)

xy<-lotkaVolterra2(a=1,0.001,minx2[10],miny2[10],10000)
par(new=TRUE)
plot(xy$x,xy$y,xlim=xlim,ylim=ylim,cex=0.1)