行列の指数関数

• 指数関数は、微分に関して、特徴的
• その指数関数を行列に拡張する
  • 参考こちら
• そうすると、連立常微分方程式を解くのに使える
• 固有値が(も)ほしい
• ケーリー・ハミルトンの定理(こちら)
• このＰＤＦもよい
• 固有値の実数解・重解虚数解で、回転が入るかどうかが決まる。虚数解はペアで回転
• 次元が奇数だと、回転成分だけにできない
• 回転のへその有無と偶奇とが関連？？？→こちら
• Rでやってみよう
  • 例題から
    • なる行列が作るなる連立微分方程式があって、その初期値がであるとき
    • 本では、と言う
    • これを固有値分解して計算することで、一致することを示そう

n<-2
A<-matrix(c(-1,-4,2,5),2,2)
# 固有値分解
eigen.out<-eigen(A)
# P=V,P^{-1}=U
V<-eigen.out[[2]]
U<-solve(V)
B<-diag(eigen.out[[1]])
# xを適当にとって
x<-seq(from=0,to=1,length=101)
# Yの値の格納庫
Ys<-matrix(0,length(x),n)
# 初期値を与え
Ys0<-c(1,3)
for(i in 1:(length(x))){
# e^{xB},e^{xA}を計算する
	Bex<-diag(exp(eigen.out[[1]]*x[i]))
	Aex<-V%*%Bex%*%U
# Yの値を計算する
	Ys[i,]<-Aex%*%Ys0
}
# 変化の様子をプロットして
matplot(Ys,type="l")
# 本の計算式でも計算してやる
Ys2<-matrix(0,length(x),n)
Ys2[,1]<--exp(x)+2*exp(3*x)
Ys2[,2]<--exp(x)+4*exp(3*x)
# ２つの方法が同じ結果であることを確認する
plot(Ys[,1],Ys2[,1])
plot(Ys[,2],Ys2[,2])

• うまく回っているので、任意の次数、固有値が複素数になったりしても、なんでもありにしてやってみる
  • 複素数が出ても、Ysの値は、実数になることは、複素成分が無視できるくらい小さいことを確認することで裏をとる

# 次数
n<-8
# 適当に正方行列を作って
A<-matrix(rnorm(n^2),n,n)
# 固有値分解
eigen.out<-eigen(A)
V<-eigen.out[[2]]
U<-solve(V)
B<-diag(eigen.out[[1]])
# 計算する
x<-seq(from=0,to=1,length=101)
Ys<-matrix(0,length(x),n)
# 初期値も適当に
Ys0<-runif(n)
# 複素成分の自乗の総和がほぼゼロであることを確認
IMS<-rep(0,length(x))
for(i in 1:(length(x))){
	Bex<-diag(exp(eigen.out[[1]]*x[i]))
	Aex<-V%*%Bex%*%U
	Ys[i,]<-Aex%*%Ys0
	ims<-sum(Im(Ys[i,])^2)
	IMS[i]<-ims
# 複素成分が大きかったら表示させるようにして確認
	if(ims>0.000001)print(ims)
}
plot(IMS)
matplot(Ys,type="l")

• 固有値を指定して行列を作って、その挙動を見てみよう(こちら)
• 行列のべき乗と指数行列

# べき乗
pow.m <- function(A,n){
	# 固有値分解
	eigen.out<-eigen(A)
	# P=V,P^{-1}=U
	V<-eigen.out[[2]]
	U<-solve(V)
	B<-diag(eigen.out[[1]]^n)
	X <- V%*%B%*%U
	return(list(matrix = X, eigen.vs <- eigen.out[[1]]))
}

# 指数行列
exp.m <- function(A,n){
	# 固有値分解
	eigen.out<-eigen(A)
	# P=V,P^{-1}=U
	V<-eigen.out[[2]]
	U<-solve(V)
	B<-diag(exp(eigen.out[[1]]*n))
	X <- V%*%B%*%U
	return(list(matrix = X, eigen.vs <- eigen.out[[1]]))
}