固定点 - ryamadaの遺伝学・遺伝統計学メモ

書きかけ
２変量のロトカ=ヴォルテラでは原点以外に１個の固定点があって、軌道はその周囲にぐるりとできる
３変量に拡張（３変量が共通の構成要素で制約し合っている状態）では、構成要素保存が曲面を作り、そこに１つの固定点（原点以外のそれ）があり、その周囲に周回曲線軌道ができる。固定点は、原点を通る直線をなすが、「構成要素保存曲面」を１点のみで貫くので、「構成要素保存曲面」にはただ一つの固定点がある
４変量に拡張すると、「構成要素保存」と「偏微分連立方程式」が定めるもう一つの「保存式」とは、３次元多様体を作る。それを、原点を通る１直線が貫き、それが、「固定点の連なりとしての固定直線」。
この固定直線を周回軌道が回ることで、結果として、３次元多様体が、「構成要素保存」と「偏微分連立方程式由来保存量の保存」が軌道によって埋め尽くされる
ある点は、ぞの座標とそのときの速度が一意に決まり、どの「軌道」に乗るかは確定するし、ある点を通る軌道はただの１本

n<-4 # 種の数
us<-sample(1:10,n) # それぞれの種の単位Pあたり個体数
# 捕食・被捕食関係の数Nrは種の数と同じ
Nr<-n
# 捕食・被捕食効率
ks<-runif(Nr)
# Inは被捕食種。１度の出会いで殺される個体数
# Playerは捕食・被捕食関係で参加する種とその個体数
# 今は１捕食個体と１被捕食個体の出会いとする
# Outは出会いの結果生じる種とその個体数
# 今は捕食者
In<-Player<-diag(n)
Out<-matrix(0,Nr,n)
for(i in 1:Nr){
	for(j in 1:n){
		if(j==i+1){
			Player[i,j]<-1
			Out[i,j]<-1
		}
	}
}
Player[Nr,1]<-Out[Nr,1]<-1
Nrepeat<-20
stockM<-NULL
tmpinit<-runif(n)
for(rr in 1:Nrepeat){
# ここから時系列シミュレーション
Niter<-1000
dt<-0.01
m<-matrix(0,Niter,n)
p<-rep(0,Niter)
if(rr==1){
	m[1,]<-tmpinit

}else{
	m[1,]<-tmpinit+0.01*(rr-1)
}

p[1]<-0
for(i in 1:n){
	p[1]<-p[1]+us[i]*m[1,i]
}
m[1,]<-m[1,]/p[1]
p[1]<-0
for(i in 1:n){
	p[1]<-p[1]+us[i]*m[1,i]
}

for(i in 2:Niter){
	m[i,]<-m[i-1,]
	for(j in 1:Nr){
		tmp<-ks[j]*dt
		pamountIn<-0
		pamountOut<-0
		for(k in 1:n){
			tmp<-tmp*m[i-1,k]^Player[j,k]
			pamountIn<-pamountIn+us[k]*In[j,k]
			pamountOut<-pamountOut+us[k]*Out[j,k]
		}

		for(k in 1:n){
			m[i,k]<-m[i,k]-In[j,k]*tmp+tmp*pamountIn/pamountOut*us[k]*Out[j,k]/pamountOut*Out[j,k]
			#m[i,k]<-m[i,k]-In[j,k]*tmp+Out[j*k]*tmp*pamountIn/pamountOut
		}
	}

	for(j in 1:n){
		p[i]<-p[i]+us[j]*m[i,j]
	}

}
plot(p,type="l",ylim=c(0,2),main="保存量")
plot(data.frame(m),cex=0.1)

#xlim<-ylim<-c(min(m),max(m))
#plot(data.frame(m),cex=0.1,xlim=xlim,ylim=ylim)
if(rr==1){
	stockM<-stockM
}else{
	stockM<-rbind(stockM,m)

}
}

plot(data.frame(stockM),cex=0.1)

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