interface( warnlevel=0 ): with( laylinalg ): QR := proc( A::Matrix ) local q, r; q,r := QRDecomposition( evalf(A), output=['Q','R'] ); return r . q end proc: # Text Reference: Section 6.4 # Lab 7: Polynomial Roots via the QR-Method for Eigenvalues # The purpose of this set of exercises is to show how the real roots of a polynomial can be calculated by finding the eigenvalues of a particular matrix. These eigenvalues will be found by the QR method described below. # # A polynomial of degree n; is a function of the form # p(t) = a[0]+a[1]*t+`. . .`+a[n-1]*t^(n-1)+a[n]*t^n; # where a[0];, a[1];, . . ., a[n-1];, and a[n]; are real numbers with a[n] <> 0;. A root of a polynomial is a value of t; for which p(t) = 0;. It is often necessary (especially in calculus-based applications) to find all of the real roots of a given polynomial. In practice this can be a difficult problem even for a polynomial of low degree. For a polynomial of degree 2, every algebra student learns that the roots of a*t^2+b*t+c; can be found by the quadratic formula # t = -b*`+-`*sqrt(-4*a*c+b^2)/(2*a); # If the polynomial is of degree 3 or 4, then there are formulas somewhat resembling the quadratic formula (but much more involved) for finding all the roots of a polynomial. However there is no general formula for finding the roots of a polynomial of degree 5 or higher. # # Example: # Consider the monic cubic polynomial p(t) = t^3-2*t^2-5*t+6;(monic means the leading coefficient is 1). This polynomial is factored rather easily to find that its roots are t = 1;, t = -2;, and t = 3;. # # Polynomial Roots using Linear Algebra # # If a polynomial cannot easily be factored, numerical techniques are used to find a polynomial's roots. There are problems with this approach as well. Algorithms such as Newton's Method may not converge to a root, or may approach the root very slowly. These methods must also be applied repeatedly to find all of the roots, and usually require a cleverly chosen starting guess for the root being sought. However, there is an algorithm from linear algebra which may be used to find all real roots of a polynomial simultaneously. # # The eigenvalues of an n; x n; matrix A are the roots of the characteristic polynomial of A, which is defined as q(lambda)=det(A-lambda I[n]).; This # polynomial of degree n, ; because A is n; x n;. So to know the eigenvalues of A is to know the roots of the monic polynomial qlambda;. # # To find the roots of any given monic polynomial p(t);, then, two problems need to be solved: # # 1. A way to construct a square matrix A whose characteristic polynomial q(lambda); equals p(lambda);. # 2. A way to find the eigenvalues of this matrix A which does not depend on finding the roots of q(lambda);. # # The first problem is solved by defining the companion matrix for a (monic) polynomial # p(t) = a[0]+a[1]*t+`. . .`+a[n-1]*t^(n-1)+t^n; # # Definition: # If p(t) = a[0]+a[1]*t+`. . .`+a[n-1]*t^(n-1)+t^n;, then the companion matrix for p is # C[p]; = [[[0,1,0,`...`,0],[0,0,1,`...`,0],[`...`,`...`,`...`,^`.` `.` [`.`],`...`],[0,0,0,`...`,1],[-a[0],-a[1],-a[2],`...`,-a[n-1]]]];. # # Example: Cp := Matrix( 3,3, [[0,1,0],[0,0,1],[-6,5,2]] ): Cp; # The companion matrix for the polynomial p(t) = t^3-2*t^2-5*t+6; is C[p]; = Matrix(3, 3, {(1, 1) = 0, (1, 2) = 1, (1, 3) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (3, 1) = -6, (3, 2) = 5, (3, 3) = 2});. # # Problems to be Submitted: # Problem 1. # Find the companion matrices for the following polynomials. # (a) # p(t) = t^2-6*t+8; # (b) # p(t) = t^3-3*t^2-10*t+24; # (c) # p(t) = t^4-2*t^3-13*t^2+14*t+24; # (d) # Find the characteristic polynomials of the matrices you just found in parts (a)-(c). A Maple command such as solve(20-10*t-3*t^2+t^3=0,t) finds the roots of each polynomial. The Maple command CharacteristicPolynomial(M,lambda) finds the characteristic polynomial from the matrix M. # # Problem 2. # Show that the characteristic polynomial of a companion matrix for the n^th; degree polynomial p(t); is det( C[p]; - I[n]; ) = (-1)^n*p(lambda); as follows. # (a) # Show that if C[p]; is the companion matrix for a quadratic polynomial p(t) = t^2+t*a[1]+a[0];, then det( C[p]; - I[2]; ) = (-1)^2;p(lambda); by direct computation. # # (b) # Use mathematical induction to show that the result holds for n; >= 2. # Hint: # Expand the necessary determinant by cofactors down the first column. # # The QR Method for Eigenvalues # # The companion matrix is a matrix A whose characteristic polynomial is p(t);. A method for finding the eigenvalues of A which does not use the characteristic polynomial is also needed. One method which accomplishes this is called the QR method because it is based on the QR factorization of A. # # The QR factorization of an m; x n; matrix A;requires the matrix to have linearly independent columns. Then A can be factored as A = QR , where Q is an m; x n; matrix with orthonormal columns and R is an n; x n; invertible upper triangular matrix with positive entries on its main diagonal. # Problem 3. # Suppose A is a n; x n; matrix. Let A = Q[0]; R[0]; be a QR factorization of A, and create A[1]; = R[0]; Q[0];. Let A[1]; = Q[1]; R[1]; be a QR factorization of A[1]; and create A[2]; = R[1]; Q[1];. # (a) # Show that A = Q[0]; A[1]; Q[0]^T;. (This is Exercise 23, Section 5.2.) # # (b) # Show that A = ( Q[0]; Q[1]; ) A[2]; ( Q[0]; Q[1]; )^T;. # # (c) # Show that Q[0]; Q[1]; is an orthogonal matrix. (This is Exercise 29, Section 6.2.) # # (d) # Show that A, A[1];, and A[2]; all have the same eigenvalues. # # # The QR Method # The QR method for finding the eigenvalues of an n; x n; matrix A extends the process in Problem 3 to create a sequence of matrices with the same eigenvalues. # # Step 1: # Let A = Q[0]; R[0]; be a QR factorization of A; create A[1]; = R[0]; Q[0];. # # Step 2: # Let A[1]; = Q[1]; R[1]; be a QR factorization of A[1];; create A[2]; = R[1]; Q[1];. # # Step m+1: # Continue this process. Once A[m]; has been created, then let A[m]; = Q[m]; R[m]; be a QR factorization of A[m]; and create A[m+1]; = R[m]; Q[m];. # Stopping Criterion: # Stop the process when the entries below the main diagonal of A[m]; are sufficiently small, or stop if it appears that convergence will not happen. # QRDecomposition( Cp ); # Example a0 := evalf( Cp ): for i from 0 to 15 do q,r := QRDecomposition( a||i, output=['Q','R'] ); a||(i+1) := r . q; end do: q0,r0 := QRDecomposition( a0, output=['Q','R'] ); a0 = q0 . r0; q1,r1 := QRDecomposition( a1, output=['Q','R'] ); a1 = q1 . r1; a2; a13; # Let A be the companion matrix for the monic cubic polynomial p(t) = t^3-2*t^2-5*t+6;; that is, A = Matrix(3, 3, {(1, 1) = 0., (1, 2) = 1., (1, 3) = 0., (2, 1) = 0., (2, 2) = 0., (2, 3) = 1., (3, 1) = -6., (3, 2) = 5., (3, 3) = 2.});. # The QR factorization of this matrix is # A = Q[0]; R[0]; = Matrix(3, 3, {(1, 1) = 0., (1, 2) = 1., (1, 3) = 0., (2, 1) = 0., (2, 2) = 0., (2, 3) = 1., (3, 1) = -1., (3, 2) = 0., (3, 3) = 0.}); Matrix(3, 3, {(1, 1) = 6., (1, 2) = -5., (1, 3) = -2., (2, 1) = 0., (2, 2) = 1., (2, 3) = 0., (3, 1) = 0., (3, 2) = 0., (3, 3) = 1.}); , so A[1]; = R[0]; Q[0]; = Matrix(3, 3, {(1, 1) = 2., (1, 2) = 6., (1, 3) = -5., (2, 1) = 0., (2, 2) = 0., (2, 3) = 1., (3, 1) = -1., (3, 2) = 0., (3, 3) = 0.});. # This operation is performed again, producing # A[2]; = R[1]; Q[1]; = Matrix(3, 3, {(1, 1) = 2.236067977, (1, 2) = 5.366563146, (1, 3) = -4.472135955, (2, 1) = 0., (2, 2) = 2.683281573, (2, 3) = -2.236067977, (3, 1) = 0., (3, 2) = 0., (3, 3) = 1.}); Matrix(3, 3, {(1, 1) = .8944271910, (1, 2) = .4472135955, (1, 3) = 0., (2, 1) = 0., (2, 2) = 0., (2, 3) = 1., (3, 1) = -.4472135955, (3, 2) = .8944271910, (3, 3) = 0.}); = Matrix(3, 3, {(1, 1) = 4.000000000, (1, 2) = -3.000000000, (1, 3) = 5.366563146, (2, 1) = .9999999998, (2, 2) = -2.000000000, (2, 3) = 2.683281573, (3, 1) = -.4472135955, (3, 2) = .8944271910, (3, 3) = 0.});. # This matrix is still far from upper triangular, so the process is continued. After 13 steps it is found that # A[13]; = Matrix(3, 3, {(1, 1) = 2.991249185, (1, 2) = 5.038287273, (1, 3) = 5.086237669, (2, 1) = 0.8671583910e-2, (2, 2) = -1.991447448, (2, 3) = -1.353492217, (3, 1) = -0.1916695814e-5, (3, 2) = 0.4401732170e-3, (3, 3) = 1.000198260});, # so the matrix is converging to an upper triangular matrix, and its diagonal elements are converging to the roots of p(t);=0, which are t = 3;, t = -2;, and t = 1;. # Maple Implementation of the QR Method QR := proc( A::Matrix ) local q, r; q,r := QRDecomposition( evalf(A), output=['Q','R'] ); return r . q end proc: # The following Maple procedure performs one iteration of the QR Method for a matrix A. # # QR := proc( A::Matrix ) # local q, r; # q,r := LinearAlgebra[QRDecomposition](evalf(A)); # return r . q # end proc: # # The first 13 iterates of the QR Method for the companion matrix for the monic cubic polynomial p(t) = t^3-2*t^2-5*t+6; can be obtained with the Maple commands: aa0 := evalf( Matrix( 3, 3, [[0,1,0],[0,0,1],[-6,5,2]] ) ); for i from 0 to 14 do aa||(i+1) := QR( aa||i ); end do; # # A0 := Matrix([[0,1,0],[0,0,1],[-6,5,2]]); # # Each iterate B will have the same eigenvalues as A0. # B:=A0:for i from 0 to 12 do # B := QR(B); # end do; # # I is the first iterate for which all entries below the diagonal are smaller than 0.01. # Decisions about an intermediate result are aided by the extra function below: # # interface(displayprecision=5);nn:=5; # F:=x->if(abs(x)<1/10^nn) then 0 else x fi; # # then for B=iterate 13, map(F,B) will display Matrix(3, 3, {(1, 1) = 2.99125, (1, 2) = -5.03829, (1, 3) = -5.08624, (2, 1) = -0.867e-2, (2, 2) = -1.99145, (2, 3) = -1.35349, (3, 1) = 0, (3, 2) = 0.44e-3, (3, 3) = 1.00020});, which has 5 display digits and numbers smaller than 1/10^5; are replaced by zero. The approximate eigenvalues of A0 are the approximate eigenvalues of B = iterate 13, which are the diagonal entries 2.99125,-1.99145,1.00020. ;The diagonal entries are printed using Maple code seq(B[j,j],j=1..3); # # # Problem 5. # Find the approximate roots according to the QR method for the following polynomials. Compare the answers using LinearAlgebra[Eigenvalues](A) where A is the companion matrix for the given polynomial. # # The companion matrix for polynomial t^2+2*t-4; can be found from maple code LinearAlgebra[CompanionMatrix(]t^2+2*t-4)^+, but the transpose operation can be ignored, due to determinant rule det(C)=det(C^(T); ), applied with C=A-; I*lambda;. # # (a) # p(t) = t^2-6*t+8; # (b) # p(t) = t^3-3*t^2+10*t+24; # (c) # p(t) = t^4-2*t^3-13*t^2+14*t+24; # Problem 6. # Define A = Matrix(6, 6, {(1, 1) = 0, (1, 2) = 0, (1, 3) = -1, (1, 4) = 4, (1, 5) = -1, (1, 6) = -6, (2, 1) = 0, (2, 2) = -2, (2, 3) = 2, (2, 4) = -5, (2, 5) = -2, (2, 6) = -5, (3, 1) = -1, (3, 2) = 2, (3, 3) = 8, (3, 4) = -4, (3, 5) = 3, (3, 6) = 2, (4, 1) = 4, (4, 2) = -5, (4, 3) = -4, (4, 4) = -6, (4, 5) = 1, (4, 6) = 0, (5, 1) = -1, (5, 2) = -2, (5, 3) = 3, (5, 4) = 1, (5, 5) = -2, (5, 6) = 7, (6, 1) = -6, (6, 2) = -5, (6, 3) = 2, (6, 4) = 0, (6, 5) = 7, (6, 6) = 10}); # # A := Matrix( # [[ 0, 0,-1, 4,-1,-6], # [ 0,-2, 2,-5,-2,-5], # [-1, 2, 8,-4, 3, 2], # [ 4,-5,-4,-6, 1, 0], # [-1,-2, 3, 1,-2, 7], # [-6,-5, 2, 0, 7,10]] ); # # Use Maple and the QR Method to make all the entries below the main diagonal of A less than 0.01. Record how many steps it takes to get this result, and then record your estimates for the eigenvalues of A. # # Answer: (1) About 100 steps. (2) Same as LinearAlgebra[Eigenvalues](A), to 4 digits. The basic code to use is B:=A: for i from 0 to 200 do B := QR( B ): if(normF(B)<1/10^2) then break; fi; end do: printf("i=%d, normF(B)=%f\n",i,normF(B));FA(B); B:=A: for i from 0 to 200 do B := QR( B ): if(normF(B)<1/10^2) then break; fi; end do: printf("i=%d, normF(B)=%f\n",i,normF(B));FA(B); B:=A: for i from 0 to 200 do B := QR( B ): if(normF(B)<1/10^2) then break; fi; end do: printf("i=%d, normF(B)=%f\n",i,normF(B));FA(B); B:=A: for i from 0 to 200 do B := QR( B ): if(normF(B)<1/10^2) then break; fi; end do: printf("i=%d, normF(B)=%f\n",i,normF(B));FA(B); B:=A: for i from 0 to 200 do B := QR( B ): if(normF(B)<1/10^2) then break; fi; end do: printf("i=%d, normF(B)=%f\n",i,normF(B));FA(B); # # Compute theta maximum entry A[i,j] below the main diagonal normF:=proc(A::Matrix) local x,i,j,n; n:=LinearAlgebra[RowDimension](A); x:=max(seq(seq(abs(A[i,j]),i=j+1..n),j=1..n-1)); RETURN ( x); end proc: # # Compute the maximum entry |A[i,j]| below the main diagonal # normF:=proc(A::Matrix) local x,i,j,n; # n:=LinearAlgebra[RowDimension](A); # x:=max(seq(seq(abs(A[i,j]),i=j+1..n),j=1..n-1)); # RETURN (x); # end proc: # # interface(displayprecision=5);nn:=2; # F:=x->if(abs(x)<1/10^nn) then 0 else x fi; # # B:=A:for i from 0 to 200 do # Compute theta maximum entry A[i,j] below the main diagonal # B := QR( B ): # if(normF(B)<1/10^nn) then break; fi # end do; # printf("i=%d, normF(B)=%f\n",i,normF(B)); # Print iterate number and error estimate # map(F,B); # Display B, but print zeros instead of small decimals # # Problem 7. # Listed below are matrices C, D and E and the corresponding Maple command which creates them. For each given matrix, do enough steps of the QR method to find a matrix B with the same eigenvalues having each entry below the main diagonal of matrix B smaller than 0.1. Record the number of steps, the final result, and give estimates for the eigenvalues of each matrix. # # (a) Bmat := Matrix( 3,3, [[1, -2, 8], [7, -7, 6], [5, 7, -8]] ): Bmat; # C = Matrix(3, 3, {(1, 1) = 1, (1, 2) = -2, (1, 3) = 8, (2, 1) = 7, (2, 2) = -7, (2, 3) = 6, (3, 1) = 5, (3, 2) = 7, (3, 3) = -8}); # CC := Matrix([[1, -2, 8], [7, -7, 6], [5, 7, -8]]); # # (b) Cmat := Matrix( 4,4, [[4, -2, 3, -7], [1, 2, 6, 8], [8, 5, 1, -5], [-5, 8, -5, 3]] ): Cmat; # D = Matrix(4, 4, {(1, 1) = 4, (1, 2) = -2, (1, 3) = 3, (1, 4) = -7, (2, 1) = 1, (2, 2) = 2, (2, 3) = 6, (2, 4) = 8, (3, 1) = 8, (3, 2) = 5, (3, 3) = 1, (3, 4) = -5, (4, 1) = -5, (4, 2) = 8, (4, 3) = -5, (4, 4) = 3}); # DD := Matrix([[4, -2, 3, -7], [1, 2, 6, 8], [8, 5, 1, -5], [-5, 8, -5, 3]]); # # (c) Dmat := Matrix( 5,5, [[2, 6, -3, 4, -9], [-1, 7, -4, -3, -7], [-6, -6, -1, 6, 5], [9, 2, 6, 2, -8], [-7, -8, 6, -9, -1]] ): Dmat; # E = Matrix(5, 5, {(1, 1) = 2, (1, 2) = 6, (1, 3) = -3, (1, 4) = 4, (1, 5) = -9, (2, 1) = -1, (2, 2) = 7, (2, 3) = -4, (2, 4) = -3, (2, 5) = -7, (3, 1) = -6, (3, 2) = -6, (3, 3) = -1, (3, 4) = 6, (3, 5) = 5, (4, 1) = 9, (4, 2) = 2, (4, 3) = 6, (4, 4) = 2, (4, 5) = -8, (5, 1) = -7, (5, 2) = -8, (5, 3) = 6, (5, 4) = -9, (5, 5) = -1}); # EE := Matrix([[ 2, 6, -3, 4, -9],[-1, 7, -4, -3, -7],[-6, -6, -1, 6, 5], # [ 9, 2, 6, 2, -8],[-7, -8, 6, -9, -1]] ); # # # #