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]] );
#
#
#
#