Numerical Methods I
G22.2420, G63.2010
Fall 2009, Thursdays 5:00-7:00pm, WWH Room 101
Instructor: Olof
B. Widlund
Coordinates
Office: WWH 612
Telephone: 212-998-3110
Office Hours: Thursdays 4:00 -- 5:00 pm and Tuesdays 1:30 -- 2:30 pm.
Or try to
stop by my office any time, or send email or call for an appointment.
Email: widlund@cims.nyu.edu
Notes
The lectures will begin at 5:00 pm instead of 5:10 pm.
There will be no class on Thursday November 26, which is Thanksgiving
Day.
The final lecture will be on Tuesday December 15, which will
run on a Thursday schedule.
Oral Exams
The course ends with individual oral exams which last about 30 minutes.
If you wish to take your exam after Christmas, I will give you an I
and then change the grade after your oral exam.
Requirements
Regular homework assignments, including programming
using Matlab or any reasonable alternative.
It is important that you do the homework yourself, but when you get stuck,
I encourage you to consult with other students
or me, to get help when necessary. However, when you get help,
it's important to acknowledge it in writing. Passing off other people's
work as your own is called plagiarism and is not acceptable.
Please staple everything together and order the problems
in the same order as given.
Homework Assignments
Homework scores will be made available through the NYU Blackboard system.
You can download the assignments in pdf format from this home
page. Some of the problems involve MATLAB programming. Those parts
of the assignments should be submitted electronically to the grader
An-Sheng Jhang (jhang@cims.nyu.edu).
Homework 1 Assigned September 12,
due September 25, at 10 am.
Homework 2 Assigned September 19,
due October 3, at 10 am.
Homework 3 Assigned September 28,
due October 10, at 10 am.
Homework 4 Assigned October 10,
due October 23, at 4:30 pm.
Homework 5 Assigned October 23,
due November 6, at 4:30 pm.
Homework 6 Assigned October 30,
due November 13, at 4:30 pm.
Homework 7 Assigned November 7,
due December 1, at 4:30 pm.
Homework 8 Assigned November 19,
due December 5 4:30 pm.
Lectures
September 10: Part I of Trefethen-Bau's book, in particular
most of lectures 4 and 5. The singular value decomposition. Reduced SVD.
Existence of the SVD.
Determining the rank of a matrix using its singular values. How to find the
singular values and singular vectors by solving certain algebraic eigenvalue
problems.
September 17:
The SVD revisited. How to identify the four principal
subspaces of a matrix A by using the SVD; these subspaces are the range
and null space of A and the range and null space of the conjugate transpose
A-star of A.
Lectures 6 - 8 and 10 of Trefethen-Bau's book:
Projectors and complementary projectors and the two subspaces
associated with them. How to define a projector given two such subspaces.
Orthogonal projectors and how to compute an orthogonal projector onto the
column space of a full rank matrix. Gram-Schmidt's algorithm in its classical
and modified forms. How to use Gram-Schmidt to compute the QR and reduced
QR factorizations of a matrix. How to use a QR factorization of a square
full rank matrix to solve linear systems of algebraic equations. Solving
triangular systems of equations straightforwardly. Householder reflectors
and how to use them to compute the QR factorization of a matrix. The geometry
of Householder transformations and the best choice of two reflectors.
September 24: Further comments on Householder transformations. How to use
the v_k matrices to compute the product of Q and Q^\star times vectors; Q
is the unitary matrix in the QR factorization of a matrix A. Linear least
squares problem. Historical remarks. Three algorithms to compute the
solution of a linear least squares problem. HANDOUT ON IEEE FLOATING POINT.
Positional number systems. IEEE single and double precision numbers. The
Toy system. Subnormal numbers. The relative precision of floating point
numbers after rounding. What can be said about the result of the standard
binary arithmetic operations if the input is a pair of floating point
numbers. Catastrophic cancellation when subtracting two close numbers
which are not exact. The intrinsic sensitivity of function evaluations.
October 1: What can we learn about the IEEE floating point numbers
by looking at two tables in the handout? +0, -0, +infinity, -infinity,
NaN. Rules of rounding, overflow and graceful underflow using the
subnormal numbers. Two examples of instability using correctly rounded
three digit decimal numbers: i) Computing the variance of a set of numbers
in two ways; one results in a most inaccurate result. ii) An unstable
recursion; here the fix is to run the recursion backwards. The conditioning
of matrix-vector multiplication, solving linear system by using the
inverse of the matrix, and the effect of changes in the system matrix
on the solution of a linear system of algebraic equations. A short introduction
to accuracy, stability, and backward stability of algorithms. A few comments
on extended precision floating point numbers under the IEEE standards.
October 8: HANDOUT ON POLYNOMIALS AND RELATED ALGORITHMS. Polynomial
in power and shifted power forms. Using several, up to n, centers in the
representation of polynomials. An algorithm to exchange the centers, one by
one and its use in demonstrating the we can factor out linear factors
of the form (x-root). Newton's form for representing polynomials. A
recursive definition of divided difference and how to compute them in
a triangular array with about n^2/2 entries. Uniqueness and existence of
interpolation polynomials. The Lagrange and Newton forms. Error estimates
for the interpolation error in terms of the product of a high order difference
and a polynomial that vanishes at the points of interpolation. A formula
for difference quotients in terms of a high order derivative. Hermite
cubic interpolation introduced using Newton formula and going to the
limit as points pairwise coincide. Piece-wise linear approximation.
The beginning of a discussion of cubic splines; to be continued.
October 15. HANDOUT ON Hermite interpolation and Gaussian quadrature.
Deriving the equations for cubic spline interpolation; it is a tridiagonal,
diagonally dominant linear system of algebraic equation. No interchanges
of rows necessary and by exploiting the sparsity (the zero-nonzero structure)
of the matrix, we obtain a linear time algorithm. Comments on interpolation
in several variables; if we have a mesh built from two orthogonal sets of
parallel lines, the problem is relatively easy. On a decomposition into
triangles, it is much harder to obtain interpolants that are continuously
differentiable. Basics on numerical integration. The rectangular, midpoint,
trapezoidal, and Simpson methods. Error bounds based on a mean-value
theorem from calculus and the error bound for polynomial interpolation.
How to use the error bound for Simpson's rule to motivate the design of
the adaptive Simpson method. A few words on Gaussian quadrature.
October 22. Orthogonal polynomials with respect to an inner
product and determined by the interval of integration and a weight function.
A proof showing that the roots of these polynomials always are simple and
in the interval of integration. Legendre, Laguerre, Hermite, Chebyshev, and
Jacobi polynomials. Three-term recursions. The orthogonal polynomials as
eigenfunctions of second order ordinary differential equations. Hermite
interpolation using n+1 points. Gaussian quadrature formulas based on
Hermite interpolation; select the nodes of integration as the zeros of
an orthogonal polynomial to make the terms dependent on the first derivatives
vanish. The backward stability of solving a linear system of algebraic
equations by using the QR factorization of Householder; if the three basic
steps of the algorithm are backward stable, then the entire algorithm is
also backward stable.
October 29. Review of backward stable algorithms. Comments on the
proof that QR-factorization of matrices using Householder's method is
backward stable. The stability of the back substitution algorithm for
linear systems with upper triangular matrices; we obtain componentwise
bounds. The conditioning of linear least squares problem; there
are four condition numbers. The under-determined case and how to find the
minimal norm solution. The experiments of Lecture 19 in the Trefethen-Bau
book; why are they well chosen? The outcome using three methods based on
Householder's method, two based on modified Gram-Schmidt, normal equations
and Cholesky factorization, and SVD factorization.
November 5. Solving linear systems of algebraic equations. Gaussian
elimination, the relation between the multipliers used in the row
operations and the elements in the unit lower triangular matrix of the
LU factorization of the matrix. The role of pivoting and the partial
pivoting rule. The possible growth of the elements of the upper triangular
factor and how it can negatively affect the stability of Gaussian
elimination. Solving banded systems of equations and a more favorable
operation count. The limited loss of sparsity for banded systems even if
partial pivoting is used. Cholesky's algorithm for Hermitian, positive
definite coefficient matrices; it is very stable. Symmetric pivoting
to attempt to preserve sparsity. How fill-in happens and how to describe
fill-in by using an undirected graph and elimination graphs. Using a
graph algorithm to determine the fill-in before factoring a matrix using
Cholesky's algorithm; this simplifies data structure issues.
November 12. Comments on iterative refinement of the solution of
linear systems of algebraic equations. Review of the algebraic eigenvalue
problem. Diagonalization of matrices and other rank-revealing
transformations. A proof of the existence of Schur's normal form and why
it is to be preferred over the Jordan normal form. Transforming matrices
into upper Hessenberg form and tridiagonal form in case the matrix is
symmetric. The Bauer-Fike theorem and what it means for symmetric
matrices. The power method and the inverse power method. The QR algorithm
and the shifted QR algorithm.
November 19. Jacobi's method to compute the eigenvalues of symmetric
matrices using two-by-two rotations imbedded is larger matrices. The classical
and the serial Jacobi methods. Symmetric tridiagonal matrices and when
tridiagonal matrices can be made symmetric by a diagonal similarity
transformation. Cauchy's interlace theorem. Sturm sequences and how to use
them to locate the eigenvalues of tridiagonal matrices by using, e.g., a
bisection algorithm. Perturbation analysis of eigenvalues and how to use the
norm of the residual of an eigenvalue problem to estimate the error of an
approximate eigenvalue. Connections between the simultaneous iteration
method and the unshifted QR method; this gives us tools for the analysis
of the rate of convergence of the QR algorithm in terms of that of the
SI algorithm. The latter can be analyzed using eigen expansions.
December 3. Computing the SVD by reducing matrices to bidiagonal form
and then applying algorithms similar to those for the algebraic eigenvalue
problem. Different alternatives for phase 1 of the computation. Iterative
solvers for linear systems and eigenvalue problems when the matrix is available
only in terms of matrix-vector multiplication. The conjugate gradient method.
Krylov spaces and Lanczos vectors reducing the restriction of the matrix
to the Krylov spaces to tridiagonal or upper Hessenberg form. Deriving the
conjugate gradient algorithm from a variational problem. HANDOUT on best
polynomial approximation in the maximum norm.
December 10. Further discussion of the conjugate gradient method:
orthogonal residuals, which can be expressed in terms of the most recent
Lanczos vector and conjugate, i.e., A-orthogonal search directions. Upper
bounds for the error of the conjugate gradient method in terms of a polynomial
approximation problem. Richardson's method and bounds for its convergence in
the positive definite, symmetric case and the case when the operator has a
positive definite symmetric part. GMRES and Arnoldi's methods and two
polynomial approximation problems which characterize the iterates.
approximation problems.
December 15. Best polynomial approximation in the maximum norm; theorems
by De La Vallee-Poussin and Chebyshev. Best approximation of x^n by polynomials
of degree n-1 and of 1 by xp_{n-1}; Chebyshev polynomials provide the answers.
Derivation of the preconditioned conjugate gradient method and a short
introduction to biconjugate gradient methods.
- Introductions to MATLAB
Cambridge University Engineering Department -"Getting
Started" and Example PDFs.
The official Matlab page - detailed with a lot of information.
Or from the same outfit
Wikipedia thread for Matlab - easy read, simple programs that are easy
to comprehend
If you find good alternatives, please e-mail me the URL.
Required Text
Numerical Linear Algebra by Lloyd N. Trefethen and David Bau, III. SIAM.
Recommended Texts
Numerical Computing with IEEE Floating Point Arithmetic by Michael L. Overton,
SIAM.
Applied Numerical Linear Algebra by James W. Demmel. SIAM
Demmel's book also contains a lot of material relevant to Numerical Methods II.