固有値計算と特異値分解　丸善出版　一般社団法人日本計算工学会編ー計算力学レクチャーコース

```r

#非線形固有値問題の解法
#ニュートン法
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)))
}


```

```r

#逆反復法

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)))
  
}

```



```r


#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]
}



```

```r

#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
}

```

```r
#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
}

```