2011-02-10

曲率が曲線を決める

こちらでローレンツアトラクタをいじった
こちらで曲線と曲線におけるMoving frame、さらには、観測曲線座標データから、Moving frameを算出することを書いた。曲率も計算で来た
両者を合わせる→ローレンツアトラクタの曲率を計算し、第一曲率の正負と第二曲率の正負でローレンツアトラクタを塗り分けてみた
第一曲率がゼロの点は、グラフの変曲点のようなもので、局所的に直線になるところ
第一曲率がゼロで、その前後の曲率の正負が入れ替われば、曲線は、曲がり向きを変える
第一曲率がゼロでその前後で符号を変えないならば、曲線は一度直線状になって、もとの曲がり向きに戻る
第二曲率は捩れのパラメタなので、それがゼロの点では、よじれがなくなり、前後で符号が変われば、捩れの向きが変わることを示し、前後で冨符号が変わらなければ、元の捩れ方に戻ることを表す。
掲載図のローレンツアトラクタでは、第一、第二曲率が（正、正）を黒、（負、負）を青で描いてある
ほぼ同じ場所で、曲がり向き（第一曲率）と捩れ向き（第二曲率）が交代していることがわかる。（正、正）の黒から、（負、正）の緑を経て、（負、負）の黒に移行している。
この緑の領域で、「振り分け」が起きている。これが「ローレンツアトラクタの八の字ねじれ」
こちらでは、速度ベクトルの成分をx,y,zの３軸方向の正負で８区分に分けて色塗りをしたが、「本質」は軸のとり方によるのではなくて、「曲線の特性」である「曲率」によってわかりやすくなる様子が、見て取れているように思える…

# 回転行列を作る
NormalBase<-function(n){ # 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
}

# Moving frameを算出するための関数


# 直交単位ベクトルを取り出す
PerpendicVector<-function(v1,v2){
	r1<-sqrt(sum(v1^2))
	r2<-sqrt(sum(v2^2))
	ret<-rep(0,length(v2))
	ct<-0
	if(r1*r2!=0){
		ip<-sum(v1*v2)
		ct<-ip/(r1*r2)
		ret<-v2/r2-v1/r1*ct
		ret<-ret/sqrt(sum(ret^2))
	}
	list(vector=ret,cos=ct)
}
# Moving frameを作る、曲率も出す
MovingFrame<-function(xs,n=NULL){
	if(is.null(n)){
		n<-length(xs[1,])
	}
	Es<-array(0,c(n,length(xs[,1]),length(xs[1,])))
	Cs<-matrix(0,n,length(xs[,1]))
	Es[1,,]<-(xs)/sqrt(apply(xs^2,1,sum))
	current<-xs
	for(i in 2:n){
		tmp<-apply(current,2,diff)
		current<-tmp
		for(j in 1:length(tmp[,1])){
			tmpout<-PerpendicVector(Es[i-1,j,],tmp[j,])
			Es[i,j,]<-tmpout$vector
			Cs[i-1,j]<-tmpout$cos
		}
	}
	list(Es=Es,Cs=Cs)
}


ct<-c(6*sqrt(2),6*sqrt(2),27)
col<-c()
	Niter<-10000
	#col<-c(col,rep(rep,Niter))
	#col<-c(col,c(1,rep(rep,Niter-1)))
dt<-0.004
xs<-matrix(0,Niter,3)

xs[1,]<-runif(3,min=-10,max=10)+ct
#xs[1,3]<-abs(xs[1,3])*2

#xs[1,]<-c(2*sqrt(2)*3,2*sqrt(2)*3,27)
vs<-runif(3)*50
vs<-rep(0,3)
for(i in 2:Niter){
	xs[i,1]<-xs[i-1,1]+(-10*xs[i-1,1]+10*xs[i-1,2]+vs[1])*dt
	xs[i,2]<-xs[i-1,2]+(28*xs[i-1,1]-xs[i-1,2]-xs[i-1,1]*xs[i-1,3]+vs[2])*dt
	xs[i,3]<-xs[i-1,3]+(-8/3*xs[i-1,3]+xs[i-1,1]*xs[i-1,2]+vs[3])*dt
}
#tmpx1<-xs[,1]/sqrt(2)-xs[,2]/sqrt(2)
#tmpy1<-xs[,1]/sqrt(2)+xs[,2]/sqrt(2)
#xs[,1]<-tmpx1
#xs[,2]<-tmpy1





xs<-Re(xs)


vs<-apply(xs,2,diff)

Es<-MovingFrame(vs)
matplot(t(Es[[2]]),type="l")

plot(as.data.frame(t(Es[[2]])))

library(rgl)
open3d()
xlim<-ylim<-zlim<-c(min(xs),max(xs))
#plot3d(xssum[,1],xssum[,2],xssum[,3],col=rainbow(1000),cex=0.1)
#plot3d(xssum[,1],xssum[,2],xssum[,3],col=gray(col/Nrep),cex=0.1,xlim=xlim,ylim=ylim,zlim=zlim)
# 第一曲率が正なら黒、負なら赤
col<-rep(0,Niter)
col[which(Es[[2]][1,]>0)]<-1
col[which(Es[[2]][1,]<0)]<-2

plot3d(xs[,1],xs[,2],xs[,3],col=col,xlim=xlim,ylim=ylim,zlim=zlim,cex=0.1)

# 第二曲率が正なら黒、負なら赤
col<-rep(0,Niter)
col[which(Es[[2]][2,]>0)]<-1
col[which(Es[[2]][2,]<0)]<-2

plot3d(xs[,1],xs[,2],xs[,3],col=col,xlim=xlim,ylim=ylim,zlim=zlim,cex=0.1)

# 第一曲率、第二曲率が、(正、正）ならcol=1（黒）、（正、負）ならcol=2（赤）、（負、正）ならcol=3（緑）、（負、負）ならcol=4（青）

# 第二曲率が正なら黒、負なら赤
col<-rep(0,Niter)
col[which(Es[[2]][1,]>0 & Es[[2]][2,]>0)]<-1
col[which(Es[[2]][1,]>0 & Es[[2]][2,]<0)]<-2
col[which(Es[[2]][1,]<0 & Es[[2]][2,]>0)]<-3
col[which(Es[[2]][1,]<0 & Es[[2]][2,]<0)]<-4

plot3d(xs[,1],xs[,2],xs[,3],col=col,xlim=xlim,ylim=ylim,zlim=zlim,cex=0.1)