2011-02-01

引き寄せて離さない

アトラクタ R

アトラクタのある空間の物体は、空間全体から「アトラクタ」に向かう
アトラクタに引き寄せられた物体はアトラクタ上を移動して「アトラクタの外」へは動かない
今、空間にスカラー場があるとしよう。それは $E=f(\mathbf{x})$ と位置の $\mathbf{x}$ の関数であって、微分可能であるとする
今、ある定数値 $C$ について $E=C$ を満足するような空間の点の集合を $K=\{\mathbf{x}|E(\mathbf{x})=C\}$ とする
この上の点の物体は上を滑るものとする
- 言い換えるとこの $K$ の上の物体が $\frac{d\mathbf{x}}{dt}$ で運動しているとすると $\frac{dE}{dt}=0$ であるから、
- $\frac{dE}{dt}=\sum \frac{dE}{dx_i} \frac{dx_i}{dt}=0$ となる
- これは、ベクトル $(\frac{dE}{dx_i})$ とベクトル $(\frac{dx_i}{dt})$ とが垂直であることも意味する
では上でない点の物体はどこへ向かって動いているだろうか
- $K$ のごく近傍では $K$ 上の点と同様に、 $K$ 上を滑るような動き成分を持つだろう（物体の運動についても空間で微分可能であると仮定している）
- しかしそれだけでなく、その動きと異なる（垂直な）成分の動きも持ち、それは $(\frac{dE}{dx_i})$ というベクトルに平行な成分であって、 $K$ に近づく方向だろう
このようにして考えてきたと外の領域について、上記の条件をまんぞくした上で、酔歩させるのが、次のソース
- $K$ に近ければ近いほど（ $E$ が $C$ に近いほど)、 $\frac{dE}{d\mathbf{x}}$ 成分が小さくなるようにして、 $\frac{dE}{d\mathbf{x}}$ と垂直な方向への動きを酔歩にしたもの
- また、アトラクタは $E=C$ は点であってもよいし、 $n-1$ 次元多様体であってもよく
- アトラクタは、非連続でもよさそうだ
- ローレンツアトラクタは、原点以外の２固定点が「アトラクティブ」な「多様体」であり、「アトラクタ」が点であるとき、そこでは動きはないので、これらは固定点となっている。また、原点は、「アトラクティブ」でない固定点である。…とこうなっているのかと思われる…
ここから「周回軌道」の話に戻さないと寄り道が長くなりすぎる
- 多様体表面の運動は、さて、無限な軌道が引けるのか、それは周回軌道の場合だけなのかそうでないのか、周回しない無限軌道で多様体表面の軌道が引けるとき、多様体には(位相的)特徴としてどうなるのか、さらに、そのような「アトラクティブ」な多様体を決める因子-因子関係を決めるのは何か、と。

d<-3

# E 
potential<-function(x,ks=NULL){
	if(is.vector(x))x<-t(as.matrix(x))
	if(is.null(ks)){
		k<-6
		ks<-matrix(rnorm(length(x[1,])*(k+1)),ncol=k+1)
	}
	
	ret<-rep(0,length(x[,1]))
	for(i in 1:length(ret)){
		for(j in 1:length(x[i,])){
			ret[i]<-ret[i]+sum(x[i,j]^(0:(length(ks[j,])-1))*ks[j,])
		}
	}

	list(ks=ks,v=ret)
}

#x<-seq(from=-10,to=10,length.out=5)
#y<-x
#z<-y
#xy<-expand.grid(x,y,z)
#out<-potential(xy)
#z<-out$v

partialdif<-function(x,ks){
	ret<-sum(ks[2:length(ks)]*(1:(length(ks)-1))*x^(0:(length(ks)-2)))
	
}

#plot3d(xy[,1],xy[,2],z)


Niter<-1000
xs<-matrix(0,Niter,d)
xs[1,]<-runif(d)
Initx<-xs[1,]
out<-potential(as.vector(xs[1,]))
InitP<-out$v
AtractP<-InitP*0.8
dt<-0.01
dp<-0.01

Nrep<-10
xsall<-NULL
col<-c()
for(rep in 1:Nrep){
	xs<-matrix(0,Niter,d)
	xs[1,]<-Initx+rnorm(d)*0.05
col<-c(col,rep(rep,Niter))


for(i in 2:Niter){
	tmpP<-potential(xs[i-1,],out$ks)$v
	print(tmpP)
	difP<-tmpP-AtractP
	tmp<-rep(0,d)
	for(j in 1:d){
		tmp[j]<-partialdif(xs[i-1,j],out$ks[j,])
	}
	v1<-runif(d)*0.01
	select<-which(abs(tmp)==max(abs(tmp)))
	tmps<-sum(tmp[-select]*v1[-select])
	v1[select]<-(-tmps)/tmp[select]
	
	v2<-v1-tmp*difP*dt*dp
		print(sqrt(sum(v1^2)))
	print(sqrt(sum(v2^2)))
	print(difP)

	xs[i,]<-xs[i-1,]+v2
}

xsall<-rbind(xsall,xs)
#plot3d(xs[,1],xs[,2],xs[,3],col=gray((length(xs[,1]):1)/length(xs[,1])))
}
plot3d(xsall[,1],xsall[,2],xsall[,3],col=col)

では空間に何点かあり、そこからの距離に反比例するポテンシャルの和が空間のポテンシャルであるとして、そのポテンシャルの値が定数である面を「アトラクティブ」な多様体とし、その多様体上では多様体表面を滑り、多様体の外では、ポテンシャルを変えない方向の速度とポテンシャルを「アトラクティブ」多様体への垂直成分の速度との和であるようにすると、次のような絵とソースになる

d<-3

Nctr<-4
library(MCMCpack)

Ctrs<-rdirichlet(Nctr,rep(1,d))
#Ctrs[1,]<-rep(0,d)
W<-runif(Nctr)

# E はすべての中心とのW/距離

calcE<-function(x,Ctrs,W){
	ret<-0
	for(i in 1:length(Ctrs[,1])){
		#tmp<-sqrt(sum((x-Ctrs[i,])^2))
		tmp<-sum((x-Ctrs[i,])^2)
		if(tmp==0){
		}else{
			ret<-ret+W[i]/tmp
		}
	}
	ret
}

E0<-1
Emean<-0
numRp<-100
randpts<-rdirichlet(numRp,rep(1,d))
for(i in 1:numRp){
	tmpE<-calcE(randpts[i,],Ctrs,W)
	Emean<-Emean+tmpE
}
E0<-Emean/numRp
#E0<-max(Emean)
Niter<-500
Nrepeats<-100
col<-c()
xsall<-NULL

for(rep in 1:Nrepeats){
	col<-c(col,rep(rep,Niter))
xs<-matrix(0,Niter,d)
xs[1,]<-runif(d)
#xs[1,]<-c(-1,rep(0,d-1))
dt<-0.02
dp<-0.5
Es<-rep(0,Niter)
Es[1]<-calcE(xs[1,],Ctrs,W)
for(i in 2:Niter){
	tmpE<-calcE(xs[i-1,],Ctrs,W)
	deltaE<-rep(0,d)
	for(j in 1:d){
		for(k in 1:length(Ctrs[,1])){
			 deltaE[j]<-deltaE[j]+W[k]*(-1)*(sum((xs[i-1,]-Ctrs[k,])^2))^(-2)*2*(xs[i-1,j]-Ctrs[k,j])
		}
	}
	print(Ctrs)
	print(deltaE)
	print(xs[i-1,])
	print(deltaE/xs[i-1,])
	v<-runif(d)
	select<-which(abs(deltaE)==max(abs(deltaE)))
	v[select]<--sum(v[-select]*deltaE[-select])/deltaE[select]
	#print(deltaE * v)
	
	v<-v-(tmpE-E0)*deltaE*dp
	v<-v/sqrt(sum(v^2))*dt
	xs[i,]<-xs[i-1,]+v
	Es[i]<-calcE(xs[i,],Ctrs,W)
}
xsall<-rbind(xsall,xs)
}


xlim<-ylim<-zlim<-c(min(xsall,Ctrs,E0),max(xsall,Ctrs,E0))
plot3d(c(xsall[,1],Ctrs[,1]),c(xsall[,2],Ctrs[,2]),c(xsall[,3],Ctrs[,3]),col=c(col,rep(2,length(Ctrs[,1]))))

plot(Es,ylim=c(min(Es,E0),max(Es,E0)))
abline(h=E0)