Scientific Computing, Spring 1999 (http://www.math.nyu.edu/faculty/goodman/teaching/sci comp 99/)

Assignment 1.

Given January 12, due February 2.

Objective: To explore computer arithmetic.

1. The fibonacci numbers, f_{k}, are defined by f_{0} = 1, f_{1} = 1, and

f_{k+1} = f_{k} + f_{k−1}  (1)

for any integer k > 1. A small perturbation of them, the “pib numbers” (“p” instead of “f” to indicate a perturbation), p_{k}, are defined by p_{0} = 1, p_{1} = 1, and

p_{k+1} = c · p_{k} + p_{k−1}

for any integer k > 1, where c = 1 + \sqrt{3}/100.

a. Make a plot of log(f_{n}) and log(p_{n}) as a function of n. On the plot, mark 1/ϵmach for single and double precision arithmetic. This can be useful in answering the questions below.

b. For various n values, compute the f_{k} for k = 2, 3, …, n using (1). Then rewrite (1) to express f_{k−1} in terms or f_{k} and f_{k+1}. Use the computed f_{n} and f_{n−1} to recompute f_{k} for k = n−2, n−3, …, 0. Make a plot of the difference between the original f_{0} = 1 and the recomputed f_{0} as a function of n. What n values result in no accuracy for the recomputed f_{0}? How do the results in single and double precision differ?

c. Repeat b. for the pib numbers. Comment on the striking difference in the way precision is lost in these two cases. Which is more typical? Extra credit: predict the order of magnitude of the error in recomputing p_{0} using what you may know about recurrence relations and what you should know about computer arithmetic.

2. The binomial coefficients, a_{n,k}, are defined by

a_{n,k} = ( \begin{matrix}n \\ k \end{matrix} ) = \frac{n!}{k!(n − k)!}  (2)

To compute the a_{n,k}, for a given n, start with a_{n,0} = 1 and then use the recurrence relation a_{n,k+1} = \frac{n−k}{k+1}a_{n,k}.

(a) For a fixed of n, compute the a_{n,k} this way, noting the largest a_{n,k} and the accuracy with which a_{n,n} = 1 is computed. Do this in single and double precision. Why is it that roundoff is not a problem here as it was in problem (1)? b. Use the algorithm of part (a) to compute

E(k) = \frac{1}{2^{n}} ∑_{k=0}^{n} ka_{n,k} = \frac{n}{2} .  (3)

In this equation, we think of k as a random variable, the number of heads in n tosses of a fair coin, with E(k) being the expected value of k. This depends on n. Write a program without any safeguards against overflow or zero divide (this time only!)^{1}. Show (both in single and double precision) that the computed answer has high accuracy as long as the intermediate results are within the range of floating point numbers. As with (a), explain how the computer gets an accurate, small, answer when the intermediate numbers have such a wide range of values. Why is cancellation not a problem? Note the advantage of a wider range of values: we can compute E(k) for much larger n in double precision. Print E(k) as computed by (4) and M_{n} = max_{k} a_{n,k}. For large n, one should be inf and the other NaN. Why?

c. For fairly large n, plot a_{n,k} as a function of k for a range of k chosen to illuminate the interesting “bell shaped” behavior of the a_{n,k} near their maximum.

d. (Extra credit, and lots of it!) Find a way to compute

S(k) = ∑_{k=0}^{n} (−1)^{k} sin(2π sin(k/n))a_{n,k}

with good relative accuracy for large n. This is very hard, so don’t spend too much time on it.

1One of the purposes of the IEEE floating point standard was to allow a program with overflow or zero divide to run and print results