ロトカ-ヴォルテラの曲率

  • ロトカ-ヴォルテラの微分方程式を以下のように書こう
    • \frac{dx}{dt}=a(x-x_1)(y-y_2)
    • \frac{dy}{dt}=b(x-x_2)(y-y_1)
  • (x_1,y_1),(x_2,y_2)が固定点になる
  • a>0,b<0の場合には、2つの固定点のうち、(x_2,y_2)が回転の中心になる
  • 一般的に(x_1,y_1)=(0,0)としている
  • 変化のベクトルv=(v_x,v_y)=(\frac{dx}{dt},\frac{dy}{dt})の長さは|v|=\sqrt{v_x^2+v_y^2}
  • 曲線のMoving frameの第一ベクトルはe_1=\frac{1}{|v|}(v_x,v_y)
  • Moving frameの第二ベクトルはe_2=\frac{1}{|v|}(-v_y,v_x)
  • 曲率はk=(\frac{d e_1}{ds},e_2):内積
  • ab(x-x_1)(y-y_1)(b(x-x_2)^2(y_1-y_2)-a(y-y_2)^2(x_1-x_2))
  • これをRで描く
    • 曲率の絶対値でプロットすると、「曲がり具合が同じところ」に等高線が引かれる
    • 円というのは、「曲率の等高線『も』円を描く場」にできる軌道(?)
# 2つの固定点
xy1<-c(0,0)
xy2<-c(0.4,0.4)
a<-2
b<--2
n<-100
minX<-0
maxX<-1
X<-Y<-seq(from=minX,to=maxX,length.out=n)
XY<-expand.grid(X,Y)
Z<-rep(0,length(XY[,1]))
for(i in 1:length(Z)){
	tmpx<-XY[i,1]
	tmpy<-XY[i,2]
	v<-c(a*(tmpx-xy1[1])*(tmpy-xy2[2]),b*(tmpx-xy2[1])*(tmpy-xy1[2]))
	Z[i]<-a*b*(tmpx-xy1[1])*(tmpy-xy1[2])*(b*(tmpx-xy2[1])^2*(xy1[2]-xy2[2])-a*(tmpy-xy2[2])^2*(xy1[1]-xy2[1]))/(sum(v^2)^(3/2))
	#Z[i]<-a*b*(tmpx-xy1[1])*(tmpy-xy1[2])*(b*(tmpx-xy2[1])^2*(xy1[2]-xy2[2])-a*(tmpy-xy2[2])^2*(xy1[1]-xy2[1]))
}
persp(X,Y,matrix(Z,length(X),length(Y)))
library(rgl)
plot3d(XY[,1],XY[,2],Z)
plot3d(XY[,1],XY[,2],Z,col=rainbow(1000))
image(X,Y,matrix(log(abs(Z)),length(X),length(Y)))