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

• ロトカ-ヴォルテラの微分方程式を以下のように書こう
  • 
  • 
• が固定点になる
• の場合には、２つの固定点のうち、が回転の中心になる
• 一般的にとしている
• 変化のベクトルの長さは
• 曲線のMoving frameの第一ベクトルは
• Moving frameの第二ベクトルは
• 曲率は:内積で
• 
• これをＲで描く
  • 曲率の絶対値でプロットすると、「曲がり具合が同じところ」に等高線が引かれる
  • 円というのは、「曲率の等高線『も』円を描く場」にできる軌道（？）

# ２つの固定点
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)))