# Assig5part2

unknown
python
a year ago
1.8 kB
2
Indexable
Never
```import numpy as np
import scipy
from scipy import linalg

def eig(A):
k = 0
#Use the householder reduction to hessenberg form
Ak = tridiagonalization(A)
m, n = np.shape(Ak)
#If it is a 1x1 matrix return
if(m < 2):
return Ak
#If the matrix is 2x2 run the algorithm for 5 times (gives good results)
#If it is bigger run untill we get 0 element in sub/super diag
max = 10000
if m == 2:
max = 5
while k < max:
#Shift
mu = Ak[m-1][m-1]
#QR factorisation
Q, R = scipy.linalg.qr(Ak-mu*np.eye(m))
#Reconstruction of the matrix for the next step.
Ak = R@Q+mu*np.eye(m)
#Test if the matrix has zero element on the sub/super diag
#If so call the function recursively
if m > 2:
for i in range(1,m-1):
if np.abs(Ak[i-1][i]) < 10**(-9):
A1 = eig(Ak[0:i,0:i])
A2 = eig(Ak[i:m,i:m])
return np.block([[A1, np.zeros((i, m-i))],[np.zeros((m-i,i)), A2]])
k += 1
return Ak

def tridiagonalization(A):
m,n = np.shape(A)
for k in range(0,m-2):
x = A[k + 1 : m, k]
size = np.size(x)
e1 = np.eye(1,size,0)
vk = np.sign(x[0])*scipy.linalg.norm(x)*e1 + x
vk = vk/scipy.linalg.norm(vk)
A[k+1:m,k:m] = A[k+1:m,k:m] - 2*vk.T@(vk@A[k+1:m,k:m])
A[0:m,k+1:m] = A[0:m,k+1:m] - 2*(A[0:m,k+1:m]@vk.T)@vk
return A

#Dim of A
n = 100
#Create a random symetric matrix
A = np.random.rand(n,n)
A = 1/2*(A+A.T)
#Call our eig to calculator
eigs = np.diag(eig(A))
eigs = np.sort(eigs)
print(eigs)
#Use scipy eig to calculate
a,_ = scipy.linalg.eig(A)
a = np.real(a)
a = np.sort(a)
print(a)
print(scipy.linalg.norm(a-eigs))```