Numerical Methods I
CSCI-GA.2420-001, MATH-GA.2010-001
Fall 2011, Thursdays 5:10-7:00pm, WWH Room 101
Instructor: Olof
B. Widlund
Coordinates
Office: WWH 612
Telephone: 212-998-3110
Office Hours: Tuesdays and Thursdays, 4:00 - 5:00 pm.
e-mail: widlund at cims.nyu.edu
Exams
The course has traditionally ended with individual oral exams lasting
about 30 minutes. Given that the enrollment is over 40, this is
unrealistic
and therefore there will be a regular in-class exam
during the exam period. It will take place on Thursday December 22,
one week after
the last lecture. You will be allowed to bring two regular sized
sheets of paper covered on both sides with notes for your private
use during the 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.
Homework 1 Assigned September 15,
due September 30, at 10 am.
Homework 2 Assigned September 26,
due October 7, at 10 am.
Homework 3 Assigned September 28,
due October 14, at 10 am.
Homework 4 Assigned September 28,
due October 21, at 10 am.
Homework 5 Assigned November 3,
due November 11, at 10 am.
Homework 6 Assigned November 8,
due November 18, at 10 am.
Homework 7 Assigned November 21,
due December 9, at 10 am.
Lectures
September 8: Singular value decomposition; cf. Lectures 4 and 5 in the
Trefethen-Bau book. This lecture was given by Professor Jonathan Goodman.
September 15: Projectors. Gram-Schmidt factorizations in its classical and
modified forms. Gram-Schmidt as a factorization in terms of an orthogonal matrix
and an upper triangular matrix. Householder reflections and how to use them to
obtain a QR factorization. A first few words on linear least squares.
September 22: Two handouts on IEEE Double Precision and polynomial
interpolation and applications. Comments on the computation of the vectors
that define Householder transformations: we need to sum the squares of
certain numbers, compute a square root, and make sure that we do not suffer
cancellation. How to study the stability of the square root function; the
relative error is preserved. Linear least squares problems and three approaches
using the normal equations, the QR or the SVD factorizations. Polynomial
interpolation, Vandermonde matrices, Lagrange's formula, and Newton's triangular
table. How to derive cubic Hermite interpolation from the Newton scheme and
taking a limit.
September 29: More on Newton's method for polynomial interpolation and
cubic Hermite interpolation. Error bound for polynomial interpolation. Runge's
example. Piecewise polynomial interpolation in particular a derivation of
cubic spline interpolation. A couple of simple numerical quadrature algorithms.
(This lecture is based on one of last week's handouts. There was also a new
handout on floating point from Overton's book.)
October 6: More on error bounds for numerical integration, in particular
for Simpson's Rule. Adaptive Simpson. Inner product spaces defined by
integrals with weights. The construction of sets of orthogonal polynomials;
they satisfy three term recurrence formulas and their zeros are all real,
simple, and they all lie in the interval of integration. (Orthogonal
polynomials are developed in many books, e.g., the one by Isaacson and Keller.)
October 13: Guest lecture by Professor Michael Overton on IEEE floating
point systems.
October 20: More on orthogonal polynomials. They can also be characterized
as eigenfunctions of certain self-adjoint differential operators. That the
operators are self-adjoint is established by integration by parts; note that
the eigenvalues of a self-adjoint operator are real, just like those of a
symmetric matrix. Gaussian quadrature: the space of polynomials which is
integrated exactly if the largest possible given the number of parameters;
a handout on these topics will be provided on October 27. Gauss-Lobatto
quadrature. Returning to part III of the text book: Wilkinson's example
of a polynomial of degree 20 the roots of which are very sensitive to
perturbations of the coefficients. This suggests that such problems should
be avoided if there is an alternative route to a solution of a numerical
problem; in particular, we should not compute eigenvalues of matrices
via their characteristic polynomials.
October 27: A further comment on Gauss-Lobatto-Legendre quadrature.
Condition numbers of function evaluation recalled. Stability and backward
stability of numerical algorithms; we succeed in separating ill-conditioning
of problems from issues related to the quality of an algorithm. Examining
the QR factorization using Householder reflections by considering the
numerical example in Trefethen-Bau's book. A proof that this factorization
can be used successfully to solve a linear system of algebraic equations
under the assumption that the three algorithms that are used each is
backward stable. A proof that the solution of solving linear systems with
upper triangular coefficient matrices is backward stable; this requires
the examination of each individual arithmetic operation.
November 3: A few further comments on the proof of backward stability of
solving linear system with upper triangular matrices. Conditioning of linear
least squares problems; there are three relevant parameters, of which one is
the condition number, and four condition numbers. Two of these formulas are
easy to derive, the other two not. How numerical experiments, carefully
selected, can be coupled with information on the conditioning of these problems
to help determine which numerical methods are backward stable. Note that the
problem used in Lecture 19 is carefully designed to demonstrate the failure
of using the normal equations to solve the least squares problem. An
introduction to Gaussian elimination and partial pivoting; watch out for growth
of the computed elements. This can never happen if we use a QR rather than a
LU factorization.
November 10: More about partial pivoting and Gaussian elimination; even
very easy small problems can fail if partial pivoting is not used. The effect
of growth of the elements of U, the upper triangular matrix of the LU
factorization; if the norm of U is not very much larger than the norm of A,
we will have backward stability. Comments on sparse matrix factorization.
How sparsity is lost and a way of describing this by using a graph of the
given matrix and those of the smaller problems obtained after eliminating
some of the variables. By running a graph algorithm we can determine, for a
given elimination order, where zeros will appear in L and U. Once this is
known, a sparse matrix representation can be set up for the matrix, with storage
allocated exactly for the nonzero elements of L and U. Note that this can be
in conflict with our requirement of backward stability when we start working
with the actual matrix elements. However, this is not a problem if the matrix
is positive definite, symmetric. Then we can use Cholesky's algorithm, which
is always backward stable. A quick review of the algebraic eigenvalue problem.
Cases when the matrix has a full set of eigenvalues. Schur factorization. How
to turn a given matrix into an upper Hessenberg matrix with the same
eigenvalues; this is a preliminary step before we use an iterative method
to find the eigenvalues.
November 17: Normal matrices and the use of the Schur factorization
to prove that they can be diagonalized by using an orthogonal set of
eigenvectors. Formulation of the Bauer-Fike theorem, which shows that the
eigenvalues of normal matrices are insensitive to perturbation of the
matrix elements. The power method and the inverse power method. Rayleigh
quotients and Rayleigh quotient iterations. The QR algorithm without and with
shifts. Jacobi's method with some indications why it converges. Eigenvalues
of symmetric tridiagonal matrices and Cauchy's interlace theorem.
November 24: Thanksgiving.
December 1: Cauchy's interlace theorem, Sturm sequences and how to combine
them with a simple dissection algorithm to find eigenvalues of symmetric
tridiagonal matrices. Inverse iteration and how combine knowledge of the
norm of the solution and a perturbation argument. Simultaneous iteration
and the QR algorithm and how to use the analysis of the former to understand
the latter. How inverse iteration fits into this picture. Some first comments
on Krylov spaces and the conjugate gradient algorithm to solve linear systems
of equations with positive definite, symmetric matrices.
December 8: Computing the SVD. Golub-Kahan biorthogonalization and
alternatives. Formulating the solution of linear systems of algebraic
equations, with a symmetric positive definite matrix, as a minimization
problem. Approximation of its solution in a Krylov
space defined by the matrix and the initial residual. Derivation of formulas
for the conjugate gradient (CG) algorithm; a handout on December 15.
A polynomial
approximation problem related to the CG algorithm. A few words about
Richardson's method; convergence can be established even for nonsymmetric
coefficient matrices provided that the symmetric part of the matrix is
positive definite.
December 15: Three handouts, one on the CG method and two on polynomials.
How best approximation problems in the maximum norm arise in error bounds
for polynomial interpolation and the conjugate gradient method; the best
choice is in terms of Chebyshev polynomials. Two theorems on the best
approximation and the oscillation of the error. Error bounds for the conjugate
gradient method in terms of the condition number of the system matrix. How to
derive a preconditioned conjugate gradient method. The conjugate residual
algorithm for systems with symmetric but indefinite matrices;
the same Krylov space as for the positive definite case can be used but we
have to work with a different inner product defined by A^2. A few words on
Arnoldi's method and GMRES which can be used for nonsymmetric problems.
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.