- 昨日までの記事で、常微分方程式が発散しないような条件を常微分方程式の階数次元空間の軌道で考えてきている
- 微分係数を指定して離散的にトラッキングすると、「理論は正しい」はずだけれども、誤差が過大になって、描図が芳しくない
- なので、軌道を決めるベクトルによる微小時間の動きを短直線ではなく、小角回転にしてみることとする
- そのためには、空間で任意になっている回転の半径ベクトルと、進行方向ベクトルとから、それに相当する回転の行列が作りたい
- ひとまずその行列を作るために次のようにする
- 半径ベクトルとそれに直交する進行方向ベクトルを(1,0,...,0),(0,1,0,...,0)ベクトル方向に回転する(そのような行列を作る)
- 回転させた後に、第1・第2軸が作る平面における回転(これは行列のうち、第1・2行、第1・2列以外は、単位行列の構成で、第1・2の行と列とはとなるような行列である)をさせ、もとの座標系に逆回転させることとする
Cartesian2AngularX<-function(x){
n<-length(x)
if(n==1){
return(list(r=abs(x),t=acos(sign(x))))
}else if(n==2){
tmpt<-0
if(x[1]!=0){
tmpt<-acos(x[1]/sqrt(sum(x^2)))*sign(x[2])
if(sign(x[2])==0){
if(x[1]>0){
tmpt<-0
}else if(x[1]<0){
tmpt<-pi
}
}
}else{
if(x[2]>0){
tmpt<-pi/2
}else if(x[2]<0){
tmpt<--pi/2
}
}
return(list(r=sqrt(sum(x^2)),t=tmpt))
}
r<-sqrt(sum(x^2))
xst<-x/r
t<-rep(0,n-1)
S<-C<-t
S[n-1]<-xst[n]
if(S[n-1]>1)S[n-1]<-1
if(S[n-1]<(-1))S[n-1]<-(-1)
t[n-1]<-asin(S[n-1])
C[n-1]<-cos(t[n-1])
cumC<-C[n-1]
for(i in (n-2):1){
if(cumC!=0){
S[i]<-xst[i+1]/cumC
if(S[i]>1)S[i]<-1
if(S[i]<(-1))S[i]<-(-1)
t[i]<-asin(S[i])
C[i]<-cos(t[i])
cumC<-cumC*C[i]
}else{
S[i]<-0
t[i]<-asin(S[i])
C[i]<-cos(t[i])
}
}
list(r=r,t=t)
}
VectorRotation<-function(P,inv=FALSE){
n<-length(P)
if(n==1){
return(matrix(c(1),1,1))
}
Q<-Cartesian2AngularX(P)
M<-diag(n)
M2<-M
for(i in 1:length(Q[[2]])){
tmp<-M
tmp[i,i]<--sin(Q[[2]][i])
tmp[i,i+1]<-cos(Q[[2]][i])
tmp[i+1,i]<-cos(Q[[2]][i])
tmp[i+1,i+1]<-sin(Q[[2]][i])
M2<-tmp%*%M2
}
M2<-t(M2)[,n:1]
S <- rep(1,length(P))
S2 <- sign(M2[,1]*P)
S[which(S2<0)] <- -1
M2<-M2*S
if(inv){
M2<-t(M2)
}
M2
}
XYRotation<-function(P,V,inv=TRUE){
R1<-VectorRotation(P,inv=TRUE)
V2<-R1%*%V
tmpV<-V2[2:length(V)]
tmpR2<-VectorRotation(tmpV,inv=TRUE)
R2<-diag(length(V))
R2[2:n,2:n]<-tmpR2
ret<-R2%*%R1
if(!inv)ret<-t(ret)
ret
}
n<-5
P<-rnorm(n)
V<-rnorm(n)
V[n]<--sum(P[1:(n-1)]*V[1:(n-1)])/P[n]
P
V
sum(P*V)
R<-XYRotation(P,V)
P2<-R%*%P
V2<-R%*%V
P2
V2