Homework 7. Assigned: Tue Apr 17. Due: Fri Apr 27 at 5 p.m.
QR Factorization, Least Squares, and Numerical Integeration
Reading: Sections 7.4, 12.1, 15.1 and 15.4
- 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.
- 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
Print the coefficients c to full machine precision. Do you see any differences?
- 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)
- 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 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.