ロトカ=ヴォルテラの方程式とその保存量を拡張してみる

ロトカ=ヴォルテラの方程式は捕食者・被捕食者の個体数が周期的変動をする様子を表した非線形偏微分方程式(Wiki記事)
被捕食者(x)は、その個体数に比例して増殖するが、捕食者と出会う確率に応じて減少する
捕食者(y)は、その個体数に比例して自然死するが、被捕食者と出会う確率に応じて増殖する
$\frac{dx}{dt} = \alpha x -\beta xy$
$\frac{dy}{dt} = \gamma xy -\kappa y$
次のx,yの関数Vは時間微分が0であることが示せる
$V(t)=\gamma x +\beta y - \kappa \log(x) -\alpha \log(y)$
$\frac{dV(t)}{dt}=\gamma \frac{dx}{dt}+\beta \frac{dy}{dt} -\frac{\kappa}{x} \frac{dx}{dt} -\frac{\alpha}{y} \frac{dy}{dt}$ を素直に式変形すれば示せる
時間微分が０なので $V(t)$ は時間不変量であって、「保存量」と呼ばれる
保存量の関数はどうやったら得られるのかを考えてみる
$V(t)=f(x(t))-g(y(t))$ であるとき
$\frac{V(t)}{dt} = \frac{f(x(t))}{dt}-\frac{g(y(t))}{dt}=\frac{f(x)}{dx}\frac{dx(t)}{dt}-\frac{g(y)}{dy}\frac{dy(t)}{dt}$ であることを利用する
今、時間に関して保存される量の関数 $V(t)$ は $\frac{V(t)}{dt}=0$ である
$\frac{dx}{dt}\frac{dy}{dt} - \frac{dy}{dt}\fra{dx}{dt}=0$ は自明であるから、 $\frac{V(t)}{dt}\propto \frac{dx}{dt}\frac{dy}{dt} - \frac{dy}{dt}\fra{dx}{dt}$ となるような $V(t)$ は保存量となるだろう
$\frac{V(t)}{dt} =\frac{f(x)}{dx}\frac{dx(t)}{dt}-\frac{g(y)}{dy}\frac{dy(t)}{dt}$ と $\frac{V(t)}{dt}\propto \frac{dx}{dt}\frac{dy}{dt} - \frac{dy}{dt}\fra{dx}{dt}$ とを比べると
$\frac{f(x)}{dx} \propto \frac{dy}{dt},\frac{g(y)}{dy} \propto \frac{dx}{dt}$ であって、 $\frac{f(x)}{dx},\frac{g(y)}{dy}$ はそれぞれx,yのみの関数である
ここで、x,yのみの関数 $u(y),v(x)$ を用いて $\frac{dx}{dt}=x^{k_x}\times u(y), \frac{dy}{dt}=y^{k_y}\times v(x)$ と表せれば、
$x \not = 0, y \not = 0$ のとき
$\frac{f(x)}{dx}=\frac{dy}{dt}/(x^{k_x}y^{k_y})=\frac{v(x)}{x^{k_x}}$
$\frac{g(y)}{dy}=\frac{dx}{dt}/(x^{k_x}y^{k_y})=\frac{u(y)}{y^{k_y}$
とできて、保存量が確保できる
保存量があるということは、「突拍子もない状態」には行かないで済む(かもしれない)と言うことなので、周期的になってくれそう
$\frac{dx}{dt}=\frac{dy}{dt}=0$ を満足する点は固定点と呼ばれるが、そこの近傍を初期値に取れば、周期的になる(なりやすい)

# dxdt=ax^k1-bx^k1y
# dydt=cxy^k2-dy^k2
#> a
#[1] 0.3796097
#> b
#[1] 0.8440414
#> c
#[1] 0.5359939
#> d
#[1] 0.4014024
k1<-2
k2<-3
a<-runif(1)
b<-runif(1)
c<-runif(1)
d<-runif(1)
#x<-c(runif(1))
#y<-c(runif(1))
# 固定点(x=d/c,y=a/b)の近傍に初期値をとる
x<-c(d/c*0.9)
y<-c(a/b*0.9)

v<-c(-a/y[1]-b*log(y[1])-(c*x[1]-d*log(x[1])))
Niter<-10000
dt<-100/Niter
for(i in 2:Niter){
	dxdt<-a*x[i-1]^k1-b*x[i-1]^k1*y[i-1]
	dydt<-c*x[i-1]*y[i-1]^k2-d*y[i-1]^k2
	x<-c(x,x[i-1]+dxdt*dt)
	y<-c(y,y[i-1]+dydt*dt)
	v<-c(v,-a/y[i]-b*log(y[i])-(c*x[i]-d*log(x[i])))
}
par(mfcol=c(1,2))
plot(x,y,type="l")

plot(v,ylim=c(-max(abs(v)),max(abs(v))),type="l")
par(mfcol=c(1,1))

par(mfcol=c(1,2))
plot(x,y,type="l")
matplot(cbind(x,y),type="l")

#plot(v,ylim=c(-max(abs(v)),max(abs(v))),type="l")
par(mfcol=c(1,1))

たくさんの固定点を持つ場合