Description
The second problem set is drawn from the text.
• (4.7.2) Use MATLAB to plot the function f(x) = (5−x)exp(x)−5 for x
between 0 and 5. (This function is associated with the Wien radiation
law, which gives a method to estimate the surface temperature of a
star.)
(a) Write a bisection routine or use routine “bisect” available
from the books web page to find a root of f(x) in the interval
[4, 5], accurate to six decimal places (that is, find an interval
of width at most 10−6
that contains the root, so that the
midpoint of this interval is within 5 × 10−7 of the root). At
each step, print out the endpoints of the smallest interval
known to contain a root. Without running the code further,
answer the following: How many steps would be required
to reduce the size of this interval to 10−12? Explain your
answer.
(b) Write a routine to use Newtons method or use routine
“newton” available from the books web page to find a root of
f(x), using initial guess x0 = 5. Print out your approximate
solution xk and the value of f(xk) at each step and run
until |f(xk)| ≤ 10−8
. Without running the code further, but
perhaps using information from your code about the rate at
which |f(xk)| is reduced, can you estimate how many more
steps would be required to make |f(xk)| ≤ 10−16 (assuming
that your machine carried enough decimal places to do this)?
Explain your answer.
(c) Take your routine for doing Newtons method and modify
it to run the secant method. Repeat the run of part (b),
using, say, x0 = 4 and x1 = 5, and again predict (without
running the code further) how many steps would be required
to reduce |f(xk)| below 10−16 (assuming that your machine
carried enough decimal places to do this) using the secant
method. Explain your answer.
1
• (4.7.19) Finding the area of 96-sided polygons that inscribe and circumscribe a circle of diameter 1, Greek mathematician Archimedes
(287-212 BC) determined that 223/71 < π < 22/7. Almost two thousand years later, John Machin (1680-1752) exploited a recent discovery
of his contemporary Brook Taylor, the Taylor series. He also employed
the following trigonometric identity that he discovered:
π = 16 arctan(1/5) − 4 arctan(1/239). (1)
(a) Find the Maclaurin series for arctan(x).
(b) Write down Pn(x), the nth degree Taylor polynomial for
arctan(x) centered at x = 0.
(c) Approximate π by using Pn(x) and (1). More specifically,
use the approximation
π ≈ Tn = 16 Pn(1/5) − 4 Pn(1/239).
Using MATLAB, find Tn and the absolute error |Tn − π| for
n = 1, 3, 5, 7, and 9. How many decimal places of accuracy
did you find for n = 9?
• (6.3.2) Let f(x) = x
1/3
.
(a) Find the absolute and relative condition numbers of f.
(b) Where is f well-conditioned in an absolute sense? In a
relative sense?
(c) Suppose x = 10−17 is replaced by ˆx = 10−16 (a small
absolute change but a large relative change). Using the
absolute condition number of f, how much of a change is
expected in f due to this change in the argument?
• (7.7.3) Let
A =
10−16 1
1 1 !
and b =
2
3
!
.
Here we will see the effect of using a tiny element as a pivot.
2
(a) By hand, solve the linear system Ax = b exactly. Write
your answer in a form where it is clear what the approximate
values of x1 and x2 are.
(b) In MATLAB, enter the matrix A and type cond(A) to
determine the 2-norm condition number of A. Would you
say that this matrix is well-conditioned or ill-conditioned in
the 2-norm?
(c) Write a MATLAB code (or use one from the text) that
does not do partial pivoting to solve the linear system
Ax = b. Compare the answer returned by this code to
the one that you determined by hand and to the one that
you get in MATLAB by typing A\b (which does do partial
pivoting).
• (7.7.10) Show that for all n-vectors v:
(a) kvk∞ ≤ kvk2 ≤
√
nkvk∞.
(b) kvk2 ≤ kvk1.
(c) kvk1 ≤ nkvk∞.
For each of the above inequalities, identify a nonzero vector v (with
n > 1) for which equality holds.
• (7.7.14) Find the polynomial of degree 10 that best fits the function
b(t) = cos(4t) at 50 equally spaced points t between 0 and 1. Set up the
matrix A and right-hand side vector b, and determine the polynomial
coefficients in two different ways:
(a) By using the MATLAB command x = A\b (which uses
a QR factorization).
(b) By solving the normal equations AT Ax = ATb. This
can be done in MATLAB by typing x = (A0 ∗ A)\(A0 ∗ b).
Print the results to 16 digits (using format long e) and comment on
the differences you see. (Note: You can compute the condition number
of A or of AT A using MATLAB function cond.)
3
