保存則と実在則を満足する世界のアトラクタを描く

  • 予定パッケージ名
    • SphereAndSimplex
  • 関数名
    • DrawAttractor
  • タイトル
DrawAttractor
  • 説明
保存則と実在則を満足する世界のアトラクタを描く
n個の複数の変数があり、それが(適当な「分子量」を持ち、その分子量で補正することにより、sum x_i=1を満足している(保存則)
また、x_i > 0(実在則)を満足している
そのようなn個の変数の値セットをを保存則があることからn-1次元空間の点と考える
n-1次元空間には、np個のm<=n-1次元球が存在しており、この球面の集合は、「アトラクティブ」である。n-1次元空間のすべての点は、このnp個の球面に向かい、かつ、np個の球面には、それ固有の回転運動があり、球面に吸い込まれつつ、その球面運動に収束していく
球面の外部の点は外側から収束し、内部の点は内側から収束する。
このようにしてできたn-1次元空間の運動をn-1正単体に縮め、平行移動させることによって、n個の保存則と実在則を満たす変数の動きとして表す
  • 使用例
DrawAttractor(n=4,np=3)
  • ソース
DrawAttractor<-function(n=3,np=2,p=NULL,Nrep=10,initx=NULL,Niter=1000,dt=0.01,Rs=NULL,dimRots=NULL,Rots=NULL,RotsSub=NULL,fracattraction=1,k2=3){
	if(is.null(p)){
		p<-matrix(rnorm(np*n),np,n)
	}

	if(is.null(Rs)){
		# 回転球の半径が、点間距離の1/2よりみじかくなるようにする
		# そうすると、球がつながった形に収束する
		dists<-c()
		for(i in 1:(np-1)){
			for(j in (i+1):np){
				dists<-c(dists,sqrt((p[i,]-p[j,])^2))
			}
		}
		Rs<-runif(np,min=0,max=max(dists)/2)
	}

	if(is.null(dimRots)){
		# 回転は最大でn次元
		dimRots<-sample(2:n,np,replace=TRUE)
	}
	if(is.null(RotsSub)){
		for(i in 1:np){
			RotsSub[[i]]<-NormalBase(n)
		}

	}
	if(is.null(Rots)){
		Rots<-NULL
		for(i in 1:np){
			Rots[[i]]<-diag(rep(1,n))
			Rots[[i]][1:dimRots[i],1:dimRots[i]]<-NormalBase(dimRots[i])
		}
	}
	if(is.null(initx)){
		initx<-matrix(rnorm(Nrep*n),Nrep,n)
	}
	

# 亜次元回転の方向を決める行列
RotsSubInv<-NULL
e.outs<-NULL
Rdts<-NULL

for(i in 1:np){
	RotsSubInv[[i]]<-solve(RotsSub[[i]])	
	e.outs[[i]]<-eigen(Rots[[i]])
	Rdts[[i]]<-(e.outs[[i]][[2]])%*%diag((e.outs[[i]][[1]])^dt)%*%solve(e.outs[[i]][[2]])
}

xssum<-NULL
col<-c()
xssum<-p 
# 集中点も描かせる
# 回転球の中心は「赤(col=2)」
col<-rep(2,np)
for(rep in 1:Nrep){
	xs<-matrix(0,Niter,n)
	xs[1,]<-initx[rep,]
	for(i in 2:Niter){
		v<-rep(0,n)
		vs<-matrix(0,np,n)
		vsRot<-vs
		for(j in 1:np){
			#vs[j,]<--(xs[i-1,]-p[j,])
			# 回転亜球の回転軸に関して回して
			# そののちに回転亜球の回転をし
			# 軸に関して戻し回転をした後
			# オリジナルの点からの移動分を求める
			tmp<-xs[i-1,]-p[j,]
			tmp<-RotsSub[[j]]%*%tmp
			tmpvs<--tmp
			tmpL<-sqrt(sum(tmpvs[1:dimRots[j]]^2))
			tmpvs[1:dimRots[j]]<-(tmpL-Rs[j])/tmpL*tmpvs[1:dimRots[j]]
			vs[j,]<-RotsSubInv[[j]]%*%tmpvs
			vsRot[j,]<-Re((RotsSubInv[[j]])%*%(Rdts[[j]]%*%(tmp)))-(xs[i-1,]-p[j,])			
		}
		# 各中心までの距離
		tmpl<-sqrt(apply(vs^2,1,sum))
		# 各中心までの距離が「与えられた半径」とどれくらい違うか
		#tmpl<-tmpl-Rs
		# 各球の表面からへと近づくが
		# 遠いと速く、近いとゆっくりと近づき、離れない
		stvs<-sign(tmpl)*vs/tmpl^k2
		tmpv<-apply(stvs,2,sum)
		tmpv<-tmpv/sqrt(sum(tmpv^2))

		v<-fracattraction*tmpv*(cumprod(tmpl)[length(tmpl)])*dt
		# 球表面に近いとその球面の回転成分が大きくなる
		vsRot<-vsRot/tmpl
		tmpvRot<-apply(vsRot,2,sum)
		v<-v+tmpvRot
#v<-tmpv*(cumprod(tmpl)[length(tmpl)])^k
		xs[i,]<-xs[i-1,]+v
	}
	
	xssum<-rbind(xssum,xs)
	col<-c(col,rep(rep,Niter))
}
plot3d(xssum[,1],xssum[,2],xssum[,3],col=col)

        
        
        
        
        
        
SimplexXs<-Sphere2Simplex(xssum)
# 正単体は中心が原点なので、全体に(1/(n-1),1/(n-1),...,1/(n-1))だけ平行移動する
SimplexXs<-SimplexXs+1/(n-1)

plot3d(SimplexXs[,1],SimplexXs[,2],SimplexXs[,3],col=col)

list(X=xssum,simplexX=SimplexXs,col=col,n=n,np=np,p=p,Nrep=Nrep,initx=initx,Niter=Niter,dt=dt,Rs=Rs,dimRots=dimRots,Rots=Rots,RotsSub=RotsSub,fracattraction=fracattraction,k2=k2)
}

  • Rdファイル
\name{DrawAttractor}
\alias{DrawAttractor}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
DrawAttractor
}
\description{
保存則と実在則を満足する世界のアトラクタを描く
n個の複数の変数があり、それが(適当な「分子量」を持ち、その分子量で補正することにより、sum x_i=1を満足している(保存則)
また、x_i > 0(実在則)を満足している
そのようなn個の変数の値セットをを保存則があることからn-1次元空間の点と考える
n-1次元空間には、np個のm<=n-1次元球が存在しており、この球面の集合は、「アトラクティブ」である。n-1次元空間のすべての点は、このnp個の球面に向かい、かつ、np個の球面には、それ固有の回転運動があり、球面に吸い込まれつつ、その球面運動に収束していく
球面の外部の点は外側から収束し、内部の点は内側から収束する。
このようにしてできたn-1次元空間の運動をn-1正単体に縮め、平行移動させることによって、n個の保存則と実在則を満たす変数の動きとして表す
}
\usage{
DrawAttractor(n=3,np=2)
}
%- maybe also 'usage' for other objects documented here.
\details{
%%  ~~ If necessary, more details than the description above ~~
}
\value{
\item{r }{Norm of vector}
\item{t }{n-1 angles}
\item{X}{matrix of No.time points x n in n-dimensional space}
\item{simplexX}{matrix of No.time points *(n-1) in (n-1) dimensional space}
\item{col}{vector of numbers indicating color of Nrep curves}
\item{n}{dimension assigned}
\item{np}{No. of attracting spheres}
\item{p}{coordinates of np centers of attracting spheres}
\item{Nrep}{No. curves}
\item{initx}{starting coordinates of Nrep curves}
\item{Niter}{No.of time points of each curve}
\item{dt}{unit time to calculate coordinate of curves}
\item{Rs}{dadia of attracting spheres}
\item{dimRots}{dimention of attracting spheres}
\item{Rots}{rotation matrices of each attracting sphere}
\item{RotsSub}{rotation matrices to determine rotating surface of each attracting spheres}
\item{fracattraction}{ratio of rotation and attraction. when being 0, no attraction but only rotation}
\item{k2}{influence of distance on velocity}
}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%%  ~~who you are~~
}
\note{
%%  ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
}
\examples{
da<-DrawAttractor()
plot3d(da$X[,1],da$X[,2],da$X[,3],col=da$col)
plot3d(da$simplexX[,1],da$simplexX[,2],da$simplexX[,3],col=da$col)
}
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
\keyword{ ~kwd1 }
\keyword{ ~kwd2 }% __ONLY ONE__ keyword per line