Description
I. We have seen that using computer arithmetic, we can evaluate polynomials exactly (up
to the rounding error). Let the polynomial p(x) of degree n be given by
p(x) = Xn
k=0
akx
k
. (1)
2
Note that p(x) can be rewritten as
p(x) = a0 + x(a1 + x(a2 + x(a3 + · · · + x(an−1 + xan)· · ·))). (2)
Explain why using the expression in (2) is more efficient over the expression in (1)for
evaluating the polynomial. Write a Matlab program for implementing the polynomial
evaluation using (2) and use it to compute π with the help of
pn(x) = Xn
k=0
(−1)k x
2k+1
2k + 1
, |x| ≤ 1
the nth order Taylor series approximation to tan−1 x. Complete the following table
n |π − 6pn(1/
√
3)|/π |π − 4pn(1)|/π
8 – –
16 – –
32 – –
64 – –
128 – –
256 – –
512 – –
1024 – –
2048 – –
4096 – –
Comment on the result.
II. Consider the function
f(x) = (e
x − 1)/x.
Use L’Hospital’s rule to show that
limx→0
f(x) = 1.
(a) Check this result empirically by writing a program to compute f(x) for x =
10−k
, k = 1, 2, . . . , 15. Do your results agree with theoretical expectations? Explain why.
(b) Perform the experiment in part (a) again, this time using the mathematically
equivalent formulation
f(x) = (e
x − 1)/ log(e
x
),
evaluated as indicated, with no simplification. Does this work better? If yes,
explain why?
3
III. Show that fj,k = sin(x0 + (j − k)π/3) satisfies the recurrence relation
fj,k+1 = fj,k − fj+1,k. (3)
We view this as a formula that computes the f values on level k + 1 from the f
values on level k. Let f
a
j,k for k ≥ 0 be the floating point numbers that come from
implementing fj,0 = sin(x0 + jπ/3) and (3) (for k > 0) in double precision floating
point. If
f
a
j,k − fj,k
≤ for all j, show that
f
a
j,k+1 − fj,k+1
≤ 2 for all j. Thus, if the
level k values are very accurate, then the level k + 1 values still are pretty good.
Write a program that computes ek =
f
a
1,k − f1,k
for 1 ≤ k ≤ 60 and x0 = 1. Print
the ek and see whether they grow monotonically. Plot the ek with respect to k and see
that the numbers seem to go bad suddenly at around k = 50.
Repeat this computation in single precision and comment on the result. Note that, if
you are using MATLAB, by default, it does computations in double precision.
IV. The Fibonacci numbers, fk, are defined by f0 = 1, f1 = 1, and
fk+1 = fk + fk−1 (4)
for any integer k > 1. A small perturbation of them, the pib numbers (“p” instead of
“f” to indicate a perturbation), pk, are defined by p0 = 1, p1 = 1, and
pk+1 = c pk + pk−1
for any integer k > 1, where c = 1 + √
3/100.
(a) Plot the fn and pn together on a log scale plot. On the plot, mark 1/mach for single
and double precision arithmetic. This can be useful in answering the questions
below.
(b) Rewrite (4) to express fk−1 in terms of fk and fk+1. Use the computed fn and
fn−1 to recompute fk for k = n − 2, n − 3, . . . , 0. Make a plot of the difference
between the original f0 = 1 and the recomputed ˆf0 as a function of n. What n
values result in no accuracy for the recomputed f0? 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?
V. The binomial coefficients, an,k, are defined by
an,k =
n
k
=
n!
k!(n − k)!.
To compute the an,k, for a given n, start with an,0 = 1 and then use the recurrence
relation
an,k+1 =
n − k
k + 1
an,k.
4
(a) For a range of n values, compute the an,k this way, noting the largest an,k and the
accuracy with which an,n = 1 is computed. Do this in single and double precision.
Why is roundoff not a problem here as it was in problem IV? Find n values for
which ˆan,n ≈ 1 in double precision but not in single precision. How is this possible,
given that roundoff is not a problem?
(b) Use the algorithm of part (a) to compute
E(k) = 1
2
n
Xn
k=0
kan,k =
n
2
. (5)
Write a program without any safeguards against overflow or zero divide. 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 (5) and
Mn = maxk an,k. For large n, one should be inf and the other NaN. Why?
(c) For fairly large n, plot an,k/Mn as a function of k for a range of k chosen to illuminate the interesting “bell shaped” behavior of the an,k near k = n/2. Combine
the curves for n = 10, n = 20, and n = 50 in a single plot. Choose the three k
ranges so that the curves are close to each other. Choose different line styles for
the three curves.