Перейти до вмісту

Чисельні методи. Лабораторний практикум/Проблеми власних значень

Матеріал з Вікіпідручника

Розв'язання повної проблеми власних значень

[ред.]

Постановка задачі

[ред.]

Знайти всі власні числа матриці 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

Посилання

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