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