#非線形固有値問題の解法
#ニュートン法
M=matrix(c(2,-1,-1,2),ncol=2)
C=matrix(c(5,6,7,8),ncol=2)
K=diag(1,2)
u=c(1,-1)
f=function(lambda){
return(M*lambda^2+C*lambda+K)
}
h=0.001
v=rep(0.1,2)
lam=1
a=10
eta=0.01
ite=10000
for(l in 1:ite){
df=(f(lam+h)-f(lam))/h
J=rbind(cbind(f(lam),df%*%v),c(t(u),0))
grad=solve(J)%*%c(f(lam)%*%v,sum(u*v)-a)
lam=lam-eta*grad[length(grad)]
v=v-eta*grad[-length(grad)]
print(sum(abs(f(lam)%*%v)))
}
#逆反復法
M=matrix(c(2,-1,-1,2),ncol=2)
C=matrix(c(5,6,7,8),ncol=2)
K=diag(1,2)
u=c(1,-1)
f=function(lambda){
return(M*lambda^2+C*lambda+K)
}
h=0.001
v=rep(1,2)
v=v/sqrt(sum(v^2))
lam=1
eta=0.001
ite=10000
for(l in 1:ite){
df=(f(lam+h)-f(lam))/h
v2=solve(f(lam))%*%t(f(lam))%*%v
lam=lam-eta*sum(u*v)/sum(u*v2)
v=v*(1-eta)+eta*v2/sqrt(sum(v2^2))
print(sum(abs(f(lam)%*%v)))
}
#QR法による高速計算アルゴリズム(ヘッセンベルグ行列)
A=matrix(c(1,2,3,4,5,6,7,8,9),ncol=3)
AA=A
n=nrow(A)
#ヘッセンベルグ行列の作成:H
#数値計算Ⅰ 朝倉書店 新谷尚義先生 p.124 ハウスホルダー法参照
for(i in 1:(n-2)){
u=c(rep(0,i),A[(i+1):n,i])
sig=sum(A[(i+1):n,i]^2)
if(sig>0){
u[(i+1)]=u[(i+1)]+sign(A[(i+1),i])*sqrt(sig)
Hi=sig+abs(A[(i+1),i])*sqrt(sig)
P=diag(rep(1,n))-t(t(c(u)))%*%t(c(u))/Hi
}else{
P=diag(rep(1,n))
}
A=P%*%A%*%P
}
H=A
t(P)%*%AA%*%P
H
#QR法(ヘッセンベルグ行列による)
A_mat=H
ite=1000
Q=diag(rep(1,nrow(A_mat)))
for(l in 1:ite){
q_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
r_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
q1=A_mat[,1]/sqrt(sum(A_mat[,1]^2))
q_mat[,1]=q1
r_mat[1,1]=sqrt(sum(A_mat[,1]^2))
for(j in 2:ncol(A_mat)){
u=A_mat[,j]-apply(t(t(q_mat[,1:(j-1)])*c(t(q_mat[,1:(j-1)])%*%A_mat[,j])),1,sum)
r_mat[j,j]=sqrt(sum(u^2))
q_mat[,j]=u/sqrt(sum(u^2))
}
r_mat2=t(A_mat)%*%q_mat
r_mat2[upper.tri(r_mat2)]=0
diag(r_mat2)=diag(r_mat)
r_mat=t(r_mat2)
A_mat=r_mat%*%q_mat
Q=Q%*%q_mat
}
#最小多項式による固有ベクトルの計算
#行列HとAAの固有値:lam
lam=sort(unique(diag(A_mat)),decreasing = T)
lambda=diag(A_mat)
count=c()
for(j in 1:length(lam)){
count=c(count,sum(lambda==lam[j]))
}
y=array(0,dim=c(nrow(A),length(lam)))
for(j in 1:length(lam)){
countvec=count
countvec[j]=countvec[j]-1
mat=diag(rep(1,nrow(A)))
for(i in 1:length(countvec)){
if(countvec[i]>0){
for(k in 1:countvec[i]){
mat=mat%*%(A-lam[i]*diag(rep(1,nrow(A))))
}}}
y[,j]=mat[,1]
}
#AAの固有ベクトル
Z=array(0,dim=dim(y))
for(j in 1:n){
Z[,j]=P%*%y[,j]
}
#Householder変換によるQR分解
a=t(matrix(c(1,5,4,2,4,-7,2,7,14),nrow=3))
n=nrow(a)
r=a
q=diag(1,n)
u=rep(0,n)
for(k in 1:(n-1)){
absx=sqrt(sum(r[k:n,k]^2))
u[k]=r[k,k]+sign(r[k,k])*absx
absu=u[k]^2
for(i in (k+1):n){
u[i]=r[i,k]
absu=absu+u[i]^2
}
H=diag(1,n)
for(i in k:n){
for(j in k:n){
H[i,j]=H[i,j]-2*u[i]*u[j]/absu
}}
r=H%*%r
q=q%*%H
}
#Householder変換によるQR分解
QR_decomp=function(A,n){
u=rep(0,n)
r=A
q=diag(1,n)
for(k in 1:(n-1)){
absx=sqrt(sum(r[k:n,k]^2))
u[k]=r[k,k]+sign(r[k,k])*absx
absu=u[k]^2
for(i in (k+1):n){
u[i]=r[i,k]
absu=absu+u[i]^2
}
H=diag(1,n)
for(i in k:n){
for(j in k:n){
H[i,j]=H[i,j]-2*u[i]*u[j]/absu
}}
r=H%*%r
q=q%*%H
}
return(list(r,q))
}
#QR法(Householder変換による)
A_mat=matrix(c(1,2,3,4,5,6,7,8,9),ncol=3)
ite=1000
Q=diag(rep(1,nrow(A_mat)))
for(l in 1:ite){
QRs=QR_decomp(A_mat,nrow(A_mat))
r_mat=QRs[[1]]
q_mat=QRs[[2]]
r_mat2=t(A_mat)%*%q_mat
r_mat2[upper.tri(r_mat2)]=0
diag(r_mat2)=diag(r_mat)
r_mat=t(r_mat2)
A_mat=r_mat%*%q_mat
Q=Q%*%q_mat
}