Heath - Scientific Computing (523150), страница 13
Текст из файла (страница 13)
. . , 20. If the programming languageyou use does not have an operator for exponentiation, you may use the equivalent formula(1 + 1/n)n = exp(n log(1 + 1/n)),where exp and log are built-in functions. Determine the error in your successive approximations by comparing them with the value ofexp(1). Does the error always decrease as nincreases? Explain your results.1.5 (a) Consider the functionf (x) = (ex − 1)/x.Use l’Hôpital’s rule to show thatlim f (x) = 1.x→0(b) Check this result empirically by writinga program to compute f (x) for x = 10−k ,32CHAPTER 1.
SCIENTIFIC COMPUTINGk = 1, . . . , 16. Do your results agree with theoretical expectations? Explain why.(c) Perform the experiment in part b again,this time using the mathematically equivalentformulationxxf (x) = (e − 1)/ log(e ),evaluated as indicated, with no simplification.If this works any better, can you explain why?1.6 Suppose you need to generate n + 1equally spaced points on the interval [a, b],with spacing h = (b − a)/n.(a) In floating-point arithmetic, which of thefollowing methods,x0 = a,xk = xk−1 + h,k = 1, .
. . , norxk = a + kh,k = 0, . . . , n,is better, and why?(b) Write a program implementing both methods and find an example, say, with a = 0 andb = 1, that illustrates the difference betweenthem.1.7 (a) Write a program to compute an approximate value for the derivative of a functionusing the finite-difference formulaf 0 (x) ≈f (x + h) − f (x).hTest your program using the function sin(x)for x = 1. Determine the error by comparing with the built-in function cos(x).
Plot themagnitude of the error as a function of h, forh = 21 , 14 , 18 , . . . . You should use a log scalefor h and for the magnitude of the error. Isthere a minimum value for the magnitude ofthe error? How does the corresponding valuefor h compare with the rule of thumbh≈√mach · |x| ?(b) Repeat the exercise using the centered difference approximationf 0 (x) ≈f (x + h) − f (x − h).2h1.8 Consider the infinite series∞X1.nn=1(a) Prove that the series is divergent. (Hint:Group the terms in sets containing terms1/(2k−1 + 1) down to 1/2k , for k = 1, 2, .
. . .)(b) Explain why summing the series infloating-point arithmetic yields a finite sum.(c) Try to predict when the partial sum willcease to change in both IEEE single-precisionand double-precision floating-point arithmetic.Given the execution rate of your computer forfloating-point operations, try to predict howlong each computation would take to complete.(d ) Write two programs to compute the sum ofthe series, one in single precision and the otherin double precision. Monitor the progress ofthe summation by printing out the index andpartial sum periodically. What stopping criterion should you use? What result is actually produced on your computer? Compareyour results with your predictions, includingthe execution time required.
(Caution: Yoursingle-precision version should terminate fairlyquickly, but your double-precision version maytake much longer, so it may not be practicalto run it to completion, even if your computerbudget is generous.)1.9 (a) Write a program to compute the exponential function ex using the infinite seriesex = 1 + x +x2x3++ ···.2!3!(b) Summing in the natural order, what stopping criterion should you use?(c) Test your program forx = ±1, ±5, ±10, ±15, ±20,and compare your results with the built-infunction exp(x).(d ) Can you use the series in this form to getaccurate results for x < 0? (Hint: e−x =1/ex .)(e) Can you rearrange the series or regroup theterms in any way to get more accurate resultsfor x < 0?COMPUTER PROBLEMS331.10 Write a program to solve the quadraticequation ax2 + bx + c = 0 using the standardquadratic formula√−b ± b2 − 4acx=2aor the alternative formulax=2c√.−b ∓ b2 − 4acYour program should accept values for the coefficients a, b, and c as input and produce thetwo roots of the equation as output.
Your program should detect when the roots are imaginary, but need not use complex arithmetic explicitly. You should guard against unnecessaryoverflow, underflow, and cancellation. Whenshould you use each of the two formulas? Tryto make your program robust when given unusual input values. Any root that is withinthe range of the floating-point system shouldbe computed accurately, even if the other isout of range. Test your program using the following values for the coefficients:a66 × 103001110−30b55 × 10301−105−4−1030c−4−4 × 1030113.9999991030ax3 + bx2 + cx + d = 0,where the coefficients are real and a 6= 0, hasat least one real root, which can be computedin closed form as follows.
Make the substitution y = x+b/(3a). Then the original equationbecomesy 3 + 3py + q = 0,p=andq=3ac − b29a227a2 d − 9abc + 2b3.27a3If we now takeα=−q +β=p4p3 + q 22−q −p4p3 + q 2,2then one real root of the original cubic equation is given byp√x = 3 α − 3 β.Write a routine using this method in real arithmetic to compute one real root of an arbitrarycubic equation given its (real) coefficients. Tryto make your routine as robust as possible,guarding against unnecessary overflow, underflow, and cancellation. What should your routine do if a = 0? Test your routine for variousvalues of the coefficients, analogous to thoseused in the previous exercise.1.12 (a) Write a program to compute themean x̄ and standard deviation σ of a finitesequence xi .
Your program should accept avector x of dimension n as input and producethe mean and standard deviation of the sequence as output. For the standard deviation,try both the two-pass formula"n1 X(xi − x̄)2σ=n − 1 i=1#1/2and the one-pass formula"1.11 A cubic equationwhereand1σ=n−1nXx2i − nx̄2!#1/2i=1and compare the results for an input sequenceof your choice.(b) Can you devise an input data sequence thatdramatically illustrates the numerical difference between these two mathematically equivalent formulas? (Caution: Beware of takingthe square root of a negative number.)1.13 If an amount a is invested at interest rater compounded n times per year, then the finalvalue f at the end of one year is given byf = a(1 + r/n)n .This is the familiar formula for compound interest.
With simple interest, n = 1. Typically,34CHAPTER 1. SCIENTIFIC COMPUTINGcompounding is done quarterly, n = 4, or perhaps even daily, n = 365. Obviously, the morefrequent the compounding, the greater the final amount, because more interest is paid onprevious interest. But how much differencedoes this frequency actually make? Write aprogram that implements the compound interest formula. Test your program using aninitial investment of a = 100, an interest rateof 5 percent (i.e., r = 0.05), and the followingvalues for n: 1, 4, 12, 365, 10,000, and 20,000.Implement the compound interest formula intwo different ways:(a) If the programming language you use doesnot have an operator for exponentiation (e.g.,C), then you might implement the compoundinterest formula using a loop that repeatedlymultiplies a by (1 + r/n) for a total of n times.Even if your programming language does havean operator for exponentiation (e.g., Fortran),try implementing the compound interest formula using such a loop and print your resultsfor the input values.
Does the final amount always grow with the frequency of compounding,as it should? Can you explain this behavior?(b) With the functions exp(x) and log(x), thecompound interest formula can also be writtenf = a exp(n log(1 + r/n)).Implement this formula using the corresponding built-in functions and compare your resultswith those for the first implementation usingthe loop, for the same input values.1.14 The polynomial (x − 1)6 has the valuezero at x = 1 and is positive elsewhere.
Theexpanded form of the polynomial, x6 − 6x5 +15x4 − 20x3 + 15x2 − 6x + 1, is mathematicallyequivalent but may not give the same resultsnumerically. Compute and plot the values ofthis polynomial, using each of the two forms,for 101 equally spaced points in the interval[0.995, 1.005], i.e., with a spacing of 0.0001.Your plot should be scaled so that the valuesfor x and for the polynomial use the full rangesof their respective axes. Can you explain thisbehavior?1.15 Write a program that sums n random,single-precision floating-point numbers xi , uniformly distributed on the interval [0, 1] (seeTable 13.1 for an appropriate random number generator). Sum the numbers in each ofthe following ways (use only single-precisionfloating-point variables unless specifically indicated otherwise):(a) Sum the numbers in the order in whichthey were generated, using a double-precisionvariable in which to accumulate the sum.(b) Sum the numbers in the order in whichthey were generated, this time using a singleprecision accumulator.(c) Use the following algorithm (due to Kahan), again using only single precision, to sumthe numbers in the order in which they weregenerated:s = x1c=0for i = 2 to ny = xi − ct=s+yc = (t − s) − ys=tend(d ) Sum the numbers in order of increasingmagnitude (this will require that the numbersbe sorted before summing, for which you mayuse a library sorting routine).(e) Sum the numbers in order of decreasingmagnitude (i.e., reverse the order of summation from part d ).Run your program for various values of n andcompare the results for methods a through e.You may need to use a fairly large value forn to see a substantial difference.
How do themethods rank in terms of accuracy, and why?How do the methods compare in cost? Canyou explain why the algorithm in part c works?1.16 Write a program to generate the first nterms in the sequence given by the differenceequationxk+1 = 2.25xk − 0.5xk−1 ,with starting valuesx1 =13and x2 =1.12Use n = 225 if you are working in single precision, n = 60 if you are working in doubleCOMPUTER PROBLEMSprecision. Make a semilog plot of the valuesyou obtain as a function of k.
The exact solution of the difference equation is given by41−k,3which decreases monotonically as k increases.Does your graph confirm this theoretically expected behavior? Can you explain your results? (Hint: Find the general solution to thedifference equation.)xk =1.17 Write a program to generate the first nterms in the sequence given by the differenceequation:xk+1 = 111 − (1130 − 3000/xk−1 )/xk ,with starting values6111and x2 =.211Use n = 10 if you are working in single precision, n = 20 if you are working in doublex1 =35precision. The exact solution is a monotonically increasing sequence converging to 6.
Canyou explain your results?1.18 The Euclidean norm of an n-dimensionalvector x is defined bykxk2 =nXx2i!1/2.i=1Implement a robust routine for computing thisquantity for any given input vector x. Yourroutine should avoid overflow and harmful underflow. Compare both the accuracy and performance of your robust routine with a morestraightforward naive implementation. Canyou devise a vector that produces significantlydifferent results from the two routines? Howmuch performance does the robust routine sacrifice?36CHAPTER 1.
SCIENTIFIC COMPUTINGChapter 2Systems of Linear Equations2.1Linear SystemsSystems of linear algebraic equations arise in almost every aspect of applied mathematicsand scientific computation. Such systems often occur naturally, but they are also frequentlythe result of approximating nonlinear equations by linear equations or differential equationsby algebraic equations. We will see many examples of such approximations throughoutthis book. For these reasons, the efficient and accurate solution of linear systems forms thecornerstone of many numerical methods for solving a wide variety of practical computationalproblems.In matrix-vector notation, a system of linear algebraic equations has the formAx = b,where A is an m × n matrix, b is a given m-vector, and x is the unknown solution n-vectorto be determined.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
