固定点を指定する - ryamadaの遺伝学・遺伝統計学メモ

こちらの記事で、複数の固定点を持つ、ロトカ=ヴォルテラを描いたが、格子状に固定点を定めていた
実際には、適当な数の固定点の位置を任意に決めたい
今、 $V_i=(x_i,y_i);i=1,2,...$ なる点を固定点にするためには、 $\frac{dx}{dt},\frac{dy}{dt}$ の両方が $V_i$ にてゼロとなり、それ以外のところではゼロにならないようにしたい
$\prod_{i=1} ((x-x_i)*f + (y-y_i)*g)$ となるような、非ゼロな $f,g$ を用いて、そのようにできる
この固定点でゼロ、それ以外で正の値を持たせ、それを、その点での速度ベクトルのノルムとする
これに微分可能なように向きをつけてやることにする
各点の速度を固定点との距離の積で定め、各点の速度ベクトルの向きは、固定点と点とを結んだ線の接線方向ベクトルの和が定める向きとする。このときに足し合わせるベクトルの重みは固定点からの距離のk乗(-1乗とか)とする

n<-5


xs<-runif(n)
ys<-runif(n)

X<-seq(from=-1,to=1,by=0.05)
Y<-X
xs<-sample(X,n)
ys<-sample(Y,n)
#xs<-c(-8,0,6)
#ys<-c(-7,6,1)
X<-seq(from=-1,to=1,by=0.05)
Y<-X

XY<-expand.grid(X,Y)
dxdt<-dydt<-rep(0,length(XY[,1]))
# 距離のk乗に比例した速度
k<--2
for(i in 1:length(dydt)){
	tmp<-1
	for(j in 1:n){
		tmpr<-sqrt((XY[i,1]-xs[j])^2+(XY[i,2]-ys[j])^2)
		tmp<-tmp*tmpr
		if(tmpr!=0){
			dxdt[i]<-dxdt[i]-tmpr^k*(XY[i,2]-ys[j])/tmpr
			dydt[i]<-dydt[i]+tmpr^k*(XY[i,1]-xs[j])/tmpr
		}
	}
	tttr<-dxdt[i]^2+dydt[i]^2
	if(tttr!=0){
		dxdt[i]<-dxdt[i]*tmp/sqrt(tttr)
		dydt[i]<-dydt[i]*tmp/sqrt(tttr)
	}
}

par(mfcol=c(2,2))
image(X,Y,matrix(dxdt,ncol=length(X)),col = terrain.colors(50),main="dxdt")

contour(X,Y,matrix(dxdt,ncol=length(X)),add=TRUE,nlevels=100)
for(i in 1:n){
	abline(v=xs[i],col="blue")
}
for(i in 1:n){
	abline(h=ys[i],col="blue")
}

image(X,Y,matrix(dydt,ncol=length(X)),col = terrain.colors(50),main="dydt")

contour(X,Y,matrix(dydt,ncol=length(X)),add=TRUE,nlevels=100)
for(i in 1:n){
	abline(v=xs[i],col="blue")
}
for(i in 1:n){
	abline(h=ys[i],col="blue")
}

image(X,Y,matrix(sqrt(dxdt^2+dydt^2),ncol=length(X)),col = terrain.colors(50),main="velocityNorm")

contour(X,Y,matrix(sqrt(dxdt^2+dydt^2),ncol=length(X)),add=TRUE,nlevels=100)
for(i in 1:n){
	abline(v=xs[i],col="blue")
}
for(i in 1:n){
	abline(h=ys[i],col="blue")
}