- 種の捕食・被捕食関係がロトカ=ヴォルテラ。昨日は、2種の捕食・被捕食関係だった
- N種を考えよう
- 単純な系として、X1->X2->...->XN->X1という循環型捕食関係とする
- 考え方は昨日の3種と同じ
n<-5
us<-sample(1:10,n)
Nr<-n
ks<-runif(Nr)
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
Niter<-100000
dt<-0.01
m<-matrix(0,Niter,n)
p<-rep(0,Niter)
m[1,]<-runif(n)
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]
}
}
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)