保存面の交線 - ryamadaの遺伝学・遺伝統計学メモ

こちらで3種の捕食・被捕食関係を見た
そこでの前提は、環境全体で「生物構成要素保存」の法則を成り立たせていた
- $\sum k_i X_i$ が一定
- これは、生物種の数Nについて、N次元空間にある、 $\sum k_i X_i=C$ なる平面上を $X=(X_i)$ は動くことを意味する
これは $N=3$ 種の場合には、平面である
今、以下のような微分方程式が成り立つとする
- $\frac{dx}{dt}=-a_x xy + b_x zx$
- $\frac{dy}{dt}=-a_y yz + b_y xy$
- $\frac{dz}{dt}=-a_z zx + b_z yz$
$a_xa_ya_z=b_xb_yb_z$ が原点以外の固定点( $\frac{dx}{dt}=\frac{dy}{dt}=\frac{dz}{dt}=0$ )を持つ条件であることは式変形から確かめられる
今、保存量 $E$ があるとすれば、 $\frac{dE}{dt}=0$ である
$\frac{dE}{dt}=\frac{dE}{dx}\frac{dx}{dt}+\frac{dE}{dy}\frac{dy}{dt}+\frac{dE}{dz}\frac{dz}{dt}=0$ である
の対称性に注意すれば
- $\frac{dE}{dx}=\frac{1}{x}$
- $\frac{dE}{dy}=\frac{a_x}{b_y}\frac{1}{y}$
- $\frac{dE}{dz}=\frac{b_x}{a_z}\frac{1}{z}$
のときに $\frac{dE}{dt}=0$ が満足されるから
$E=\log(x) + \frac{a_x}{b_y}\log(y) +\frac{b_x}{a_z}\log(z)$ が保存量であることがわかる
以下のソースはそれを確かめるソースであり、３種ロトカ=ヴォルテラの軌道は、この保存量と、 $\sum k_i X_i=C$ 平面との交線となる

library(lattice)
as<-runif(3)
bs<-runif(3)
bs[3]<-cumprod(as)[3]/(bs[1]*bs[2])
cumprod(c(as,1/bs))
dt<-0.01
Niter<-1000
da<-matrix(0,Niter,3)
da[1,]<-runif(3)

p<-q<-rep(0,Niter)
q[1]<-log(da[1,1])+as[1]/bs[2]*log(da[1,2])+bs[1]/as[3]*log(da[1,3])

for(i in 2:Niter){
	xyz<-cumprod(da[i-1,])[3]
	da[i,1]<-da[i-1,1]+xyz*(-as[1]/da[i-1,2]+bs[1]/da[i-1,3])*dt
	da[i,2]<-da[i-1,2]+xyz*(-as[2]/da[i-1,3]+bs[2]/da[i-1,1])*dt
	da[i,3]<-da[i-1,3]+xyz*(-as[3]/da[i-1,1]+bs[3]/da[i-1,2])*dt
	q[i]<-log(da[i,1])+as[1]/bs[2]*log(da[i,2])+bs[1]/as[3]*log(da[i,3])

}
plot(as.data.frame(da))

plot(q)

plot3d(da[,1],da[,2],da[,3])


mmmm<-1/100
x<-y<-seq(from=mmmm,to=1,length.out=100)
xy<-expand.grid(x,y)
z<-z2<-rep(0,length(xy[,1]))

for(i in 1:length(xy[,1])){
	z[i]<-exp(as[3]/bs[1]*(-log(xy[i,1])-as[1]/bs[2]*log(xy[i,2])))
}

for(i in 1:length(xy[,1])){
	z2[i]<-1/(bs[2]*bs[3])*as[2]*(-xy[i,1]*as[1]-xy[i,2]*bs[2])

}
C<-min(z-z2)+1
z2<-z2+C

xx<-rep(xy[,1],2)
yy<-rep(xy[,2],2)
zz<-c(z,z2)


gr<-c(rep(1,length(xy[,1])),rep(2,length(xy[,1])))

wireframe(zz ~ xx * yy, groups = gr,screen = list(x = 0, y =0,z=0))
wireframe(zz ~ xx * yy, groups = gr,screen = list(x = 80, y =0,z=10))