この世の偏微分方程式を作る

• この日の続き
• n次元空間にnp個の半径Rsの球が相互に重なりあいを持ちながら存在している
• 空間のすべての点は、この球の重なりあいが作る球面に向かう
• この球は、回転運動を持ち、その表面と周囲にその回転運動を起こさせる
• したがって、球の表面でも周囲でも、球が定める空間回転運動をしつつ、「複合球面」に向かう
  • ただし、nが偶数のときの回転は、軸がなく、nが奇数のときの回転は、軸がある
  • それは、こちらで述べた通りで、回転を表す行列がが、を満足するようななるを持つか持たないか(このような点は、回転による不動点)で決まるが、nが偶数のときはこのような点はなく、nが奇数のときにはあることから言える

• ソースはこちら

library(rgl)
 NormalBase<-function(n){
I<-X<-diag(rep(1,n))

thetas<-runif(n*(n-1)/2)*2*pi
T<-matrix(0,n,n)
T[lower.tri(T)]<-thetas

for(i in 1:(n-1)){
for(j in (i+1):n){
R<-I
R[i,i]<-R[j,j]<-cos(T[j,i])
R[i,j]<-sin(T[j,i])
R[j,i]<--R[i,j]
X<-R%*%X
}
}
X
}

n<-4
np<-3
Niter<-10000
dt<-0.01

p<-matrix(rnorm(np*n),np,n)


# 回転球の半径が、点間距離の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)

Rots<-NULL
e.outs<-NULL
Rdts<-NULL
for(i in 1:np){
	Rots[[i]]<-NormalBase(n)
	e.outs[[i]]<-eigen(Rots[[i]])
	Rdts[[i]]<-(e.outs[[i]][[2]])%*%diag((e.outs[[i]][[1]])^dt)%*%solve(e.outs[[i]][[2]])
}
Nrep<-1
xssum<-NULL
col<-c()
xssum<-p 
# 集中点も描かせる
# 回転球の中心は「赤(col=2)」
col<-rep(2,np)
#k<--runif(1)*10
k2<-3

for(rep in 1:Nrep){
	xs<-matrix(0,Niter,n)
	xs[1,]<-rnorm(n)
	xs[1,]<-xs[1,]/sqrt(sum(xs[1,]^2))*5
	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,])
			vsRot[j,]<-Re(Rdts[[j]]%*%(xs[i-1,]-p[j,]))-(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<-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))
}
cex<-rep(0.1,length(col))
cex[1:Nrep]<-3
plot3d(xssum[,1],xssum[,2],xssum[,3],col=col)