-
The function \(f(x)=\sqrt{x^2+1}-x\) is to be evaluated for \(x=10^{8}\) in double precision. Find its relative condition number at that point, show that the formula as written returns \(0\) , and give a formula that returns the value to full precision.
Solution
Conditioning. \(f'(x)=x/\sqrt{x^2+1}-1=-f(x)/\sqrt{x^2+1}\) , so by Proposition 1.8 \(\kappa_f(x)=\lvert xf'(x)/f(x)\rvert=x/\sqrt{x^2+1}\) , which is less than \(1\) for every \(x\) and essentially \(1\) at \(x=10^{8}\) . The problem is well conditioned.
The formula as written. \(x^2=10^{16}\) , and \(10^{16}+1\) is not a double, since floats near \(10^{16}\) are spaced \(2\) apart; it rounds to \(10^{16}\) . The square root is then exactly \(10^{8}\) , and the subtraction gives \(0\) . The true value is about \(5\times10^{-9}\) : the relative error is \(100\) percent, from an unstable algorithm for a well-conditioned problem, by the cancellation mechanism of Proposition 1.6.
Repair. Multiply and divide by the conjugate: \(f(x)=1/\bigl(\sqrt{x^2+1}+x\bigr)\) . Now the only addition is of positive numbers. In double precision this gives \(1/(2\times10^{8})=5\times10^{-9}\) , correct to full precision; the lost \(1\) under the root matters only in the seventeenth digit of the sum.
-
A cubic polynomial has the values \(-3\) , \(-1\) , \(5\) , \(21\) at \(x=0,1,2,3\) . Find the divided-difference form of the interpolant, use inverse linear interpolation between the two values of opposite sign to estimate its root, and then locate the root to three decimals by any method you can justify.
Answer
First differences \(2,6,16\) ; second \(2,5\) ; third \(1\) . So \(p(x)=-3+2x+2x(x-1)+x(x-1)(x-2)=x^3-x^2+2x-3\) . The sign change is between \(x=1\) and \(x=2\) ; inverse linear interpolation gives \(1+\tfrac{0-(-1)}{5-(-1)}=1.167\) . Since \(p'(x)=3x^2-2x+2>0\) the root is unique, and Newton’s method from \(1.167\) gives \(1.2840\) , \(1.2757\) , \(1.2757\) : the root is \(1.276\) , with the bracket \([1,2]\) as safeguard.
-
The integral \(\int_0^1\sqrt{x}\,e^{x}\,dx\) is needed to six digits. Explain why the composite Simpson rule converges slowly on it, transform it into an integral with a smooth integrand, and state which rule you would then use and how you would confirm the accuracy.
Solution
The integrand \(\sqrt x\,e^{x}\) has an unbounded derivative at \(0\) , so Theorem 6.4 does not apply; as in Example 6.6 the composite Simpson rule converges only like \(h^{3/2}\) . Indeed its errors with \(8\) , \(16\) , \(32\) panels are \(3.5\times10^{-3}\) , \(1.3\times10^{-3}\) , \(4.5\times10^{-4}\) , a factor of \(2.8\) per doubling.
Substitute \(x=t^2\) , \(dx=2t\,dt\) :
\[ \int_0^1\sqrt x\,e^{x}\,dx=\int_0^12t^2e^{t^2}\,dt . \]The new integrand is smooth on \([0,1]\) , and the composite Simpson rule now converges with order four: with \(8\) , \(16\) , \(32\) panels its errors are \(4.8\times10^{-4}\) , \(3.1\times10^{-5}\) , \(2.0\times10^{-6}\) , factors of \(16\) . A Gauss–Legendre rule with a handful of points would do better still.
To confirm six digits without knowing the answer, compute \(S_{32}\) and \(S_{64}\) of the transformed integral, check that the differences \(S_{16}-S_{32}\) and \(S_{32}-S_{64}\) are in the ratio \(16\) , and use \((S_{64}-S_{32})/15\) as the error estimate, as in Section 6.5. The value is \(1.255630\) .
-
A linear system with \(\kappa_2(A)\approx10^{10}\) is solved in double precision by a backward stable method, and the computed solution has relative residual \(10^{-15}\) . State what can and cannot be concluded about the relative error, and what would change if the entries of \(A\) were known only to six digits.
Answer
The small residual certifies only that the solver was backward stable. By Corollary 7.7 the relative error may be as large as \(\kappa\times10^{-15}=10^{-5}\) , so about five digits can be trusted, and no more can be claimed. If the entries of \(A\) are known to six digits, Theorem 7.6(b) allows a relative error of \(\kappa\times10^{-6}=10^{4}\) : the data do not determine the solution at all, whatever the solver.
-
Newton’s method applied to a function produces the successive steps \(\lvert x_{k+1}-x_k\rvert=10^{-1},\,5\times10^{-2},\,2.5\times10^{-2},\,1.25\times10^{-2}\) . Diagnose what is happening, say what accuracy can ultimately be expected in double precision, and give a modification of the iteration.
Answer
The steps halve, so the convergence is linear with ratio \(\tfrac12\) : by Proposition 3.6 the root is a double root. Near it the computed function values are dominated by rounding within about \(\sqrt u\approx10^{-8}\) of the root, as in Example 3.4, so about eight digits are attainable. Use \(x_{k+1}=x_k-2f(x_k)/f'(x_k)\) , or apply Newton’s method to \(f/f'\) .
-
Ten thousand data points \((t_i,y_i)\) are to be fitted by a polynomial of degree \(8\) in \(t\) on the interval \([1000,1001]\) . Explain what goes wrong if the monomial basis and the normal equations are used, and give two changes that make the computation reliable.
Answer
On \([1000,1001]\) the monomials \(1,t,\dots,t^8\) are nearly proportional to one another, so the columns of \(A\) are almost dependent and \(\kappa_2(A)\) is far beyond \(1/u\) ; the normal equations square it. Change the variable to \(s=2(t-1000.5)\in[-1,1]\) and use Chebyshev polynomials \(T_k(s)\) as the basis, which makes \(\kappa_2(A)\) modest; and solve by QR, not by the normal equations.
-
The initial-value problem \(y'=-200\,(y-\sin t)\) , \(y(0)=1\) , is integrated to \(t=10\) by the classical fourth-order Runge–Kutta method. Find the largest stable step, the resulting number of steps, and an alternative method that needs far fewer, with the reason.
Answer
\(\lambda=-200\) , so stability requires \(h<2.785/200=0.0139\) , at least \(719\) steps, although after \(t\approx0.05\) the solution is as smooth as \(\sin t\) . The problem is stiff. An A-stable implicit method such as the trapezoidal method or backward Euler can use \(h\approx0.1\) , about \(100\) steps, with a linear equation to solve at each; accuracy, not stability, then sets the step.
-
For the fixed-point iteration \(x_{k+1}=\tfrac12(x_k+\cos x_k)\) , show that it converges for every real starting value, find the asymptotic rate at the fixed point \(0.7391\) , and compare with the rate of \(x_{k+1}=\cos x_k\) .
Solution
Let \(g(x)=\tfrac12(x+\cos x)\) . Its fixed points satisfy \(x=\cos x\) , and since \(x-\cos x\) has derivative \(1+\sin x\ge0\) , vanishing only at isolated points, it is strictly increasing and has exactly one zero, \(r=0.7391\) .
Theorem 3.2 does not apply on the whole line, because \(g'(x)=\tfrac12(1-\sin x)\) equals \(1\) wherever \(\sin x=-1\) . Argue by monotonicity instead. Since \(g'\ge0\) , \(g\) is nondecreasing. Suppose \(x_0>r\) . Then \(g(x_0)\ge g(r)=r\) , and \(g(x_0)-x_0=\tfrac12(\cos x_0-x_0)<0\) because \(x-\cos x\) is positive to the right of its only zero. So \(r\le x_1<x_0\) , and by induction the iterates decrease and are bounded below by \(r\) . They converge, and the limit is a fixed point of the continuous function \(g\) , hence equal to \(r\) . The case \(x_0<r\) is symmetric, with increasing iterates.
By Theorem 3.3 the asymptotic rate is \(g'(r)=\tfrac12(1-\sin r)=\tfrac12(1-0.6736)=0.163\) . The iteration \(x_{k+1}=\cos x_k\) of Example 3.1 has rate \(\sin r=0.674\) . Averaging the iterate with its image has cut the rate fourfold, from one digit every six steps to one digit every \(1.3\) steps, and has replaced an alternating approach by a monotone one.
-
A function is known only through a routine that returns its values with absolute error up to \(10^{-6}\) . Its derivative at a point is needed. Recommend a difference formula and a step size, estimate the attainable accuracy, and explain why halving the step further would be harmful.
Answer
Use the central difference. With \(\epsilon_f=10^{-6}\) and \(f'''\) of order \(1\) , the best step is \((3\epsilon_f/M_3)^{1/3}\approx0.014\) and the attainable accuracy is about \(1.5\,\epsilon_f/h\approx10^{-4}\) : four digits. Halving the step doubles the rounding term \(\epsilon_f/h\) while the truncation term is already smaller than it, so the total error grows.
-
The \(n\times n\) system arising from a natural cubic spline is solved for \(n=10^{6}\) . State the cost, explain why no pivoting is needed, and say what would be different if the spline were replaced by a single interpolating polynomial through all the points.
Answer
The system is tridiagonal, so the Thomas algorithm solves it in about \(8n=8\times10^{6}\) operations. It is strictly diagonally dominant (Theorem 5.4), so elimination needs no pivoting and is stable. A single interpolating polynomial through \(10^{6}\) equally spaced points would have a Lebesgue constant beyond any floating-point range and would oscillate uncontrollably; it is not a usable object.
-
The trapezoid rule with \(16\) panels gives \(2.000013\) for an integral and with \(32\) panels gives \(2.000003\) . Estimate the error of the second value and give an improved value, stating the assumption involved; then name a class of integrands for which the assumption fails and the trapezoid rule is nevertheless extremely accurate.
Solution
Assume the trapezoid error has the form \(Ch^2\) to leading order, as the Euler–Maclaurin formula says for a smooth non-periodic integrand. Then halving \(h\) divides the error by \(4\) , and by the reasoning of Exercise 6.14
\[ I-T_{32}\approx\frac{T_{32}-T_{16}}{3}=\frac{2.000003-2.000013}{3}=-3.3\times10^{-6}. \]So the second value is too large by about \(3.3\times10^{-6}\) , and the improved value is \(T_{32}+(T_{32}-T_{16})/3=2.000003-0.0000033=1.9999997\) , which is Simpson’s value \(S_{32}\) . With only two computed values the assumed order cannot be tested; a third, \(T_{64}\) , would allow the ratio of differences to be compared with \(4\) .
The assumption fails for smooth periodic integrands integrated over a full period, where every coefficient in the Euler–Maclaurin expansion vanishes. There the error decreases faster than any power of \(h\) , the difference \(T_{32}-T_{16}\) is no guide to the error of \(T_{32}\) , which is far smaller, and extrapolation would make the result worse, as in Example 6.4.
-
Power iteration on a symmetric matrix gives Rayleigh quotients \(5.2000\) , \(5.9010\) , \(5.9870\) , \(5.9983\) at successive steps. Estimate the ratio \(\lvert\lambda_2/\lambda_1\rvert\) , extrapolate the dominant eigenvalue, and describe a way to accelerate the computation once this estimate is in hand.
Answer
The successive differences are \(0.7010\) , \(0.0860\) , \(0.0113\) , with ratios \(0.123\) and \(0.131\) . For a symmetric matrix the Rayleigh quotient error contracts by \((\lambda_2/\lambda_1)^2\) , so \(\lvert\lambda_2/\lambda_1\rvert\approx\sqrt{0.131}=0.36\) . Summing the remaining geometric series, \(\lambda_1\approx5.9983+0.0113\times0.131/(1-0.131)=6.0000\) . To accelerate, switch to inverse iteration with shift \(\mu=6.0\) , whose convergence ratio is tiny because the shift is so close to \(\lambda_1\) , or to Rayleigh quotient iteration.
-
The boundary-value problem \(-y''+100\,y=100\) , \(y(0)=y(1)=0\) is to be solved to four digits. Compare finite differences with shooting for this problem, with reasons, and say what grid evidence you would present.
Answer
The equation \(y''=100y-100\) has solutions growing like \(e^{10x}\) , and the shooting sensitivity is \(\sinh(10)/10\approx1100\) : shooting would work but loses three digits. Finite differences are the safer choice: with \(q=100\ge0\) the matrix is strictly diagonally dominant and the maximum principle holds. The solution has layers of width about \(0.1\) at both ends, so \(h\) must be well below \(0.1\) . Present values at a fixed point on three grids, for example \(n=50,100,200\) , show that the differences fall by a factor of four, and quote the Richardson estimate.
-
Jacobi’s iteration and the conjugate gradient method are applied to a symmetric positive definite system with \(\kappa_2=10^{4}\) arising from a strictly diagonally dominant matrix with \(\rho(G_J)=0.9998\) . Estimate the number of iterations each needs to reduce the error by \(10^{-8}\) , and say how the stopping test should be chosen in view of the condition number.
Answer
Jacobi: \(\ln10^{-8}/\ln0.9998\approx92\,000\) sweeps. Conjugate gradients: \(2(99/101)^{k}\le10^{-8}\) gives \(k\ge956\) . Because the relative error can exceed the relative residual by the factor \(\kappa=10^{4}\) , a residual tolerance of \(10^{-12}\) is needed to guarantee an error of \(10^{-8}\) .
-
A colleague reports that a computed solution of an initial-value problem “converged, because halving the step changed the result by only \(10^{-9}\) ”. Describe one situation in which this evidence is sound and one in which it is misleading, and say what additional computation would distinguish them.
Hint
Sound: the method has known order \(p\) , the step is in the asymptotic regime, and the change is then about \((2^{p}-1)\) times the remaining error. Misleading: both steps lie above a stability threshold or fail to resolve a feature, so both results share the same defect; or the change is small because the method has stalled at a rounding floor. A third computation, giving a ratio of successive changes to compare with \(2^{p}\) , distinguishes the cases.
-
Explain, with one example from each of three different chapters, the difference between an ill-conditioned problem and an unstable algorithm, and in each case state what a user should do about it.
Hint
Pair each ill-conditioned problem with a stable algorithm that cannot help, and each unstable algorithm with a well-conditioned problem it ruins. Candidates: the logarithm near \(1\) against \(1-\cos x\) (Chapter 1); a triple root against Newton far from a root (Chapters 2 and 3); the Hilbert system against elimination without pivoting (Chapter 7); least squares with nearly dependent columns against the normal equations (Chapter 5); a stiff equation, which is well conditioned, against explicit Euler (Chapter 9); a boundary-value problem with a growing mode against shooting (Chapter 10). For conditioning the user changes the problem or the data; for stability, the algorithm.