Чисельні методи. Лабораторний практикум/Проблеми власних значень
Розв'язання повної проблеми власних значень
ред.Постановка задачі
ред.Знайти всі власні числа матриці 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[0]
#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[0])):
print "Eigenvalue: ", res[0][i]," Eigenvector: ", res[1][:,i]
print "| v*A - v*alpha | = ",numpy.linalg.norm(numpy.dot(res[1][:,i],a) - res[1][:,i]*res[0][i])
print