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.

  • 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 (

    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.

    1. 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.