2010-12-28

位置変化(物理)と構成変化(化学)

R 拡散反応拡散系偏微分

こちらで拡散を扱っている
- 拡散とは、ある地点にある物質が、そこでの「濃度」に比例して「移動可能な方向」へ位置変化すること
- なので、状態推移行列を使って表すと、行列には拡散のしやすさを表す係数が入って、それにある時点の濃度をかけるような時系列変化式になる
- ある地点から「移動可能な方向」への位置変化のしやすさがどの地点でも同じとき、「ある地点からの出」と「ある地点への入り」とは、濃度の勾配として差引することになる。これが、濃度分布の微分(濃度勾配)を使って拡散を記述する理由である
- 濃度勾配を格子モデルで表したのがこちらの記事の８周囲方向・４周囲方向の話
- 全方向についての出入りをまとめたものがラプラシアン(こちら)を作用させた結果
さて、拡散は、ある物質の位置変化が、その物質の濃度(の１乗)に比例した現象だった
化学反応に代表される、物質の「構成変化」について考える
- $A+B \to C+D$ は、Aの濃度とBの濃度の積に比例して進行する。拡散が着目物質の濃度にのみ比例していたのと様子が異なっている
- $A+B^2 \to C+D$ という構成変化の場合は、Aの濃度と「Bの濃度の２乗」の積に比例して進行する。
- 今、一般に、 $N_S$ 種類の物質があって、それ以外の物質はないとする。
- このような世界の「構成変化」は $\sum_{i=1}^{N_S} n_i S_i \to \sum_{i=1}^{N_S} m_i S_i$ という式で表される。
- 拡散のときには、拡散しやすさを表す拡散係数があったが、構成変化の場合のそれを、「反応係数」と呼ぶことにすれば、話がしやすくなる
さらに一般的に考えるなら、ある物質の移動のしやすさが、その物質の濃度ではなく、他の物質の濃度に比例するような現象(誘引・排斥とか？）も考えられるし、「構成変化」の進行が、構成変化の左辺の濃度の積に比例するとは限らず、構成しない要素の濃度に比例させたりすることもできる。また、濃度によらず、消滅する(胎生致死)とか、定期的な供給であるとか、を入れることもできる。
どこまでを、モデルとして入れ込むかが考えどころだが、一般的に記述したうえで、そこに制約を入れると、モデルがなんなのかの見通しが良くなることは確か。
反応をこの方式で進めるソースが以下。ここでは反応の前後で分子数が変化しないような反応を複数発生させている

# X,Y,Z
library(MCMCpack)
k<-4
xyz<-rdirichlet(1,rep(1,k))

Nr<-6
Nmol<-rpois(Nr,6)+1
pre<-matrix(0,Nr,k)
post<-matrix(0,Nr,k)
for(i in 1:Nr){
	prepre<-sample(1:k,Nmol,replace=TRUE)
	postpost<-sample(1:k,Nmol,replace=TRUE)
	for(j in 1:k){
		pre[i,j]<-length(which(prepre==j))
		post[i,j]<-length(which(postpost==j))
	}
}

#pre<-matrix(sample(0:3,Nr*k,replace=TRUE),Nr,k)
#post<-matrix(sample(0:3,Nr*k,replace=TRUE),Nr,k)
Rfw<-runif(Nr)
Rrev<-runif(Nr)


Niter<-1000
d<-matrix(0,Niter,k)

d[1,]<-xyz
for(i in 2:Niter){
	tmppred<-d[i-1,]
	for(j in 1:Nr){
		lnfreq<-log(tmppred)
		tmpd<-rep(0,k)

		kfw<-exp(sum(lnfreq*pre[j,]))*Rfw[j]
		krev<-exp(sum(lnfreq*post[j,]))*Rrev[j]
		tmpd<-tmpd-kfw*pre[j,]+kfw*post[j,]-krev*post[j,]+krev*pre[j,]
		tmppred<-tmppred+tmpd
	}
	d[i,]<-tmppred
}

matplot(d,type="l")