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