Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit

CS 314 Numerical Methods

Fall 2022

Homework 5

1.  10 points

In this problem, we will do a simple rounding error analysis for oating point multiplication, obtaining expressions for both forward error and backward error. Let x and y be real numbers so that

f1(x) = x(1 + δ ),    f1(y) = y(1 + δy ),

f1(x * y) = f1(x) * f1(y)(1 + δm ),                                         (1)

where δ is the rounding error in representing x, δy the rounding error in representing y, and δm is the rounding error made by oating point multiplication. Each δ is bounded by u, the unit roundoff equal to ∈/2 when round to nearest is used, and ∈ is what the textbook calls the machine epsilon.

(a) First we obtain an upper bound on the relativeforward error. Simplify Equation (1) to obtain the bound

If1(x * y) - x * yI Ix * yI

< 3u + O(u2 ).

(b) For the backward error analysis, we need the Taylor series for ^1 + x for IxI < 1.

Show that

^1 + x = 1 + x/2 - x2 /8 + O(x3 ).

Now simplify Equation (1), distributing the (1 + δm ) term equally to both x and y using the identity (1 + δm ) = (1 + δm )1/2(1 + δm )1/2 . Then use the Taylor series given above to show that

f1(x * y) = x * y ,    where     <  + O(u2 ),    and

Iy - yI      3u

IyI          2

Thus the oating point product is the exact product of slightly changed input data.

2. 20 points

Now we consider the problem of computing the roots of a quadratic equation ax2 +bx+c = 0 in double precision oating point arithmetic. The formula that you have learned in algebra

class

-b ±^b2  - 4ac

2a          ,

could suffer numerical cancellation for one of these roots.  (We will call this the school formula.)

(a) Let a = 1, b = -230 , and c = 4. Compute the two roots x1  and x2  using the school formula in double precision arithmetic. Compute the correctly rounded values of the roots using the roots function in Matlab.  Since b < 0 we have -b > 0, and the solution using the negative square root suffers cancellation error of two nearly equal quantities. Thus one root is incorrectly computed. Report the roots you have computed using the school formula and the roots function.

b +^b2 - 4ac

we can compute x2 by an alternative method. By multiplying x1 and x2 from the school formula, show analytically that x1 * x2  = c/a. Since we have computed x1 already, we can obtain x2 from the formula x2  = c/(a * x1 ). Verify that the correctly rounded value of x2 is now obtained.

The moral here is that we avoid the catastrophic cancellation that happens in one of the roots when floating point arithmetic is employed, by rewriting the formula to avoid that cancellation in oating point arithmetic.

(c) Now suppose that a = 1, b = 229 , and c = 4. How would you correctly compute the two roots using the school formula to obtain one root? (You should not use the roots function for this problem.)

(d) Finally consider a = 2-540 , b = 3 * 2-540 , and c = 2-540 . What happens when you compute b2  in double precision arithmetic?  The numerical difficulty that you have run into is entirely avoidable, since if ax2  + bx + c = 0, then we should be able to factor any constant out of all three coefficients a, b, and c. (In this problem we could factor 2-540  out without changing the problem.) More generally, we could write m = max{IaI, IbI, IcI}, and then a\  = a/m, b\  = b/m and c\  = c/m. The new coefficients lie between -1 and 1.  Then compute the roots of the equation a\ x2  + b\ x + c\  = 0. Report the roots you have computed.

(The issue we see here is called unnecessary numerical underflow.  There is a corre- sponding unnecessary numerical overflow problem as well. This problem also shows that a lot of care is needed to correctly compute the roots of a quadratic equation using floating point arithmetic.)

3. 20 points

For the following functions of x and the intervals given, calculate the absolute and relative condition numbers. Indicate the values of the argument x where the condition number tends to infinity, and thus the problem becomes ill-conditioned.

(a) x-n , where x is any real number and n is an integer > 1.

(b) a , where a is a positive real number. Here x is the independent variable.

(c) tan x, where x e [0 π/2].

(d) ln x, where x > 0.

(e) arccos(x), for -1 < x < 1.

4.  15 points

(a) In lecture we learned how to compute the absolute and relative condition numbers of an inverse function in terms of the condition numbers of the function.  The root finding problem of nding a root r of a function f (x) corresponds to an inverse function problem: find the value of the function f-1 (0) = r. We can use this idea to nd the condition numbers of a root nding problem. Hence find the absolute condition number of nding a root of the function f (x) = x + 5 sin x. This function has a root equal to 0, and two other roots near 4 and 5. Find the absolute condition number of the problem of computing the root at 0. Is this a well-conditioned problem or ill-conditioned problem for this root in the absolute error sense?

(b) Consider the problem of nding the repeated root of the function x2  - 2x + 1 = 0. Compute the absolute condition number of computing the root at x = 1.  Is this ill- conditioned or well-conditioned in the absolute error sense?

(c) Now we can consider the condition number for a xed point computation x = g(x), since this corresponds to nding the root of the function x - g(x) = 0. Compute the absolute condition numbers for the problem of computing the two xed points of the function g(x) = (x2 + 4)/5.

5. 20 points

In this problem we investigate how to rewrite expressions involving subtraction of nearly equal quantities to avoid cancellation error as far as possible. In each case, choose a value of x that shows the difference in accuracy between the original and rewritten expressions, and compute the values obtained by them using double precision arithmetic.

(a) Rewrite sinh(x) = , near x = 0. (Use the Taylor series expansions

of exp(x) and exp(-x) about the point 0, and keep the rst two nonzero terms. Report

the three values you compute for x = 1e -5: (1) ; (2) the expression

you obtained from Taylor series; and (3) sinh(x) computed by Matlab.

(b) Rewrite sin x - cos x, near x = π/4. The problem here is that cos x = sin(x) = 1/^2 at x = π/4, causing cancellation errors for nearby x. Obtain a more accurate formula by writing x = π/4 + y, and then use the cosine and sine addition formulae, to obtain a simplified expression.

sin(a + b) = sin a cos b + cos a sin b,    cos(a + b) = cos a cos b - sin a sin b.

Report two values you compute for sin(x) - cos(x) for x = π/4 + 1e - 15: (1) the value obtained from Matlab for sin(x) - cos(x), and (2) the value computed from the expression you obtained.

6. 20 points

Write a script to compute the time it takes to do a matrix-vector product computed three different ways. For this problem, generate a 1000 x 1000 random matrix A and a 1000 x 1 random vector x. In Matlab, you can use the rand command to do this. We will do 100 repetitions of the matrix vector product y = A * x and nd the time needed to do the computation.  The commands tic and toc will let you compute the time elapsed from the rst command to the second command.  (In order to get accurate timings, you need a script here rather than a collection of commands from the command line, since tic starts a stop-watch, and toc stops it.)

We will compute the matrix vector product using three different algorithms. First, we will use the built-in matrix multiplication operation by typing y  =  A *x  . Second, we compute the product by writing a for loop over the row indices i.  For each i, compute the inner product of the i-th row of A with the vector x. Recall we call this the sdot operation since it uses inner products. Third, we compute the matrix vector product as a linear combination of the columns of A multiplied by coefficients in x; i.e., initialize y to the vector of zeros, and then compute y = y + A(:, j) * x(j), for each column A(:, j). Recall this is the saxpy operation.

(a) First run these experiments (do not forget to repeat the multiplication 100 times). Re- port the run times for the three options.

(b) For the second and third option, count the number of multiplications and additions performed when A is n x n and x is n x 1. Your answer should be given in the form cne + O(ne-1 ), where we need the values of c and e. Show your work.

(c) Rank the three algorithms in the rst part in increasing order of run times. Explain the cause for the differences in the times.