曲率・捩率…

• ２年ぶりに再読する『じっくりと学ぶ曲線と曲面』

じっくり学ぶ曲線と曲面―微分幾何学初歩

• ２次元曲線は曲率で定義
• ３次元曲線は曲率と捩率(れいりつ)で定義
• n次元曲線は、曲がりに関するn-1個の変数(第一が曲率、第二が捩率、第三以降は、一般化されたもの)で定義
• n次元の曲がりを表わすn-1個の変数はnxn行列を以下のように作り、それはフルネ=セレの一般次元の公式
• 
• こちら(Wiki)
• フルネ=セレの行列成分が時間に関して(曲線上の位置に関して)一定であれば、それは、らせんになる

# フルネ=セレ

n<-5

M<-matrix(0,n,n)

c<-runif(n-1)
#c<-rep(1,n-1)
#c<-1^((n-1):1)
#c<-rnorm(n-1)
c<-1.1^((n-1):1)
for(i in 1:length(c)){
	M[i,i+1]<-c[i]
	M[i+1,i]<--c[i]
}
#M<-M+diag(rep(1,n))
Niter<-10000
dt<-0.001
xs<-matrix(0,Niter,n)
xs[1,]<-runif(n)*10
e.out<-eigen(M)
u<-e.out[[2]]
v<-solve(u)
s<-diag(e.out[[1]])

# どうひょうこうを作る
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
}

dhk<-NormalBase(n)
v<-dhk[1,]

for(i in 2:Niter){
	#xs[i,]<-M%*%xs[1,]
	xs[i,]<-xs[i-1,]+v*dt
	#xs[i,]<-(u%*%s^(i*dt)%*%v)%*%xs[1,]
	#xs[i,]<-(u%*%v^(i*dt)%*%v)%*%xs[1,]
	#dhk<-(u%*%diag(s^(i*dt))%*%v)%*% dhkInit
	dhk<-(M %*% dhk ) *dt + dhk
	# 単位ベクトル補正
	ll<-apply(dhk^2,1,sum)
	dhk<-((dhk)/sqrt(ll))
	v<-dhk[1,]
	#print(sum(v^2))
}

library(rgl)
plot3d(xs[,1],xs[,2],xs[,3],col=rainbow(Niter))

matplot(xs,type="l")