Homework 7. Assigned: Tue Apr 17. Due: Fri Apr 27 at 5 p.m.

1. QR factorization
   • Code the function qrrot listed in Van Loan's book, p.251, to compute the QR factorization of an m by n matrix A, where you can assume that m &ge n.
   • Compare your output with that returned by [Q,R] = qr(A). Is there any significant difference, and if so, does it make sense?
   • Add the feature that if there is a second input argument, 0 (test the value of nargin), the function returns the "economy-size" version of the QR factorization, just as qr does.
2. Least Squares.
   • Load this x, y data into Matlab and plot the points as asterisks. There are 41 points, so the interpolating polynomial would have degree 40, which will pass through all the points but oscillate wildly in between. A better idea is to approximate the data with a polynomial of much smaller degree. Experiment with several possible choices of degree d, all less than 10. To do this, you can use vander to set up the Vandermonde matrix, from which you set A to consist of the final d+1 columns (or first d+1 columns if you fliplr the Vandermonde matrix first). Then you can use A\y to automatically solve the least squares problem, which is an overdetermined system of 41 equations in d+1 variables, to get the coefficients of the best approximating polynomial of degree d. Then use polyval to evaluate the polynomial at the x data points and, say, 10 equally spaced points between each pair of x data points (401 points total). Plot the result, displaying both the approximating polynomial and the original data. What is the smallest choice of d that gives you a reasonably good approximation of the original data? There is more than one possible answer as it is a matter of individual judgment what "good" means here. If you get wildly bad approximations, you have made a mistake.
   • Once you've chosen d, try computing the least squares solution in several other ways:
     • explicitly solving the normal equations using \ (see Sections 7.4 and 12.1)
     • explicitly solving the normal equations using chol (see Section 5.5)
     • using the economy-sized version of the built-in qr, solving Rc = Q'y (see Section 7.4 but note that the notation is different there, in particular b is used for our y and c is used for our Q'y)
     • same but using your qrrot
     • using the SVD (see Section 7.4)
     Print the coefficients c to full machine precision. Do you see any differences?
3. Numerical Integration (Ascher-Greif, Chapter 15, Sections 1, 2 and 4)
   • Write an adaptive quadrature routine [intval, numevals] = myquad(fun, a, b, fa, fb, tol) as follows
     • the first input argument, fun, is an anonymous function as in Homework 3
     • the next two arguments are assumed to satisfy a < b
     • the next two arguments are assumed to be the values of the function fun at a and b respectively
     • the last argument is a tolerance
     • the routine should compute T1, the value of the simple trapezoidal rule on [a,b], and T2, the sum of the two trapezoidal rules on the left and right halves of the interval: for this you need to evaluate the function at (a + b)/2
     • for reasons that we will discuss in class, the a posteriori error estimate for T2 is (T2 - T1)/3
     • if the absolute value of this quantity is less tol, return the value T2 for intval and 1 for numevals
     • otherwise, the routine should recursively call itself on the two half-intervals, returning the sum of the two intval values as the final intval and the sum of the two numevals values, plus 1, as the final numevals, with the value of tol divided by 2 in each of the 2 recursive calls
     • this program is similar to the one in Section 4.3 and the quadtx code on the Moler web page, but much simpler than either. It will probably more confusing to try to understand the other codes than it is to just write this one from scratch
     • start by testing your code with fun set to the sin function on [a,b]=[0,pi], for which you know the exact answer
     • If the code goes into an infinite recursion, you will need to debug it. Try larger tolerances first.
     • then test your code on Example 15.12 of Ascher-Greif and compare the number of function evaluations required for the various tolerances shown on p.359. Your code will need more evaluations, but that's because it's a simpler code.
     • For what degree polynomials does myquad return the exactly correct answer? Check this experimentally by calling it to integrate the polynomial x ⁿ for n=0,1,2,... Explain why this happens.
     • How many function evaluations are required in each of these cases? Why?
   • Finally, make myquad more accurate with a simple change: when it returns because the tolerance is satisfied, return the extrapolated value (4*T2 - T1)/3 instead of T2 (again to be explained in class). How do the answers to the previous three questions change? How does the actual final accuracy of the solution change?
   • Another possibility would be to return Simpson's rule (p.335) on [a,b] instead of T2. How does this compare with using (4*T2 - T1)/3? Explain.