指数行列2

  • この絵は、固有値がすべて実部0の虚数である行列を連立微分方程式の係数行列として持つ初期値問題の3因子の値変化の3次元プロット
  • 描き方は以下の通り
  • 昨日の記事の続き
  • 行列の固有値を実・虚まぜこぜで指定して、そのような固有値を持つ行列を作ってやろう
    • そうすることで、固有値のパターンごとに、常微分方程式の初期値問題の結果をプロットしてやることにする
  • 周期性が出るのは、虚数固有値がある場合
  • 指数関数性質がなくなって周期性のみになるのは、すべての固有値が実部が0の虚数の場合
  • 行列の固有値を指定して、行列を作る
    • 実数固有値はそのまま1x1行列とし、虚数固有値は、共役複素数による2x2行列を作り、このようにして作った、1x1行列と2x2行列とを対角に並べた行列を作る
    • それに、任意の実数行列とその逆行列とをサンドイッチさせてやると、その行列の固有値は、与えた固有値になっている
####
# 固有値は実数のものと、複素数のものがあり
# 複素数のものは共役なペアを持つものとする
# 実数の固有値は、k-重根であることもあり
# 複素数のペアも、重複があることがあるとする

# 今、異なる値の実数値固有値の数をnR
# その重複度をdupRという長さnRの自然数値ベクトルとし
# 異なる値の虚数値固有値ペアの数をnI
# その重複度をdupIという長さnIの自然数値ベクトルとする
# このとき、固有値の総数nは、n=nR+nI*2

nR<-0
Rs2<-c()
if(nR>0){
	Rs<-c(1,runif(nR-1)*10)
	dupR<-rep(2,nR)

	Rs2<-c()
	for(i in 1:nR){
		Rs2<-c(Rs2,rep(Rs[i],dupR[i]))
	}

}
nI<-3
Is3<-c()
if(nI>0){
	dupI<-rep(1,nI)
# そのノルム
norI<-rep(1,nI)
norI<-runif(nI)
#Is<-complex(real=rnorm(nI),imaginary=rnorm(nI))
Is<-complex(real=rep(0,nI),imaginary=rnorm(nI))

Is<-Is/Mod(Is)*norI
Mod(Is)
Is2<-c()
for(i in 1:nI){
	Is2<-c(Is2,rep(Is[i],dupI[i]))
}
Is3<-c(Is2,Conj(Is2))

}
Is3

Ss<-c(Rs2,Is3)

n<-nR+nI*2

# 
S<-matrix(0,n,n)
for(i in 1:nR){
	S[i,i]<-Rs[i]
}
for(i in 1:nI){
	st<-(i-1)*2+1+nR
	S[st,st]<-Re(Is[i])
	S[st,st+1]<--Im(Is[i])
	S[st+1,st]<-Im(Is[i])
	S[st+1,st+1]<-Re(Is[i])
}
S
library(GPArotation)
V2<-Random.Start(n)
V2<-matrix(rnorm(n^2),n,n)
U2<-solve(V2)

A<-V2%*%S%*%U2

eigen.out<-eigen(A)
eigen.out[[1]]
Ss
> eigen.out[[1]]
[1] 0+0.9917628i 0-0.9917628i 0+0.8862974i 0-0.8862974i 0+0.1676730i 0-0.1676730i
> Ss
[1] 0-0.8862974i 0+0.1676730i 0+0.9917628i 0+0.8862974i 0-0.1676730i 0-0.9917628i
  • こんな風にして作った行列による変化の様子を見る
    • すべて実数解、すべて異なる値

nR<-4
Rs2<-c()
if(nR>0){
	Rs<-c(1,runif(nR-1)*10)
	dupR<-rep(1,nR)

	Rs2<-c()
	for(i in 1:nR){
		Rs2<-c(Rs2,rep(Rs[i],dupR[i]))
	}

}
nI<-0
Is3<-c()
if(nI>0){
	dupI<-rep(1,nI)
# そのノルム
norI<-rep(1,nI)
norI<-runif(nI)
#Is<-complex(real=rnorm(nI),imaginary=rnorm(nI))
Is<-complex(real=rep(0,nI),imaginary=rnorm(nI))

Is<-Is/Mod(Is)*norI
Mod(Is)
Is2<-c()
for(i in 1:nI){
	Is2<-c(Is2,rep(Is[i],dupI[i]))
}
Is3<-c(Is2,Conj(Is2))

}
Is3

Ss<-c(Rs2,Is3)
n<-length(Ss)

# 
S<-matrix(0,n,n)
for(i in 1:nR){
	S[i,i]<-Rs[i]
}
for(i in 1:nI){
	st<-(i-1)*2+1+length(Rs2)
	S[st,st]<-Re(Is[i])
	S[st,st+1]<--Im(Is[i])
	S[st+1,st]<-Im(Is[i])
	S[st+1,st+1]<-Re(Is[i])
}

S
library(GPArotation)
V2<-Random.Start(n)
V2<-matrix(rnorm(n^2),n,n)
U2<-solve(V2)

A<-V2%*%S%*%U2

eigen.out<-eigen(A)
eigen.out

V<-eigen.out[[2]]
U<-solve(V)
B<-diag(eigen.out[[1]])

x<-seq(from=0,to=10,length=101)

Ys<-matrix(0,length(x),n)
Ys0<-runif(n)
Ys0<-rep(0,n)
Ys0[1]<-1
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")
  • すべて同じ値で実数解

nR<-1
Rs2<-c()
if(nR>0){
	Rs<-c(runif(nR)*2)
	dupR<-rep(4,nR)

	Rs2<-c()
	for(i in 1:nR){
		Rs2<-c(Rs2,rep(Rs[i],dupR[i]))
	}

}
nI<-0
Is3<-c()
if(nI>0){
	dupI<-rep(1,nI)
# そのノルム
norI<-rep(1,nI)
norI<-runif(nI)
#Is<-complex(real=rnorm(nI),imaginary=rnorm(nI))
Is<-complex(real=rep(0,nI),imaginary=rnorm(nI))

Is<-Is/Mod(Is)*norI
Mod(Is)
Is2<-c()
for(i in 1:nI){
	Is2<-c(Is2,rep(Is[i],dupI[i]))
}
Is3<-c(Is2,Conj(Is2))

}
Is3

Ss<-c(Rs2,Is3)

n<-length(Ss)

# 
S<-matrix(0,n,n)
for(i in 1:nR){
	S[i,i]<-Rs[i]
}
for(i in 1:nI){
	st<-(i-1)*2+1+length(Rs2)
	S[st,st]<-Re(Is[i])
	S[st,st+1]<--Im(Is[i])
	S[st+1,st]<-Im(Is[i])
	S[st+1,st+1]<-Re(Is[i])
}
S
library(GPArotation)
V2<-Random.Start(n)
V2<-matrix(rnorm(n^2),n,n)
U2<-solve(V2)

A<-V2%*%S%*%U2

eigen.out<-eigen(A)
eigen.out

V<-eigen.out[[2]]
U<-solve(V)
B<-diag(eigen.out[[1]])

x<-seq(from=0,to=10,length=101)

Ys<-matrix(0,length(x),n)
Ys0<-runif(n)
Ys0<-rep(0,n)
Ys0[1]<-1
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")

nR<-0
Rs2<-c()
if(nR>0){
	Rs<-c(runif(nR)*2)
	dupR<-rep(4,nR)

	Rs2<-c()
	for(i in 1:nR){
		Rs2<-c(Rs2,rep(Rs[i],dupR[i]))
	}

}
nI<-4
Is3<-c()
if(nI>0){
	dupI<-rep(1,nI)
# そのノルム
norI<-rep(1,nI)
norI<-runif(nI)
Is<-complex(real=rnorm(nI),imaginary=rnorm(nI))
#Is<-complex(real=rep(0,nI),imaginary=rnorm(nI))

Is<-Is/Mod(Is)*norI
Mod(Is)
Is2<-c()
for(i in 1:nI){
	Is2<-c(Is2,rep(Is[i],dupI[i]))
}
Is3<-c(Is2,Conj(Is2))

}
Is3

Ss<-c(Rs2,Is3)

n<-length(Ss)

# 
S<-matrix(0,n,n)
for(i in 1:nR){
	S[i,i]<-Rs[i]
}
for(i in 1:nI){
	st<-(i-1)*2+1+length(Rs2)
	S[st,st]<-Re(Is[i])
	S[st,st+1]<--Im(Is[i])
	S[st+1,st]<-Im(Is[i])
	S[st+1,st+1]<-Re(Is[i])
}
S
library(GPArotation)
V2<-Random.Start(n)
V2<-matrix(rnorm(n^2),n,n)
U2<-solve(V2)

A<-V2%*%S%*%U2

eigen.out<-eigen(A)
eigen.out

V<-eigen.out[[2]]
U<-solve(V)
B<-diag(eigen.out[[1]])

x<-seq(from=0,to=10,length=101)

Ys<-matrix(0,length(x),n)
Ys0<-runif(n)
Ys0<-rep(0,n)
Ys0[1]<-1
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")
  • すべて虚数解ですべて実部が0

nR<-0
Rs2<-c()
if(nR>0){
	Rs<-c(runif(nR)*2)
	dupR<-rep(4,nR)

	Rs2<-c()
	for(i in 1:nR){
		Rs2<-c(Rs2,rep(Rs[i],dupR[i]))
	}

}
nI<-4
Is3<-c()
if(nI>0){
	dupI<-rep(1,nI)
# そのノルム
norI<-rep(1,nI)
norI<-runif(nI)
#Is<-complex(real=rnorm(nI),imaginary=rnorm(nI))
Is<-complex(real=rep(0,nI),imaginary=rnorm(nI))

Is<-Is/Mod(Is)*norI
Mod(Is)
Is2<-c()
for(i in 1:nI){
	Is2<-c(Is2,rep(Is[i],dupI[i]))
}
Is3<-c(Is2,Conj(Is2))

}
Is3

Ss<-c(Rs2,Is3)

n<-length(Ss)

# 
S<-matrix(0,n,n)
for(i in 1:nR){
	S[i,i]<-Rs[i]
}
for(i in 1:nI){
	st<-(i-1)*2+1+length(Rs2)
	S[st,st]<-Re(Is[i])
	S[st,st+1]<--Im(Is[i])
	S[st+1,st]<-Im(Is[i])
	S[st+1,st+1]<-Re(Is[i])
}
S
library(GPArotation)
V2<-Random.Start(n)
V2<-matrix(rnorm(n^2),n,n)
U2<-solve(V2)

A<-V2%*%S%*%U2

eigen.out<-eigen(A)
eigen.out

V<-eigen.out[[2]]
U<-solve(V)
B<-diag(eigen.out[[1]])

x<-seq(from=0,to=10,length=101)

Ys<-matrix(0,length(x),n)
Ys0<-runif(n)
Ys0<-rep(0,n)
Ys0[1]<-1
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")
plot(Re(Ys[,1:2]),type="l")
library(rgl)
plot3d(Re(Ys[,1:3]),type="l")