Homework 7. Assigned: Tue Apr 17. Due: Fri Apr 27 at 5 p.m.
QR Factorization, Least Squares, and Numerical Integeration
 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 "economysize" version
of the QR factorization, just as qr does.
 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 economysized version of the builtin 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?
 Numerical Integration (AscherGreif, 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 halfintervals,
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 AscherGreif 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 ^{n} 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.
Reading: Sections 7.4, 12.1, 15.1 and 15.4