
	Trending ▼ ICSE CBSE 10th ISC CBSE 12th CTET GATE UGC NET Vestibulares ResFinder

numerical analysis

31 pages, 0 questions, 0 questions with responses, 0 total responses,

	anup16	+Fave Message
	Profile Timeline Uploads

Home > anup16 >

Formatting page ...

Computing Eigenvalues and/or Eigenvectors;Part 1, Generalities and symmetric matrices Tom Lyche Centre of Mathematics for Applications, Department of Informatics, University of Oslo November 9, 2008 Today Given a matrix A Cn,n . Finding the eigenvalues using the characteristic polynomial? Perturbation theory Reduction to Hessenberg form Sylvester s inertia theorem Find one or more selected eigenvalues of a symmetric, tridiagonal matrix Find one or more selected eigenvectors (next time) Find all eigenvalues and eigenvectors (next time) Eigenvalues and Characteristic Polynomial The eigenvalues of A Cn,n are the n roots of the characteristic polynomial A ( ) := det(A I) = 0. A ( ) is of exact degree n Except for some special matrices the eigenvalues must be found numerically. Characteristic Polynomial Possible method: Compute the characteristic polynomial A ( ) and apply a numerical method like Newton s method to nd one or more of its roots. Not suitable as an all purpose method. Reason: A small change in one of the coe cients of A ( ) can lead to a large change in the roots of the polynomial Example: A ( :) = 16 . q ( ) = 16 10 16 . Roots of A are all equal to zero. Roots of q are j = 10 1 e 2 ij /16 , j = 1, . . . , 16. The roots of q have absolute value 0.1 Computed roots can be very inaccurate. Need to work directly with the matrix. Gerschgorins circle theorem Where are the eigenvalues? Theorem Suppose A Cn,n . De ne for i , j = 1, 2, . . . , n n Ri = { C : | aii | ri }, ri := |aij |, j =1 j =i n Cj = {z C : |z ajj | cj }, cj := |aij |. i =1 i =j Then any eigenvalue of A lies in R C where R = R1 R2 Rn and C = C1 C2 Cn . If AH = A then Ci = Ri = [aii ri , aii + ri ]. Examples Locate eigenvalues for A = 2 1 1 2 A is symmetric. R1 = R2 = [2 1, 2 + 1] = [1, 3] = R , so [1, 3]. 1 = 3 and 2 = 1. In this case the smallest interval possible. Let T = tridiag( 1, 2, 1) Rm,m be the second derivative matrix. R1 = Rm = [1, 3], Ri = [0, 4], i = 2, 3, . . . , m 1, so R = [0, 4]. j j = 4 sin 2(m +1) 2 , j = 1, 2, . . . , m. j [ , 4 ], where = 1/(2(m + 1)) Gerschgorins theorem gives a remarkably good estimate. proof of Gerschgorin Proof. Suppose ( , x) is an eigenpair for A. We claim that Ri , where i is such that |xi | = x . Indeed, Ax = x implies that j aij xj = xi or ( aii )xi = j =i aij xj . Dividing by xi and taking absolute values we nd | aii | = | j =i aij xj /xi | j =i |aij ||xj /xi | ri Thus Ri . Since is also an eigenvalue of AT , it must be in one of the row disks of AT . But these are the column disks Cj of A. Hence Cj for some j . Distinct circles Sometimes some of the Gerschgorin disks are distinct and we have Corollary If p of the Gerschgorin row disks are distinct from the others, the union of these disks contains precisely p eigenvalues. The same result holds for the column disks. 112 A = 3 2 4 where | j | 10 10 3 5 6 | j j | 2 10 10 for j = 1, 2, 3 Perturbation Analysis Recall linear systems Ax = b and Ay = b + e e y x p Kp (A) b p , where Kp (A) := A xp p p A 1 p . This means that the relative error y xp p in y as an x approximation to x can possibly be Kp (A) as large as the e relative error b p in the right hand side b. p Consider now the eigenvalue problem. We restrict the discussion to nondefective matrices. Absolute errors Theorem Suppose A Cn,n has linearly independent eigenvectors {x1 , . . . , xn } and let X = [x1 , . . . , xn ] be the eigenvector matrix. If ( , x) is an eigenpair for A + E, then we can nd an eigenvalue of A such that | | Kp (X) E p , 1 p . (1) If A is symmetric then | | E 2 . (2) Two observations It is di cult or sometimes impossible to compute accurate eigenvalues and eigenvectors of matrices with linearly dependent eigenvectors or almost linearly dependent eigenvectors. The eigenvalue problem for symmetric matrices is well conditioned. Upper Hessenberg Before attempting to nd eigenvalues and eigenvectors of a matrix (exceptions are made for certain sparse matrices), it is often advantageous to reduce it by similarity transformations to a simpler form. Orthogonal similarity transformations are particularly important since they are insensitive to noise in the elements of the matrix. Recall that a matrix A Rn,n is upper Hessenberg if ai ,j = 0 for j = 1, 2, . . . , i 2, i = 3, 4, . . . , n. x x 0 0 0 0 x x x 0 0 0 x x x x 0 0 x x x x x 0 x x x x x x x x x x x x A matrix H Rn,n of the form H := I uuT , where u Rn and uT u = 2 is called a Householder transformation. Zero out entries in a vector x Find a Householder transformation H := I uuT such that Hx = e1 . u := x / e1 1 x1 / 2e1 , if x = 0, , := otherwise. x +x 2 2 if x1 > 0 otherwise, H = diag( 1, 1, . . . , 1) if = 0. Assume = 0. uT u = (x/ e1 )T (x/ e1 ) 1 x1 / T = ( 2 / 2x1 / +1 1 x1 / =2 Hx = x (u x)u x/ uT x = ( e1 ) Tx 1 x1 / T x = x/ x1 = x1 1 x1 / 1 x1 / = Hx = x (uT x)u = x (x/ e1 ) = e1 . 1 x1 / . Computing u u := x / e1 1 x1 / 2e1 , = x 2 If = 0 then u = 2e1 exit v = x/ e1 u = v/ v (1) if x = 0, otherwise. Recall Algorithm housegen To given x Rn the following algorithm computes a = and the vector u so that (I uuT )x = e1 . Algorithm function [u,a]=housegen(x) a=norm(x); u=x; if a==0 u(1)=sqrt(2); return; end if u(1)>0 a=-a; end u=u/a; u(1)=u(1)-1; u=u/sqrt(-u(1)); Reduction to upper Hessenberg we de ne A1 = A. Suppose for k 1 that Ak is upper Hessenberg in its rst k 1 columns. Bk Ck Ak = Dk Ek Bk Rk ,k is upper Hessenberg and Dk = [0, 0, . . . , 0, dk ] Rn k ,k . I0 Let Hk = 0 Vk , where Vk = I vk vT Rn k ,n k is a k Householder transformation such that Vk dk = k e1 , 2 where k = dT dk . k Reduction 2 Ak +1 = Hk Ak Hk = = Ik 0 0 Vk Bk Ck Vk Vk Dk Vk Ek Vk Bk Ck Dk Ek Ik 0 0 Vk . Now Vk Dk = [Vk 0, . . . , Vk 0, Vk dk ] = (0, . . . , 0, k e1 ) and the matrix Bk is not a ected by the Hk transformation. Therefore the matrix Ak +1 will be upper Hessenberg in its rst k columns, and the reduction is carried one step further. The reduction stops with An 1 which is upper Hessenberg. O (10n3/3) algorithm Algorithm function [L,B] = hesshousegen(A) n=length(A); L=zeros(n,n); B=A; for k=1:n-2 [v,B(k+1,k)]=housegen(B(k+1:n,k)); L(k+1:n,k)=v; B(k+2:n,k)=zeros(n-k-1,1); C=B(k+1:n,k+1:n); B(k+1:n,k+1:n)=C-v*(v *C); C=B(1:n,k+1:n); B(1:n,k+1:n)=C-(C*v)*v ; end This algorithm uses Householder similarity transformations to reduce a matrix A Rn,n to upper Hessenberg form B. Details of the transformations are stored under the diagonal in a matrix L. The entries of L can be used to assemble an orthogonal matrix Q such that B = QT AQ. Algorithm housegen is used in each step of the reduction. The symmetric case If A1 = A is symmetric, the matrix An 1 will also be symmetric since AT = Ak implies k AT+1 = (Hk Ak Hk )T = Hk AT Hk = Ak +1 . k k Since An 1 is upper Hessenberg and symmetric, it must be tridiagonal. Thus the algorithm above reduces A to symmetric tridiagonal form if A is symmetric. Symmetric tridiagonal AT = A Rn,n eigenvalues 1 2 n . Using Householder similarity transformations we can assume that A is symmetric and tridiagonal. d1 c1 c1 d2 c2 .. .. .. A= . . . . cn 2 dn 1 cn 1 cn 1 dn (3) Split tridiagonal A into irreducible components Recall that A is reducible if ci = 0 for at least one i . Example: Suppose n = 4 and c2 = 0 d1 c1 A= 0 0 c1 d2 0 0 0 0 d3 c3 0 0 = A1 0 . 0 A2 c3 d4 The eigenvalues of A are the union of the eigenvalues of A1 and A2 . Thus if A is reducible then the eigenvalue problem can be split into smaller irreducible problems. So assume that A is irreducible. Theorem: An irreducible, symmetric, tridiagonal matrix has distinct eigenvalues. The inertia theorem We say that two matrices A, B Cn,n are congruent if A = EH BE for some nonsingular matrix E Cn,n . Let (A), (A) and (A) denote the number of positive, zero and negative eigenvalues of A. If A is Hermitian then (A) + (A) + (A) = n. Theorem (Sylvester s Inertia Theorem) If A, B Cn,n are Hermitian and congruent then (A) = (B), (A) = (B) and (A) = (B). LDLT factorization If A = LT DL is an LDLT-factorization of A then A and D are congruent. (A) = (D), (A) = (D) and (A) = (D). 13 10 = 34 31 10 0 5 13 01 one positive and one negative eigenvalue. Corollary Suppose for some x R that A x I has an LDLT-factorization A x I = LDLT . Then the number of eigenvalues of A strictly less than x equals the number of negative diagonal entries in D. (A I) = (D) If Ax = x then (A I)x = ( )x, (A I) equals the number of eigenvalues of A which are less than . Counting eigenvalues in an interval Suppose AT = A Rn,n Using for example Gerschgorin s theorem we can nd an interval [a, b ) containing the eigenvalues of A. For x [a, b ) let (x ) be the number of negative diagonal entries in D in an LDLT-factorization of A x I. (x ) is the number of eigenvalues of A which are strictly less than x . (a) = 0, (b ) = n (e ) (d ) is the number of eigenvalues in [d , e ). Approximating m 1 2 n . Suppose 1 m n. Find m using interval bisection. Let c = (a + b )/2 and k := (c ). If k m then m c and m [a, c ], while if k < m then m c and m [c , b ]. Continuing with the interval containing m we generate a sequence {[aj , bj ]} of intervals, each containing m and bj aj = 2 j (b a). Fixing a possible failure The method will fail if one of the diagonal entries in D is zero or very close to zero We replace such an entry by a suitable small number, say k = |ck | M , where the negative sign is used if ck < 0, and M is the Machine epsilon, typically M = 2 10 16 for Matlab. This replacement is done if |dk ( )| < | k |. function k=count(c,d,x) Algorithm Suppose A = tridiag(c, d, c) is symmetric and tridiagonal with entries d1 , . . . , dn on the diagonal and c1 , . . . , cn 1 on the neighboring diagonals. For given x this function counts the number of eigenvalues of A strictly less than x. We use the replacement described above if one of the dj (x ) is close to zero. function lambda= ndeigv(c,d,m) Algorithm Suppose A = tridiag(c, d, c) is symmetric and tridiagonal with entries d1 , . . . , dn on the diagonal and c1 , . . . , cn 1 on the neighboring diagonals. We rst estimates an interval [a, b ] containing all eigenvalues of A and then generates a sequence {[ak , bk ]} of intervals each containing m . We iterate until bk ak (b a) M , where M is Matlab s machine epsilon eps. Typically M 2.22 10 16 . Example Given T := tridiag( 1, 2, 1) of size 100. Estimate l5 5 . Using findeigv we nd l5 = 0.024139120518486 Using Matlab s eig we nd 5 = 0.024139120518486 Which is most accurate? Exact 5 = 4 sin ( 5/202)2 = 0.024139120518487 5 5 = 8.6e 016 not bad! l5 5 = 3.4e 016 better!

Formatting page ...

anup16 chat

Horizontal lines at:
Guest Horizontal lines at:
AutoRM Data:
Box geometries: