力 - ryamadaの遺伝学・遺伝統計学メモ

位置、その微分の速度、その微分の力、についてこちらで書いた
惑星の軌道が楕円（や放物線）を描くのは
$\frac{d^2x_i}{dt^2}=F_{xi}=-k\times x_i$ というように、力の場ができているから
この場は原点でゼロで、それ以外では位置座標に比例して大きく、向きは原点方向、そして、すべての軸で係数 $k$ と等しい（→これがすべての軸で等周期になる理由）
これと同じことをロトカ=ヴォルテラの偏微分方程式に当てはめれば
$\frac{d^2x}{dt^2}=(x-x_0)(\beta^2(y-y_1)^2-\beta\gamma(y-y_0)(x-x_1))$
- ただし $x_0=y_0=0$
これは、 $x=x_0,y=y_0$ と $x=x_1,y=y_1$ で速度変化がなく、速度もゼロ
$x=x_0,y\not = y_0$ ではy方向に速度は増加、 $x\not = x_0,y=0$ でもx方向に速度が増加
$F_x,F_y$ の正負の変化の様子を描くと掲載図の通り

# ロトカ
k<-sort(runif(4))
k<-rep(1,4)
X<-Y<-seq(from=0,to=1,length.out=100)*3
XY<-expand.grid(X,Y)
Fx<-Fy<-rep(0,length(XY[,1]))
for(i in 1:length(Fx)){
	x<-XY[i,1]
	y<-XY[i,2]
	Fx[i]<-x*(k[2]^2*(y-k[1]/k[2])^2-k[2]*k[3]*y*(x-k[4]/k[3]))
	Fy[i]<-y*(k[3]^2*(x-k[4]/k[3])^2-k[2]*k[3]*x*(y-k[1]/k[2]))
}
Fxy<-sqrt(Fx^2+Fy^2)
Fpm<-sign(Fx)*sign(Fy)

par(mfcol=c(2,2))
image(X,Y,matrix(Fx,length(X),length(Y)),col=gray((1:20)/20))
points(k[4]/k[3],k[1]/k[2])

image(X,Y,matrix(Fy,length(X),length(Y)),col=gray((1:20)/20))
points(k[4]/k[3],k[1]/k[2])
image(X,Y,matrix(Fxy,length(X),length(Y)),col=gray((1:20)/20))
points(k[4]/k[3],k[1]/k[2],col=2)

par(mfcol=c(1,1))
image(X,Y,matrix(Fpm,length(X),length(Y)),col=gray((1:20)/20))
points(k[4]/k[3],k[1]/k[2],col=2)
points(k[4]/k[3],k[1]/k[2],col=2,cex=4)
# 万有引力
k<-runif(4)
X<-Y<-seq(from=-1,to=1,length.out=100)*10
XY<-expand.grid(X,Y)
Fx<-Fy<-rep(0,length(XY[,1]))
for(i in 1:length(Fx)){
	x<-XY[i,1]
	y<-XY[i,2]
	Fx[i]<-x
	Fy[i]<-y
}

Fxy<-sqrt(Fx^2+Fy^2)
Fxy<-sign(Fx)*sign(Fy)

par(mfcol=c(2,2))
image(X,Y,matrix(Fx,length(X),length(Y)),col=gray((1:20)/20))
points(k[4]/k[3],k[1]/k[2])

image(X,Y,matrix(Fy,length(X),length(Y)),col=gray((1:20)/20))
points(k[4]/k[3],k[1]/k[2])
image(X,Y,matrix(Fxy,length(X),length(Y)),col=gray((1:20)/20))
points(k[4]/k[3],k[1]/k[2])

par(mfcol=c(1,1))

初期位置と初期速度を与えて、「力」を与えてもロトカ=ヴォルテラの周回曲線は描ける
初期位置と初期速度が不適当だと、すぐに周回しなくなる

# ロトカ
ns<-2
k<-rnorm(4)
#k<-rep(1,4)
X<-Y<-seq(from=0,to=1,length.out=100)*3
XY<-expand.grid(X,Y)
Fx<-Fy<-rep(0,length(XY[,1]))
Niter<-10000
dt<-1/1000

Xs<-matrix(0,Niter,ns)
deviation<-0.6
Xs[1,]<-c(k[4]/k[3]*deviation,k[1]/k[2])
Vs<-c(k[1]*Xs[1,1]-k[2]*Xs[1,1]*Xs[1,2],k[3]*Xs[1,1]*Xs[1,2]-k[4]*Xs[1,2])
for(i in 2:Niter){
	Xs[i,]<-Xs[i-1,]+Vs*dt
	x<-Xs[i-1,1]
	y<-Xs[i-1,2]
	Fs<-rep(0,2)
	Fs[1]<-x*(k[2]^2*(y-k[1]/k[2])^2-k[2]*k[3]*y*(x-k[4]/k[3]))
	Fs[2]<-y*(k[3]^2*(x-k[4]/k[3])^2-k[2]*k[3]*x*(y-k[1]/k[2]))
	Vs<-Vs+Fs*dt
}

plot(Xs,type="l")
matplot(Xs,type="l")