MATH-UA.0252 (Aleks Donev, NYU Courant)

Spring 2021: Numerical Analysis

Assignment 2 (due Monday March 8th 2pm)

1. [Compound interest, 8pts] For a yearly interest rate 0 < r < 1 compounded over n intervals, an amount of money C grows to be

after one year. Let C = 1 and r = 0.025. If n is large, there may be loss of digits when evaluating this using finite-precision arithmetic.

(a) [2pts] If n is extremely large, say n = 10¹⁶, in IEEE double precision arithmetic (try it in Matlab),

when in fact,

What happened?

(b) [2pts] To compute f(C, r, n) without roundoff problems in Matlab, compute first ln f using the (magic!) built-in Matlab function log1p which computes ln (1 + x) without loosing digits even for very small x, and then compute f from its logarithm. Write down the formulas used. Try this for n = 10¹⁶. From now on take n = 10⁸.

(c) [2pts] Using the result from part (b), how many digits of accuracy do you get for f with direct evaluation of (1)?

(d) [1pts] For large n, we can just use the approximation f(C, r, n) = Ce^r . How many digits of accuracy do you get with this approximation?

(e) [1pts] An improved approach for large n is to compute a few terms in the Taylor series expansion (not a trivial calculation per se),

and then use this approximation for small x. How many digits of accuracy do you get using this approach?

Don’t just report answers, explain how you computed this.

2. [Backward substitution implementation, 5pts] [3pts] Write a code for backward sub-stitution to solve systems of the form Ux = b, i.e., write a function x = backward(A,b), which expects as inputs an upper triangular matrix U ∈ R^{n x n}, and a right hand side vector b ∈ Rⁿ , which returns the solution vector x ∈ Rⁿ . The function should find the size n from the vector b and also check if the matrix and the vector sizes are compati-ble before it starts to solve the system. Apply your program for the computation of for x ∈ R⁴ , with

[2pts] How do you know that your code is working correctly?

3. [LU factorization of tridiagonal matrix, 6pt] Given is a tridiagonal matrix, i.e., a matrix with nonzero entries only in the diagonal, and the first upper and lower subdiagonals:

Assuming that A has an LU decomposition A = LU with

derive iterative expressions for di , ei and fi (i.e., how to compute the value for i+ 1 from the values at i, and how to start for i = 1).

Hint: You could check your answer by implementing the formulas in code and checking that LU = A in Matlab for some specific example.

4. [Inverse matrix computation, 8pts] Let us use the LU-decomposition to compute the inverse of a matrix¹.

(a) [2pts] Describe an algorithm that uses the LU-decomposition of an n × n matrix A for computing A^-1 by solving n systems of equations (one for each unit vector).

      ¹This also illustrates that computing a matrix inverse is significantly more expensive than solving a linear system. That is why to solve a linear system, you should never use the inverse matrix!

(b) [2pts] Calculate the floating point operation count of this algorithm. It is OK to use estimates from class/worksheets but write them down so the grader knows what you are doing.

(c) [4pts] Improve the algorithm by taking advantage of the structure (i.e., the many zero entries) of the right-hand side. What is the new algorithm’s floating point operation count?

     [Hint: Consider splitting the solution vector for the k-th equation from part (a) into two pieces, and solve for each piece separately, on paper or using forward/back substitution.]

5. [Stability of the Gaussian elimination, 8pts]

Consider the linear system

Ax = b,

where A is an n×n matrix that has ones on the diagonal, minus ones below the diagonal, and ones in the last column, with all other entries zero. For example, when n = 5, we have

(a) [3pts] Prove that A is invertible for any n, by induction. [Hint: Perform a column operation on A to eliminate the reduce it to a smaller matrix of size n - 1 and ask whether that smaller matrix is invertible under the induction hypothesis.]

(b) [3pts] Now consider the matrix A for some unspecified (arbitrary) n. Perform Gaussian elimination on A to obtain the upper triangular matrix U appearing in the LU factorization A = LU. What is as a function of n?

(c) [2pts] For large n, e.g., n = 2000, what problems can you envision if you try to solve (4) using Gaussian elimination on a computer? Explain.

     [Note: This is one rare example matrix for even Matlab will fail to solve a linear system correctly even though the matrix is well-conditioned, see discussion in Section 7.5 of Practice textbook.]

6. [Matrix square root, 6pts] Newton’s method for finding roots can be extended to matrix-valued functions as well. Here you will devise a Newton method (i.e., generalize the Babylonian method) to compute the square root of a matrix. If it exists, the square root of a real symmetric n×n matrix A is another real square symmetric matrix X such that

X^TX = A

Just like the square root of even a positive number is not unique, the matrix square is not unique (one can roughly think of having to choose n signs, as we will revisit in a future homework once we cover eigenvalue decompositions).

(a) [2pts] In class/worksheet we used derivatives to obtain Newton’s method. Instead of computing derivatives of matrix-valued functions, however, it is useful to think of computing derivatives from a linearization of the function around a given value (this allows to generalize the notion of a derivative and makes it easier to compute in some cases). Set X = X + δX in (5) and keep only the terms that are linear in the ’perturbation” δX. Use this to write down an equation for δX.

(b) [2pts] The equation you obtained in part (a) can be solved explicitly for any n – can you explain why? It is OK if you assume a unique solution exists. Take n = 2 and write down the solution explicitly. [Hint: It is always a good idea to check by plugging in specific numbers.]

(c) [2pts] It would be nice to write down an explicit formula for the solution of the equation you got from part (a) for any n. Do this by assuming that the matrices X and δX commute, i.e., that

X (δX) = (δX) X.

[Hint: Recall that X is symmetric.].

Note: One can prove (6) holds at all iterations if the initial guess commutes with A; if interested, look at the paper “Newton’s Method for the Matrix Square Root” by Nicholas Higham, freely available on the web.