M.I.T., Department of Mathematics

18.335: Numerical Methods of Applied Mathematics -- I, Fall 2004

Where and when: 1-390, MW 3-4:30

Introduction: This course will consist of two parts. During the first two thirds of the course we will concentrate on Numerical Linear Algebra. We will study the solutions of linear systems of equations, least square problems, eigenvalue problems and singular value problems. Techniques for dense, sparse and structured problems will be covered. Students should still come to appreciate many state-of-the-art techniques and recognize when to consider applying them. We will also learn basic principles applicable to a variety of numerical problems and learn how to apply them. These principles include (1) matrix factorizations, (2) perturbation theory and condition numbers, (3) effect of roundoff on algorithms, including properties of floating point arithmetic, (4) analyzing the speed of an algorithm, (5) choosing the best algorithm for the mathematical structure of your problem, and (6) engineering numerical software. In addition to discussing established solution techniques, open problems will also be presented. The course assumes familiarity with linear algebra and will involve a reasonable amount of programming in MATLAB. During the second part of the course we will concentrate on numerical methods for solving ordinary differential equations. These methods are usually introduced in undergraduate numerical analysis courses such as 18.330. Such courses, however, are not prerequisites for 18.335. This graduate-level exposition will be self contained and lecture notes will be provided.

Textbooks: Numerical Linear Algebra, Trefethen and Bau, SIAM 1997
(Optional) Applied Numerical Linear Algebra, Demmel, SIAM 1997.
Lecture notes for Numerical ODEs will be provided

Instructor: Plamen Koev, office: 2-376, phone: 253-5013, e-mail: plamen (at) math (dot) mit (dot) edu
Office hours:
Monday 4:30-6:30

Exam Dates: November 3rd. Practice midterm: pdf.

Grading: About eight homework assignments (80%), and one in class midterm on November 3, 2004

Lectures
09/08 Introduction, Examples, Matrix-vector and Matrix-Matrix products
09/13 Orthogonal matrices, norms of vectors and matrices
09/15 The singular value decomposition
09/20 QR decomposition
09/22 Givens rotations and Householder reflectors
09/27 Least Squares Problems
09/29 Floating point arithmetic
10/04 Accuracy, Stability and Conditioning
10/06 Backward Stability
10/13 Conditioning of Least Squares Problems
10/18 Gaussian Elimination
10/20 Cholesky Factorization
10/25 Eigenvalue Problems
10/27 Hessenberg reduction
11/01 QR Iteration I
11/03 Midterm PDF
11/08 QR Iteration II -- Simultaneous iteration; Jacobi algorithm
11/10 Bisection and divide-and-conquer algorithms
11/15 SVD algorithms, LR interation, dqds
11/17 Iterative algorithms, Arnoldi
11/22 Lanczos algorithm
11/24 Conjugate gradients
11/29 ODEs Handout
12/01 Runge Kutta Methods
12/06 Solution to stiff ODEs

Demonstration software used in class: Householder QR, Givens QR, Clown, Secular Equation, Golub-Kahan bidiagonalization, Jacobi method, Arnoldi, Arnoldi Driver, Conjugate gradient preconditioners.
Homework:
  1. Due: September 22. Solutions.
    1. Let A be an orthogonal matrix. Prove that |det(A)|=1. Show that if B is also orthogonal and det(A)=-det(B) then A+B is singular.
    2. Trefethen 2.5, 3.2, 3.3
    3. Prove that ||xy*||F=||xy*||2=||x||2||y||2 for any x,y in Cn
  2. Due: September 29 Solutions.
    1. Count the number of floating point operations required to compute the QR decomposition of an m-by-n matrix using (a) Householder reflectors (b) Givens rotations.
    2. Trefethen 5.4
    3. If A=R+uv*, where R is upper triangular matrix and u and v are (column) vectors, describe an algorithm to compute the QR decomposition of A in O(n2) time.
    4. Given the SVD of A, compute the SVD of (A*A)-1, (A*A)-1A*, A(A*A)-1, A(A*A)-1A* in terms of U, Sigma and V.
  3. Due October 6. Solutions.
    1. Trefethen 10.1
    2. Let B be an n-by-n upper bidiagonal matrix. Describe an algorithm for computing the condition number of B measured in the infinity norm in O(n) time.
  4. Due October 20. Solutions. Solve at least 4 of the problems below.
    1. Trefethen 11.1
    2. Trefethen 13.3
    3. Trefethen 13.4
    4. Prove that (13.7) in Trefethen is valid for complex arithmetic (all four arithmetic operations) with |epsilon| now bounded by a modest multiple of epsilon_machine.
    5. Prove that in IEEE binary floating point arithmetic sqrt(x^2) returns x exactly.
    6. Let a and b be positive IEEE binary floating point numbers such that a < b < 2a. Prove that fl(b-a)=b-a exactly.
  5. Due October 27. Solutions. Solve 4 of the following 5 problems:
    1. Trefethen 20.1
    2. Trefethen 21.6
    3. Trefethen 22.1
    4. Let A be symmetric and positive definite. Show that |aij|2 < aii ajj.
    5. Let A and A-1 be given real n-by-n matrices. Let B=A+xyT be a rank-one perturbation of A. Find an O(n2) algorithm for computing B-1. Hint: B-1 is a rank-one perturbation of A-1.
  6. Due November 3. Solutions.
    1. Let A be skew Hermitian, i.e. A*=-A. Show that (I-A)-1(I+A) is unitary.
    2. Trefethen 25.1
  7. Due November 17. Solutions.
    1. Compute the smallest eigenvalue of the 100-by-100 Hilbert matrix H=1/(i+j-1). (Hint: The Hilbert matrix is also Cauchy. The determinant of a Cauchy matrix C(i,j)=1/(xi+yj) is
      det C = [\prodi< j(xj-xi)(yj-yi)]/ [\prodi,j(xi+yj)].
      Any submatrix of a Cauchy matrix is also Cauchy. You can use Cramer's rule in order to compute accurate formulas for H-1 and then compute its largest eigenvalue.
    2. Trefethen 30.2
  8. Due November 24. The writeups should be brief (preferably < 1 page) and include 1) the answer 2) the method you used and why you think it works 3) one paragraph on how you implemented the method. You need to compute the answer with at least 9 correct decimal digits.
    1. Compute the smallest eigenvalue of the 100-by-100 matrix H=1/(i+j).
      Solution: See the solution to Problem 7.1 above.
    2. The infinite matrix A has entries A(1,1)=1, A(1,2)=1/2, A(2,1)=1/3, A(1,3)=1/4, A(2,2)=1/5, A(3,1)=1/6, etc. Compute ||A||2.
      Solution: See this website, and problem 3 (and its solution) there.
  9. Due December 1.
    1. A is a square matrix of size 19,000. All entries of A are zero except for the primes 2,3,5,7,...,212369 along the main diagonal and the number 1 in all the positions aij with |i-j|=1,2,4,8,...,16384. Compute the (1,1) entry of A-1.
      Solution: See this website, and problem 7 (and its solution) there.
    2. The 200-by-200 (diagonally dominant) matrix A has offdiagonal entries -1/i-1/j and row sums s(i)=sum(A(i,:))=22i. Compute the smallest in magnitude eigenvalue of A. (Hint: The diagonal entries A(i,i) can then be computed as sum of positive (why?) numbers, thus to high relative accuracy. Abandoning the row sums s(i), however, robs you of any chance of computing the smallest eigenvalue accurately). You will encounter no (true) subtractions when running Gaussian elimination, obtaining the LU decomposition of A to high relative accuracy componentwise. Using the LU factors to compute the inverse of A one column at a time will not result in any subtractions either, yielding a positive A-1.)
      Solution: PDF, hw6.m, InverseMM.m. Answer: 1.207993236710136e+002
  10. Due December 8. Compute the smallest eigenvalue of the 200-by-200 Pascal matrix. Details on the Pascal matrix can be obtained here.
    Solution: While there are many approaches to solve this problem, one way would be to compute the largest eigenvalue of the inverse. The inverse (which has checkerboard sign pattern) can (once again) be computed without performing any subtractions if one takes the correct approach. Instead of eliminating the matrix in the typical Gaussian elimination fashion, try to eliminate it by using only ADJACENT rows and columns. This process is called Neville elimination. Once you eliminate the first row and first column, you will see that the Schur complement is also a Pascal matrix of one size less. In matrix form this elimination can be written as

    L * Pn * LT = P'n-1

    where L is lower bidiagonal matrix with ones on the main diagonal and -1 on the first subdiagonal and P'n-1 is an n-by-n matrix with zeros in the first row and column, except for the (1,1) entry (which equals one), and the matrix Pn-1 in the lower right hand corner. You can now observe (no need to prove) that if you have (Pn-1)-1 you can compute (Pn)-1 using the above equality without performing any subtractions.