この世の偏微分方程式を作る - ryamadaの遺伝学・遺伝統計学メモ

昨日の記事で正の変量に「質量保存の法則」が成り立つときの偏微分方程式が定める曲線とその上の運動について書いた(こちら)
今日は、場合分けに対応する偏微分方程式の例を作成してみる
$n$ ：変量の数、次元
$\sum_{i=1}^n m_i x_i=C$ を、単純にするために、 $\sum_{i=1}^n x_i=1$ と標準化する
$D=\{\mathbf{x}|\sum_{i=1}^n x_i =1, x_i>0\}$ を領域 $D$ と呼ぶことにする
領域内の点に向かって移動し、そこで静止する運動の例(以下、断らない限りは、例であるものとする)
- $\frac{d\mathbf{x}}{dt}=-(\mathbf{x}-\mathbf{p})$

library(rgl)
library(MCMCpack)
# 領域[tex:D]内の点[tex:\mathbf{p}=(p_i)]に向かって移動し、そこで静止する運動の例
n<-5
p<-rdirichlet(1,rep(1,n))
Niter<-100
dt<-0.01
Nrep<-100
xssum<-NULL
col<-c()
for(rep in 1:Nrep){
	xs<-matrix(0,Niter,n)
	xs[1,]<-rdirichlet(1,rep(1,n))
	for(i in 2:Niter){
		v<--(xs[i-1,]-p)
		xs[i,]<-xs[i-1,]+v*dt
	}
	xssum<-rbind(xssum,xs)
	col<-c(col,rep(rep,Niter))
}
xlim=ylim=zlim=c(0,1)
plot3d(xssum[,1],xssum[,2],xssum[,3],col=col,xlim=xlim,ylim=ylim,zlim=zlim)
matplot(xssum,type="l",ylim=ylim)

領域 $D$ 内の点 $\mathbf{p}=(p_i)$ に向かって移動し、 $\mathbf{p}$ に収束するが、到達しない運動の例。速さが収束点からの距離のk(正の数)乗であるとしてある。

# 領域[tex:D]内の点[tex:\mathbf{p}=(p_i)]に向かって移動し、Pに到達せずに収束する運動の例

Niter<-1000
dt<-0.1

Nrep<-10
xssum<-NULL
col<-c()
k<-runif(1)*10
for(rep in 1:Nrep){
	xs<-matrix(0,Niter,n)
	xs[1,]<-rdirichlet(1,rep(1,n))
	for(i in 2:Niter){
		v<--(xs[i-1,]-p)
		
		v<-v*(sqrt(sum(v^2)))^k
		xs[i,]<-xs[i-1,]+v*dt
	}
	xssum<-rbind(xssum,xs)
	col<-c(col,rep(rep,Niter))
}
xlim=ylim=zlim=c(0,1)
plot3d(xssum[,1],xssum[,2],xssum[,3],col=col,xlim=xlim,ylim=ylim,zlim=zlim)
matplot(xssum,type="l",ylim=ylim)

複数の点からの距離に-k2乗比例して進行方向が決まり、各点への距離の積に比例して速度が決まると、複数の点のどれかに無限に近づく

# 領域[tex:D]内の複数の点[tex:\mathbf{p}=(p_i)]に向かって移動し、速さは距離に比例させる運動の例
n<-5
np<-6
p<-rdirichlet(np,rep(1,n))

Niter<-1000
dt<-0.1

Nrep<-8
xssum<-NULL
xssum<-p # 集中点も描かせる
col<-rep(1,np)
k<--runif(1)*10
k2<-runif(1)*10
for(rep in 1:Nrep){
	xs<-matrix(0,Niter,n)
	xs[1,]<-rdirichlet(1,rep(1,n))

	for(i in 2:Niter){
		v<-rep(0,n)
		vs<-matrix(0,np,n)
		for(j in 1:np){
			vs[j,]<--(xs[i-1,]-p[j,])
		}
		tmpl<-sqrt(apply(vs^2,1,sum))

		stvs<-vs/tmpl^k2
		tmpv<-apply(stvs,2,sum)
		tmpv<-tmpv/sqrt(sum(tmpv^2))

		v<-tmpv*(cumprod(tmpl)[length(tmpl)])
#v<-tmpv*(cumprod(tmpl)[length(tmpl)])^k
		
		xs[i,]<-xs[i-1,]+v*dt
	}
	xssum<-rbind(xssum,xs)
	col<-c(col,rep(rep+1,Niter))
}
xlim=ylim=zlim=c(0,1)
plot3d(xssum[,1],xssum[,2],xssum[,3],col=col,xlim=xlim,ylim=ylim,zlim=zlim)
matplot(xssum,type="l",ylim=ylim)