速度ベクトルの動きを回転行列にしてシミュレーション

• 前の記事を用いて、離散的取扱いのために「飛び出さ」ないようにしてシミュレーションしてみた
• 飛び出さないよう、高速になると、１ステップの時間を短くしている(スローモーション化)ので、本当は「さーっ」と通り過ぎるだろうあたりで収束してしまっているようなプロットになっているけれども。。。

n<-3

Niter<-1000
dx<-exp(1)/1000

Ys<-matrix(0,Niter,n+1)

loop<-TRUE
while(loop){
	Ys[1,]<-c(rnorm(n),0)*0.001
	if(sum(Ys[1,1:(n-2)]*Ys[1,2:(n-1)])<0)loop<-FALSE
}
#Ys[1,]<-c(-1,1,0.5,0)
shogen<-sign(Ys[1,])
TmpFx<-function(ys,Equal=FALSE,a=0.1){
	tmp<-0
	ret<-0
	for(i in 1:(length(ys)-1)){
		tmp<-tmp+ys[i]*ys[i+1]
	}
	ret<--tmp/ys[length(ys)]
	if(!Equal){
		ret<-ret+sign(ys[length(ys)])*a
	}
	ret
}


Ys[1,n+1]<-TmpFx(Ys[1,1:n])
alpha<-0.01
for(i in 2:Niter){
	#V<-c(Ys[i-1,2:(n+1)],0)
	V<-Ys[i-1,2:(n+1)]
	#V<-Vs*dx
	P<-Ys[i-1,1:n]
	R<-XYRotation(P,V)
	Rinv<-solve(R)
	Vnorm<-sqrt(sum(V^2))
	Pnorm<-sqrt(sum(P^2))
	theta<-atan(Vnorm/Pnorm)
	# 半径が大きく、移動ベクトルも大きいと
	# １ステップの移動距離が大きくなりすぎるので
	# １ステップあたりの移動距離を等しくする
	theta<-theta/Vnorm*dx
	#print(theta)
	tmpRot<-diag(n)
	tmpRot[1,1]<-cos(theta)
	tmpRot[1,2]<--sin(theta)
	tmpRot[2,1]<-sin(theta)
	tmpRot[2,2]<-cos(theta)
	
	tmpP<-R%*%P
	tmpP<-tmpRot%*%tmpP
	tmpP<-Rinv%*%tmpP
	Ys[i,1:n]<-tmpP
	
	#Vst<-V/Vnorm
	#Ys[i,]<-Ys[i-1,]+dx*Vst
	Ys[i,n+1]<-TmpFx(Ys[i,1:n],Equal=TRUE)
	#zeroq<-sum(Vst[1:n]*Ys[i-1,1:n])
	#print(zeroq)
	#if(abs(zeroq)>0.0001)print(zeroq)
	tmpshogen<-sign(Ys[i,])
	#if(sum((shogen-tmpshogen)^2)!=0){
	#	print("chage shogen")
	#	tmpshogen2<-abs(shogen-tmpshogen)*(-0.5)+1
	#	tmpshogen3<-(tmpshogen2-1)*(-1)
	#	Ys[i,]<-Ys[i,]*tmpshogen2+c(rep(alpha,n)*tmpshogen3[1:n],0)
	#}
	
}
#matplot(Ys[,1:n],type="l")

plot(as.data.frame(Ys))

library(rgl)
plot3d(Ys[,1:3],col=rainbow(1000),xlim=range(Ys[,1:3]),ylim=range(Ys[,1:3]),zlim=range(Ys[,1:3]))

matplot(Ys[,1:n],type="l")