Чисельні методи. Лабораторний практикум Розв'язання повної проблеми власних значень Постановка задачі Знайти всі власні числа матриці A та власні вектори методом обертання Якобі, де

Аналіз Читаємо статтю: Метод обертання Якобі, там все непогано описано.

Код Писати коментарі в коді краще англійською, менше проблем з кодуванням, і взагалі так прийнято. Короткий словничок: eigenvalue - власне значення eigenvector - власний вектор condition number - число обумовленості identity matrix - одинична матриця. (Відповідає перетворенню ідентичності).

from __future__ import division import numpy

def jacobi(M): """ Determines eigenvalues and eigenvectors for a matrix @param M: matrix to determine eigenvalues and eigenvectors for (numpy array) @return [E, V, s]: numpy array of eigenvalues and numpy array of eigenvectors, as well as number of sweeps used """

s = 0 # number of iterations n = numpy.size(M, 0) # size of matrix eps = 1E-12 # precision

#make identity matrix # here we will store eigenvectors V = numpy.identity(n, float) M = M.astype('float') #rotate until accuracy is reached rotated = True while(rotated): rotated = False for r in range(n): for c in range(r+1, n): #if(numpy.abs(M[c,r]) > eps*(numpy.abs(M[c,c])+numpy.abs(M[r,r]))): #it means if nondiagonal elements aren't far smaller than diagonal we must continue our rotations. if eps*numpy.abs(M[r,r])<numpy.abs(M[r,c]) or eps*numpy.abs(M[c,c])<numpy.abs(M[r,c]): rotated = True rotate(r, c, M, V) s += 1 return [numpy.diag(M), V, s]

def rotate(p, q, A, V): #p is always the lesser one if q < p: p,q = q,p

#calculate "rotations angle" and corresponding coefficients phi = 0.5*numpy.arctan2(2*A[p, q], A[q, q]-A[p, p]) s = numpy.sin(phi) c = numpy.cos(phi)

#work the diagonal app = A[p,p] aqq = A[q,q] apq = A[p,q] A[p,p] = c*c*app+s*s*aqq-2*s*c*apq A[q,q] = s*s*app+c*c*aqq+2*s*c*apq A[p,q] = 0

n = A.shape

#do rotation stuff for i in range(p): aip = A[i,p] aiq = A[i,q] A[i,p] = c*aip-s*aiq A[i,q] = s*aip+c*aiq

for i in range(p+1, q): api = A[p,i] aiq = A[i,q] A[p,i] = c*api-s*aiq A[i,q] = s*api+c*aiq

for i in range(q+1, n): api = A[p,i] aqi = A[q,i] A[p,i] = c*api-s*aqi A[q,i] = s*api+c*aqi

for i in range(n): egvpi = V[i,p] egvqi = V[i,q] V[i,p] = c*egvpi-s*egvqi V[i,q] = s*egvpi+c*egvqi

return

def genarray(size): # create our task data a=numpy.zeros((size,size)) for i in range(size): for j in range(i+1,size): a[i,j]=a[j,i]=size-j a[i][i]=size return a

if __name__ == "__main__": size=input("Input size:") a=genarray(size) print a print "Condition number: ", numpy.linalg.cond(a) res=jacobi(a) for i in range(len(res)): print "Eigenvalue: ", res[i]," Eigenvector: ", res[:,i] print "| v*A - v*alpha | = ",numpy.linalg.norm(numpy.dot(res[:,i],a) - res[:,i]*res[i]) print

Посилання Функція NumPy для обчислення власних значень та векторів матриці NumPy для обчислення числа обумовленості