Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 8
Текст из файла (страница 8)
This analysissuggests a remedy: for small θ, expand the coefficients in Taylor series and deal withthe cancellation and small divisor analytically. The resulting series are1.3 CASE STUDY 127It might be remarked that it was easy to compute these expansions by means of thesymbolic capabilities of the Student Edition of MATLAB. In the program used to compute the integral of Case Study 3, these expressions were used for θ < 0.1. Because theterms decrease rapidly, nested multiplication is not only an efficient way to evaluatethe expressions but is also accurate.As a numerical illustration of the difficulty we evaluated both forms of a for arange of θ in single precision in FORTRAN.
Reference values were computed usingthe trigonometric form and double precision. This must be done with some care. Forinstance, if T is a single precision variable and we want a double precision copy DTfor computing the reference values, the lines of codeT = 0.lE0DT = 0. 1D0are not equivalent toT = 0.1E0DT=TThis is because on a machine with binary or hexadecimal arithmetic, 0. 1E 0 agrees with0.1 D0 only to single precision. For the reference computation we require a double precision version of the actual machine number used in the single precision computations,hence we must use the second code.
As we have remarked previously, most computerstoday perform intermediate computations in higher precision, despite specification ofthe precision of all quantities. With T, S, and C declared as single precision variables,we found remarkable differences in the result ofS = SIN(T)C = COS(T)ALPHA=(T**2+T*S*C-2EO*S**2)/T**3andALPHA=(T**2+T*SIN(T)*COS(T)-2EO*SIN(T)**2)/T**3differences that depended on the machine and compiler used.
On a PC with a Pentiumchip, the second code gave nearly full single precision accuracy. The first gave thepoor results that we expect of computations carried out entirely in single precision.The coefficient a was computed for a range of θ using the trigonometric definition and single precision arithmetic and its relative error computed using a referencevalue computed in double precision. Similarly the error of the value computed in single precision from the Taylor series was found. Plotted against θ in Figure 1.3 is therelative error for both methods (on a logarithmic scale). Single precision accuracycorresponds to about seven digits, so the Taylor series approach gives about all theaccuracy we could hope for, although for the largest value of θ it appears that anotherterm in the expansion would be needed to get full accuracy.
Obviously the trigonometric definition leads to a great loss of accuracy for “small” θ. Indeed, θ is not verysmall in an absolute sense here; rather, it is small considering its implications for thecost of evaluating the integral when the parameter ω is moderately large.28CHAPTER 1ERRORS AND FLOATING POINT ARITHMETICFigure 1.3 Error in series form (•) versus trig form (*) for a Filon coefficient.1.4 FURTHER READINGA very interesting and readable account of the interaction of the floating point numbersystem with the solution of quadratic equations has been given by Forsythe [6]. Henrici[8] gives another elementary treatment of floating point arithmetic that introduces theuseful idea of a statistical treatment of errors. To pursue the subject in depth, consultthe book Rounding Errors in Algebraic Processes by J.
H. Wilkinson [10]. Wilkinson’s books are unmatched for their blend of theoretical advance, striking examples,practical insight, applications, and readability. For more information on the practicalevaluation of special functions, see the books by Cody and Waite [3] or Fike [5]. Otherinteresting discussions on floating point arithmetic are the books of Goldberg [7] andHigham [9].REFERENCES1.
ANSI/IEEE, IEEE Standard for Binary Floating Point Arithmetic, Std 754-1985, New York,1985.2. R. Bate, D. Miller, and J. White, Fundamentals of Astrophysics, Dover, New York, 1971.3. W. Cody and W. Waite, Software Manual for the Elementary Functions, Prentice Hall, Englewood Cliffs, N.J., 1980.4. J. Dongarra, J. Bunch, C. Moler, and G. Stewart, LINPACK Users’ Guide, SIAM, Philadelphia,1979.5.
C. Fike, Computer Evaluation of Mathematical Functions, Prentice Hall, Englewood Cliffs,N.J., 1968.6. G. Forsythe, “What is a satisfactory quadratic equation solver?,” in Constructive Aspects of theFundamental Theorem of Algebra, B. Dejon and P. Henrici, eds., Wiley, London, 1969.REFERENCES 297. D. Goldberg, “What every computer scientist should know about floating-point arithmetic,”ACM Computing Surveys, 23 (1991), pp.
5-48.8. P. Henrici, Elements of Numerical Analysis, Wiley, New York, 1964.9. N. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, Philadelphia, 1996.10. J. Wilkinson, Rounding Errors in Algebraic Processes, Dover, Mineola, N.Y., 1994.MISCELLANEOUS EXERCISES FOR CHAPTER 11.15 Use three-digit decimal chopped arithmetic with m =- 100 and M = 100 to construct examples for whichare of great practical value. It appears to be necessaryto evaluate a large number of sines and cosines if wewish to evaluate such a series, but this can be donecheaply by recursion. For the specific x of interest, forn = 1, 2,...letsn = sinnx and cn = cosnx.You are allowed to use negative numbers. Examplescan be constructed so that either one of the expressions cannot be formed in the arithmetic, or both canbe formed but the values are different.1.16 For a set of measurements xl ,x2,. .
. ,xN, the samplemeanis defined to beThe sample standard deviation s is defined to beShow that for n = 2, 3,. . .sn =S 1 C n-I+ c 1 s n-l and c n = c 1 c n-l -s 1 s n-l .After evaluating sl = sinx and c1 = cosx with the intrinsic functions of the programming language, thisrecursion can be used to evaluate simply and inexpensively all the sinnx and cosnx that are needed.To see that the recursion is stable, suppose that forsome m > 1, sm and cm are computed incorrectly as= sm + εm and= cm +If no further arithmetic errors are made, the errors ε m andwill propagate in the recurrence so that we computeAnother expression,for n = m + l,.... Let εn andso that, by definition,is often recommended for hand computation of s.Show that these two expressions for s are mathematically equivalent.
Explain why one of them may provide better numerical results than the other, and construct an example to illustrate your point.Prove that for all n > mbe the errors inandwhich implies that for all n > m1.17 Fourier series,In this sense, errors are not amplified and the recurrence is quite stable.Previous Home NextCHAPTER 2SYSTEMS OF LINEAR EQUATIONSOne of the most frequently encountered problems in scientific computation is that ofsolving n simultaneous linear equations in n unknowns. If we denote the unknowns byx1,x2, .
. . xn, such a system can be written in the formThe given data here are the right-hand sides bi , i = 1, 2,. . . , n, and the coefficients aijfor i, j = 1, 2,..., n. Problems of this nature arise almost everywhere in the applicationsof mathematics (e.g., the fitting of polynomials and other curves through data andthe approximation of differential and integral equations by finite, algebraic systems).Several specific examples are found in the exercises for this chapter (see also [12]or [13]).
To talk about (2.1) conveniently, we shall on occasion use some notationfrom matrix theory. However, we do not presume that the reader has an extensivebackground in this area. Using matrices, (2.1) can be written compactly asAx=b,(2.2)whereIf a110, the equation has a unique solution, namely x1 = b1/a11. If a11 = 0, thensome problems do not have solutions (b1 0) while others have many solutions (ifb1 = 0, any number x1 is a solution). The same is true for general n. There are twokinds of matrices, nonsingular and singular.
If the matrix A is nonsingular, there is aunique solution vector x for any given right-hand side b. If A is singular, there is no3031solution for some right-hand sides b and many solutions for other b. In this book weconcentrate on systems of linear equations with nonsingular matrices.Example 2.1.The problem2 x15 x1+ 3x2 = 8+ 4x2 = 13has a nonsingular coefficient matrix.
The linear system has the unique solutionxExample 2.2.1=1,x2=2 or x=The problem2x14x1+ 3x2+ 6x2==47has a singular coefficient matrix. Ifthere is no solution, for if x1 and x2 were numbers such that 4 = 2x1 + 3x2, then wewould have 8 = 2 × 4 = 2 × (2x1+ 3x 2 ) = 4x1+ 6x2, which is impossible because ofthe second equation.
Ifthere are many solutions, namelynfor all real numbers c.In the nonsingular case there exists a matrix called the inverse of A, denoted byA-1, such that the unique solution of (2.2) is given byx = A-1 b.For n = 1, A-1 = (1/a 11). Should we compute A-1 and then form the product A-1 bto solve (2.2)? We shall see that the answer is generally no even if we want to solve(2.2) with the same matrix A and many different right-hand sides b.32CHAPTER 22.1GAUSSIAN ELIMINATION WITH PARTIAL PIVOTINGSYSTEMS OF LINEAR EQUATIONSThe most popular method for solving a nonsingular system of linear equations (2.1) iscalled Gaussian elimination.
It is both simple and effective. In principle it can be usedto compute solutions of problems with singular matrices when they have solutions, butthere are better ways to do this. The basic idea in elimination is to manipulate theequations of (2.1) so as to obtain an equivalent set of equations that is easy to solve.An equivalent set of equations is one that has the same solutions. There are three basicoperations used in elimination: (1) multiplying an equation by a nonzero constant,(2) subtracting a multiple of one equation from another, and (3) interchanging rows.First, if any equation of (2.1) is multiplied by the nonzero constant a, we obtain anequivalent set of equations. To see this, suppose that we multiply the k th equation bya to get(2.3)If x1, x2, .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
