2011-01-20

和が保存、積が保存

変数がNx個
ある基本単位があって、それがいくつか集まって、変数X_iの１個ができるとする。そのユニット数をns
X_iが反応するときには、ある数単位psで反応する。１個でもよいし複数でもよい。この数は、反応の速さを決めるときに[X_i]^数、とべき乗で効いてくるものとする
すべてのX_iはX_{i-1}と反応して増え、X_{i+1}と反応して減るものとする
増える反応、減る反応では、「基本単位」が変わらないものとする（和保存の法則）
- $S=\sum_i^{Nx} n_i X_i$ は一定
そのうえで、反応は次の速度で進むものとする
$\frac{dX_i}{dt}=-a_i X_i^{p_i}X_{i+1}^{p_{i+1}} + b_i X_i^{p_i}X_{i-1}^{p_{i-1}}$
和保存の法則から $a_i ns_i = b_{i+1}ns_{i+1}$ が成り立つ。これから $\prod_i^{Nx} a_i = \prod_i^{Nx} b_i$ が成り立つ
今、この反応系に自明な固定点 $\forall i X_i=0$ 以外の固定点があるとすれば $\forall i \frac{dX_i}{dt}=0$
これを満足するのは、 $Nx$ が奇数なら $a_i ns_i = b_{i+1}ns_{i+1}$
偶数なら、 $\prod_{i mode 2 =0} a_i = \prod_{i mode 2 =0} b_i$ かつ $\prod_{i mode 2 =1} a_i = \prod_{i mode 2 =1} b_i$
このとき次の量が保存される
- $\sum_i^{Nx} u_i(X_i)$
- ただし、
  - $u_i(X_i)=\frac{c_i}{1-p_i}X_i^{1-p_i}; p_i \not =1$
  - $u_i(X_i)=c_i \log(X_i); p_i =1$
  - $c_i b_i = c_{i-2} a_{i-2}$
特に、のとき
- $\sum_i^{Nx} u_i(X_i)=\sum_i^{Nx} c_i \log(X_i)=\log(\prod_i^{Nx} X_i^{p_i})$ が保存される
さらに $\forall i a_i=b_i=1$ のとき $\forall i,j c_i=c_j$ となって $\log(\prod_i^{Nx} X_i)$ が、つまり $\prod_i^{Nx} X_i$ が保存される（積の保存）

# No. variables
Nx<-4
# Units of Xs
ns<-sample(1:3,Nx,replace=TRUE)
#ns<-rep(1,Nx)
# No. involved individuals
ps<-sample(1:5,Nx,replace=TRUE)
#ps<-rep(1,Nx)
# Coefficiency of reactions
#Nx が奇数のときには
#prod(as)=prod(bs)かつ
#ns[i]as[i]=ns[i+1]bs[i+1]を満足する
#Nx が偶数のときには
#prod(as)=prod(bs)かつ
#ns[i]as[i]=ns[i+1]bs[i+1]かつ
#prod(as(偶数))=prod(bs(偶数))かつ
#prod(as(奇数))=prod(bs(奇数))を満足する

as<-runif(Nx)
bs<-as

if(Nx%%2==1){
	for(i in 1:Nx){
		pre<-i-1
		if(pre==0)pre<-Nx
		bs[i]<-ns[pre]/ns[i] * as[pre]
	}

}else{
	for(i in 2:Nx){
		pre<-i-1
		if(pre==0)pre<-Nx
		bs[i]<-ns[pre]/ns[i] * as[pre]
	}
	acumprod<-bcumprod<-1
	
	for(i in 1:((Nx/2)-1)){
		acumprod<-acumprod*as[i*2]
		bcumprod<-bcumprod*bs[i*2]
	}
	as[Nx]<-bcumprod*bs[Nx]/acumprod
	bs[1]<-ns[Nx]/ns[1]*as[Nx]

}
cumprod(as)[Nx]/cumprod(bs)[Nx]
cumprod(as[which((1:Nx)%%2==0)])/cumprod(bs[which((1:Nx)%%2==0)])
cumprod(as[which((1:Nx)%%2==1)])/cumprod(bs[which((1:Nx)%%2==1)])




# 固定点は？

ys<-rep(runif(1)*1,Nx)
if(Nx%%2==1){
	inds<-c(1,2,3)
	while(inds[3]!=1){
		ys[inds[3]]<-bs[inds[2]]/as[inds[2]]*ys[inds[1]]
		inds<-inds+2
		inds[which(inds>Nx)]<-inds[which(inds>Nx)]-Nx
	}
}else{
	inds<-c(1,2,3)
	inds2<-c(2,3,4)
	while(inds[3]!=1){
		ys[inds[3]]<-bs[inds[2]]/as[inds[2]]*ys[inds[1]]
		inds<-inds+2
		inds[which(inds>Nx)]<-inds[which(inds>Nx)]-Nx
		ys[inds2[3]]<-bs[inds2[2]]/as[inds2[2]]*ys[inds2[1]]
		inds2<-inds2+2
		inds2[which(inds2>Nx)]<-inds2[which(inds2>Nx)]-Nx
		
	}

}
fixxs<-ys^1/ps






#as<-rep(1,Nx)
#bs<-rep(1,Nx)
# 保存量の係数
cs<-rep(1,Nx)
cs[1]<-runif(1)*10
cs[2]<-runif(1)*10
cs
#counter<-1
if(Nx%%2==1){
	prei<-1
	i<-3
	while(i!=1){
	#while(counter<100){
		#counter<-counter+1
		cs[i]<-as[prei]/bs[i]*cs[prei]
		prei<-i
		i<-i+2
		if(i>Nx)i<-i-Nx
		if(prei>Nx)prei<-prei-Nx
	}
}else{
	prei<-1
	prei2<-2
	i<-3
	i2<-4
	#cs[2]<-runif(1)
	while(i!=1){
	#while(counter<100){
		#counter<-counter+1
		cs[i]<-as[prei]/bs[i]*cs[prei]
		cs[i2]<-as[prei2]/bs[i2]*cs[prei2]
		prei<-i
		i<-i+2
		prei2<-i2
		i2<-i2+2
		if(i>Nx)i<-i-Nx
		if(i2>Nx)i2<-i2-Nx
		if(prei>Nx)prei<-prei-Nx
		if(prei2>Nx)prei2<-prei2-Nx
	}

}
cs

Niter<-10000
dt<-1
xs<-matrix(0,Niter,Nx)
xs[1,]<-runif(Nx)
xs[1,]<-fixxs
dev<-rnorm(Nx)/10
dev<-dev-mean(dev)
xs[1,]<-xs[1,]+runif(1)/10
library(MCMCpack)
#xs[1,]<-rdirichlet(1,rep(1,Nx))

E1<-E2<-rep(0,Niter)
E1[1]<-sum(ns*xs[1,])
for(i in 1:Nx){
	tmp<-0
	if(ps[i]!=1){
		tmp<-cs[i]/(1-ps[i])*xs[1,i]^(1-ps[i])
	}else{
		tmp<-cs[i]*log(xs[1,i])
	}
	E2[1]<-E2[1]+tmp
}
for(i in 2:Niter){
	for(j in 1:Nx){
		pre<-j-1
		post<-j+1
		if(pre<1)pre<-Nx
		if(post>Nx)post<-1
		xs[i,j]<-xs[i-1,j]-(as[j]*xs[i-1,j]^ps[j]*xs[i-1,post]^ps[post]-bs[j]*xs[i-1,j]^ps[j]*xs[i-1,pre]^ps[pre])*dt
	}
	E1[i]<-sum(ns*xs[i,])
	for(j in 1:Nx){
		tmp<-0
		if(ps[j]!=1){
			tmp<-cs[j]/(1-ps[j])*xs[i,j]^(1-ps[j])
		}else{
			tmp<-cs[j]*log(xs[i,j])
		}
		E2[i]<-E2[i]+tmp
	}

}

plot(as.data.frame(xs))
library(rgl)
plot3d(xs[,1],xs[,2],xs[,3])
matplot(xs,type="l")