#function [lamda_1, v_1, success] = inverse_power_method... # (lamda, q0, A, tol, maxitr) # Accompanying program for the text # # Classical and Modern Numerical Analysis: # Theory, Methods and Practice # by Azmy S. Ackleh, Edward J. Allen, # R. Baker Kearfott, and Padmanabhan Seshaiyer # # (Taylor and Francis / CRC Press, 2009) # # [lamda_1, v_1, success] = inverse_power_method(lamda, q0, A, tol, maxitr) # computes an eigenvalue and eigenvector of the matrix A, according to the # inverse power method as described in Section 5.3 (starting on page 303) # of the text. The eigenvalue is approximated using a Rayleigh quotient # (described on page 300 of the text). Iteration stops when either maxitr # iterations of (5.19) have been reached, or the magnitude of the # difference between successive eigenvalue approximations is less than tol. # 'success' is set to true if the difference between successive # approximations is less than tol upon return, and is set to 'false' # otherwise. In any case, lamda_1 is set to the most recent approximation # to the eigenvalue and v_1 is set to the most recent approximation to the # eigenvector, normalized to infinity norm 1. import numpy as np def inverse_power_method(lamda, q0, A, tol, maxitr): # lambda is a reserved word in Python n = A.shape[0] a = np.linalg.inv(A-lamda*np.eye(n)) # inefficient to actually compute the inverse q_nu = q0 alam = 2*np.linalg.norm(q_nu,np.inf) check = 1 success = False for k in range(maxitr): # k = k+1 Why was this here? alam2 = alam q_nu = np.dot(a,q_nu) q_nu = q_nu/np.linalg.norm(q_nu,np.inf) alam = np.dot( q_nu, np.dot(a,q_nu) ) / np.dot(q_nu,q_nu) check = np.linalg.norm(alam-alam2) if (check < tol): success = True break lamda_1 = lamda + 1/alam v_1 = q_nu print('{:9.0f} {:15.4f} {:15.5f}'.format(k,lamda,lamda_1) ) return lamda_1, v_1, success if __name__ == '__main__': A = np.array( [[1,0,0],[0,.375,0],[0,0,-2]] ) print( inverse_power_method(.3, np.random.rand(3), A, 1.e-7, 100) ) ''' $ py inverse_power_method.py 4 0.3000 0.37500 (0.37500000000216793, array([ 5.68986953e-06, 1.00000000e+00, -1.61848527e-08]), True) '''