Lab Session 8 Solution

MA-423 : Matrix Computations Lab







√




1. The roots of a quadratic polynomial p(x) := ax2 +bx+c is given by (−b± b2 − 4ac)/2a. Write a matlab function that implements the above formula to compute the roots. Your function will look like this:







function [x1, x2] = quadroot1( a, b, c)




d = sqrt( b^ 2 - 4 * a* c );




x1 = (-b + d) / (2*a);




x2 = (-b - d) / (2*a);




√




The largest (in magnitude) root of p can be computed as x1 = (−b−sign(b) b2 − 4ac)/2a and the second root x2 can be computed from the identity x1x2 = c/a. Write a matlab function that implements the modified method to compute the roots. Your function will look like this:







function [x1, x2] = quadroot2(a, b, c)




d = sqrt( b^ 2 - 4 * a* c );




x1 = (-b - sign(b) * d) / (2*a);




x2 = c/( a * x1 );




Find the roots of x2 −(107 +10 7)x+1 using quadroot1 and quadroot2. Do you observe any difference? Which method is better and why?

]
2. Consider the matrix A := f1 f1 , where f = √eps and eps is matlab’s machine epsilon.







Let T be the result of 10 basic QR steps (without shift) performed on A. Determine T. Now compute eigenvalues of A using matlab command eig. Can the diagonals of T be considered as eigenvalues of A? Justify your answer. Now use the functions quadroot1 and quadroot2 in Question 1 and compute the eigenvalues of A from its characteristic polynomial. Which method is better?




Consider the matrix given by the matlab command A = gallery(5). Compute A5. What are the the eigenvalues of A? Now compute eigenvalues of A using matlab com-mand eig. What are the eigenvalues? Now plot the eigenvalues with the following commands:



A = gallery(5) e = eig(A) plot(real(e),imag(e),’r*’,0,0,’ko’) axis(.1*[-1 1 -1 1])




axis square




What do you observe? Next, repeat the experiment with a matrix where each element is perturbed by a single roundoff error. The elements of gallery(5) vary over four orders of magnitude, so the correct scaling of the perturbation is obtained with







1






e = eig(A + eps*randn(5,5).*A)




Put this statement, along with the plot and axis commands, on a single line and use the up arrow to repeat the computation several times. You will see that the pentagon flips orientation and that its radius varies between 0.03 and 0.07, but that the computed eigenvalues of the perturbed problems behave pretty much like the computed eigenvalues of the original matrix.




This experiment provides evidence for the fact that the computed eigenvalues are the exact eigenvalues of a matrix A + E where the elements of E are on the order of roundoff error compared to the elements of A. This is the best we can expect to achieve with floating-point computation.




The matlab command [V, D] = eig(A) computes eigenvalues and eigenvectors of A. Type help eig for more information. You can compute condition numbers of the eigenvalues by using the command [V, D, s] = condeig(A). Here V and D contain eigenvectors and eigenvalues of A, respectively, and the vector s contains the condition numbers of the eigenvalues, that is, s(j) is the condition number of the eigenvalue D(j,j). Type help condeig for more information. The condition number of λ is
defined by cond(λ) := ∥x∥ ∥y∥ , where Ax = λx and y A = λy .




|y x|




Consider the Wilkinson’s matrix W. It is a 20-by-20 matrix whose diagonal entries are 20, 19, . . . , 1, supper diagonal (just above diagonals) entries are 20 (fixed for all) and rest of the entries are zero. This matrix can be generated as follows: W = zeros(20); W = diag([20:-1:1])+ diag( 20 * ones(1,19), 1).



What are the eigenvalues of W ? Compute condition number of each of the eigenvalues of W. Now perturb W slightly as follows. Set W 1 := W and W 1(20, 1) := ϵ. For




:= 7.8×10 10, 7.5×10 12, 7.8×10 14, compute eigenvalues of W 1. Do these eigenvalues satisfy the perturbation bounds |λ(W ) − λ(W 1)| ≤ cond(λ)ϵ + O(ϵ2)?



Now compute [V, D] = eig(W) and cond(V ). Do you observe some sort of rela-tionship between cond(V ) and the condition numbers of the eigenvalues of W ? Which eigenvalues are most sensitive to perturbations? (look at the results you have computed above)



Next, for 500 random perturbations Ei with ∥Ei∥ ≤ 10 12, plot (real and imagi-nary parts) of the eigenvalues of W + Ei and W (in a single plot). The distribution of eigenvalues illustrate geometrically the sensitivity of the eigenvalues of W.



The matlab command jordan(A), computes jordan canonical form of a small matrix A with integers entries. Type help jordan for more information. Try to compute jordan canonical forms of W and W 1 considered above. What do you observe?



From all the results above, can you conclude that the distance of W from the set of de-fective matrices is O(10 14)? [Exact distance is 6.13×10 14.] As an illustration, compute eigenvalues of W 1 for ϵ := 10 15. Then W 1 is away from defective matrices and so W should have real and simple eigenvalues as W does. Does your experiment confirm this? Do the eigenvalues of W 1 now satisfy the perturbation bounds given above?










*************End************










2