5.6 Quadratic and Cubic Equations 183 The solution in such cases is to use an alternative Clenshaw recurrence that incorporates ck's in an upward direction.The relevant equations are y-2=y-1=0 (5.5.25) 1 张=k+1,可lk-2-ak,躲-1-c (k=0,1,,N-1) (5.5.26) f(x)=cvFv(z)-(N,x)FN-1(x)N-1-Fw(z)yN-2(5.5.27) The rare case where equations (5.5.25)-(5.5.27)should be used instead of 81 equations (5.5.21)and (5.5.23)can be detected automatically by testing whether the operands in the first sum in(5.5.23)are opposite in sign and nearly equal in magnitude.Other than in this special case,Clenshaw's recurrence is always stable, 餐 independent of whether the recurrence for the functions F is stable in the upward or downward direction. RECIPESI CITED REFERENCES AND FURTHER READING: 2 Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- matics Series,Volume 55 (Washington:National Bureau of Standards;reprinted 1968 by Dover Publications,New York).pp.xii,697.[1] Press. Gautschi,W.1967,S/AM Review,vol.9,pp.24-82.[2] Lakshmikantham,V.,and Trigiante,D.1988,Theory of Difference Equations:Numerical Methods and Applications (San Diego:Academic Press).[3] Acton,F.S.1970,Numerica/Methods That Work,1990,corrected edition (Washington:Mathe- matical Association of America),pp.20ff.[4] SCIENTIFIC Clenshaw,C.W.1962,Mathematica/Tables,vol.5,National Physical Laboratory (London:H.M. Stationery Office).[5] 6 Dahlquist,G..and Bjorck,A.1974,Numerica/Methods (Englewood Cliffs,NJ:Prentice-Hall). 84.4.3,p.111. Goodwin,E.T.(ed.)1961,Modern Computing Methods,2nd ed.(New York:Philosophical Li- brary),p.76. 10-521 5.6 Quadratic and Cubic Equations Numerical Recipes 43106 The roots of simple algebraic equations can be viewed as being functions of the (outside equations'coefficients.We are taught these functions in elementary algebra.Yet. North Software. surprisingly many people don't know the right way to solve a quadratic equation with two real roots,or to obtain the roots of a cubic equation. There are two ways to write the solution of the quadratic equation az2+bx+c=0 (5.6.1) with real coefficients a,b,c,namely x=-b±vP-4ac (5.6.2) 2a
5.6 Quadratic and Cubic Equations 183 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The solution in such cases is to use an alternative Clenshaw recurrence that incorporates ck’s in an upward direction. The relevant equations are y−2 = y−1 =0 (5.5.25) yk = 1 β(k + 1, x) [yk−2 − α(k, x)yk−1 − ck], (k = 0, 1,...,N − 1) (5.5.26) f(x) = cN FN (x) − β(N, x)FN−1(x)yN−1 − FN (x)yN−2 (5.5.27) The rare case where equations (5.5.25)–(5.5.27) should be used instead of equations (5.5.21) and (5.5.23) can be detected automatically by testing whether the operands in the first sum in (5.5.23) are opposite in sign and nearly equal in magnitude. Other than in this special case, Clenshaw’s recurrence is always stable, independent of whether the recurrence for the functions Fk is stable in the upward or downward direction. CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), pp. xiii, 697. [1] Gautschi, W. 1967, SIAM Review, vol. 9, pp. 24–82. [2] Lakshmikantham, V., and Trigiante, D. 1988, Theory of Difference Equations: Numerical Methods and Applications (San Diego: Academic Press). [3] Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), pp. 20ff. [4] Clenshaw, C.W. 1962, Mathematical Tables, vol. 5, National Physical Laboratory (London: H.M. Stationery Office). [5] Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §4.4.3, p. 111. Goodwin, E.T. (ed.) 1961, Modern Computing Methods, 2nd ed. (New York: Philosophical Library), p. 76. 5.6 Quadratic and Cubic Equations The roots of simple algebraic equations can be viewed as being functions of the equations’ coefficients. We are taught these functions in elementary algebra. Yet, surprisingly many people don’t know the right way to solve a quadratic equation with two real roots, or to obtain the roots of a cubic equation. There are two ways to write the solution of the quadratic equation ax2 + bx + c =0 (5.6.1) with real coefficients a, b, c, namely x = −b ± √ b2 − 4ac 2a (5.6.2)
184 Chapter 5. Evaluation of Functions and 2c T= (5.6.3) -b士V2-4ac If you use either(5.6.2)or(5.6.3)to get the two roots,you are asking for trouble: If either a or c(or both)are small,then one of the roots will involve the subtraction of b from a very nearly equal quantity (the discriminant);you will get that root very inaccurately.The correct way to compute the roots is 1 q三一2 b+sgn(b)v62 -4ac (5.6.4) Then the two roots are 乃÷9 8合 and 工2=一 (5.6.5) If the coefficients a,b,c,are complex rather than real,then the above formulas RECIPES still hold,except that in equation(5.6.4)the sign of the square root should be 、远 2 chosen so as to make Re(b*Vb2-4ac)≥0 (5.6.6) Press. where Re denotes the real part and asterisk denotes complex conjugation. 9 Programs Apropos of quadratic equations,this seems a convenient place to recall that IENTIFIC( the inverse hyperbolic functions sinh-and cosh-1 are in fact just logarithms of solutions to such equations, 6 sinh-(x)=In(x+Vz2+1) (5.6.7) cosh-(x)=+In(x+Vz2-1) (5.6.8) Equation(5.6.7)is numerically robust for z >0.For negative z,use the symmetry Numerical 10621 sinh(-)=-sinh().Equation (5.6.8)is of course valid only for 21. 43126 For the cubic equation z3+az2+bx+c=0 (5.6.9) with real or complex coefficients a,b,c,first compute Q .a2-36 andR= 2a3-9ab+27c (5.6.10) 9 54 If Q and R are real (always true when a,b,c are real)and R2<Q3,then the cubic equation has three real roots.Find them by computing 0 arccos(R/VQ3) (5.6.11)
184 Chapter 5. Evaluation of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). and x = 2c −b ± √ b2 − 4ac (5.6.3) If you use either (5.6.2) or (5.6.3) to get the two roots, you are asking for trouble: If either a or c (or both) are small, then one of the roots will involve the subtraction of b from a very nearly equal quantity (the discriminant); you will get that root very inaccurately. The correct way to compute the roots is q ≡ −1 2 b + sgn(b) b2 − 4ac (5.6.4) Then the two roots are x1 = q a and x2 = c q (5.6.5) If the coefficients a, b, c, are complex rather than real, then the above formulas still hold, except that in equation (5.6.4) the sign of the square root should be chosen so as to make Re(b* b2 − 4ac) ≥ 0 (5.6.6) where Re denotes the real part and asterisk denotes complex conjugation. Apropos of quadratic equations, this seems a convenient place to recall that the inverse hyperbolic functions sinh−1 and cosh−1 are in fact just logarithms of solutions to such equations, sinh−1(x) = ln x + x2 + 1 (5.6.7) cosh−1 (x) = ± ln x + x2 − 1 (5.6.8) Equation (5.6.7) is numerically robust for x ≥ 0. For negative x, use the symmetry sinh−1(−x) = − sinh−1(x). Equation (5.6.8) is of course valid only for x ≥ 1. For the cubic equation x3 + ax2 + bx + c =0 (5.6.9) with real or complex coefficients a, b, c, first compute Q ≡ a2 − 3b 9 and R ≡ 2a3 − 9ab + 27c 54 (5.6.10) If Q and R are real (always true when a, b, c are real) and R2 < Q3, then the cubic equation has three real roots. Find them by computing θ = arccos(R/Q3) (5.6.11)
5.6 Quadratic and Cubic Equations 185 in terms of which the three roots are E1=-2VQcos a 3-3 E2=-2VQcos 0+2π - (5.6.12) 3 3 =-2VQcos 8-2m 、3 - (This equation first appears in Chapter VI of Francois Viete's treatise "De emen- datione,"published in 1615!) Otherwise,compute 菲 A--R+VR-Q35 (5.6.13) where the sign of the square root is chosen to make 2d州 Re(R*VR2-Q3)≥0 (5.6.14) (asterisk again denoting complex conjugation).If Q and R are both real,equations (5.6.13)-(5.6.14)are equivalent to A=-sgn(R) (5.6.15) 9 where the positive square root is assumed.Next compute B={8A (A≠0) (A=0) (5.6.16) IENTIFIC( in terms of which the three roots are 1=4+B)-号 (5.6.17) (the single real root when a,b,c are real)and =1+-营+ 1 2(A-B) 10521 (5.6.18) Numerical 有=一4+-景-9A- 1 431 (outside Recipes (in that same case,a complex conjugate pair).Equations(5.6.13)-(5.6.16)are arranged both to minimize roundofferror,and also(as pointed out by A.J.Glassman) to ensure that no choice of branch for the complex cube root can result in the North spurious loss of a distinct root. If you need to solve many cubic equations with only slightly different coeffi- cients,it is more efficient to use Newton's method(89.4). CITED REFERENCES AND FURTHER READING: Weast,R.C.(ed.)1967.Handbook of Tables for Mathematics,3rd ed.(Cleveland:The Chemical Rubber Co.),pp.130-133. Pachner,J.1983.Handbook of Numerical Analysis Applications(New York:McGraw-Hill),$6.1. McKelvey,J.P.1984,American Journal of Physics,vol.52,pp.269-270;see also vol.53,p.775. and vol..55,pp.374-375
5.6 Quadratic and Cubic Equations 185 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). in terms of which the three roots are x1 = −2 Q cos θ 3 − a 3 x2 = −2 Q cos θ + 2π 3 − a 3 x3 = −2 Q cos θ − 2π 3 − a 3 (5.6.12) (This equation first appears in Chapter VI of Fran¸cois Viete’s treatise “De emen- ` datione,” published in 1615!) Otherwise, compute A = − R + R2 − Q3 1/3 (5.6.13) where the sign of the square root is chosen to make Re(R* R2 − Q3) ≥ 0 (5.6.14) (asterisk again denoting complex conjugation). If Q and R are both real, equations (5.6.13)–(5.6.14) are equivalent to A = −sgn(R) |R| + R2 − Q3 1/3 (5.6.15) where the positive square root is assumed. Next compute B = Q/A (A = 0) 0 (A = 0) (5.6.16) in terms of which the three roots are x1 = (A + B) − a 3 (5.6.17) (the single real root when a, b, c are real) and x2 = −1 2 (A + B) − a 3 + i √3 2 (A − B) x3 = −1 2 (A + B) − a 3 − i √3 2 (A − B) (5.6.18) (in that same case, a complex conjugate pair). Equations (5.6.13)–(5.6.16) are arranged both to minimize roundoff error, and also (as pointed out by A.J. Glassman) to ensure that no choice of branch for the complex cube root can result in the spurious loss of a distinct root. If you need to solve many cubic equations with only slightly different coeffi- cients, it is more efficient to use Newton’s method (§9.4). CITED REFERENCES AND FURTHER READING: Weast, R.C. (ed.) 1967, Handbook of Tables for Mathematics, 3rd ed. (Cleveland: The Chemical Rubber Co.), pp. 130–133. Pachner, J. 1983, Handbook of Numerical Analysis Applications (New York: McGraw-Hill), §6.1. McKelvey, J.P. 1984, American Journal of Physics, vol. 52, pp. 269–270; see also vol. 53, p. 775, and vol. 55, pp. 374–375.
186 Chapter 5.Evaluation of Functions 5.7 Numerical Derivatives Imagine that you have a procedure which computes a function f(),and now you want to compute its derivative f'().Easy,right?The definition of the derivative,the limit as h-0 of f'()≈f+)-f) (5.7.1) h practically suggests the program:Pick a small value h;evaluate f(x+h);you probably have f()already evaluated,but if not,do it too;finally apply equation (5.7.1).What more needs to be said? 鱼 Quite a lot,actually.Applied uncritically,the above procedure is almost 二 guaranteed to produce inaccurate results.Applied properly,it can be the right way to compute a derivative only when the function f is fiercely expensive to compute, when you already have invested in computing f(z),and when,therefore,you want RECIPES to get the derivative in no more than a single additional function evaluation.In such a situation,the remaining issue is to choose h properly,an issue we now discuss: 言O。 9 There are two sources of error in equation(5.7.1).truncation error and roundoff error.The truncation error comes from higher terms in the Taylor series expansion, fe+)=f回+hf回+fa回)+f”a)+ (5.7.2) 。2 whence fe+)-f@=f+hf”+… 1 h (5.7.3) The roundoff error has various contributions.First there is roundoff error in h: Suppose,by way of an example,that you are at a point x=10.3 and you blindly choose h =0.0001.Neither x 10.3 nor x+h =10.30001 is a number with an exact representation in binary;each is therefore represented with some fractional error characteristic of the machine's floating-point format,em,whose value in single Recipes 10621 precision may be~10-7.The error in the effective value of h,namely the difference Numerica between x+h and x as represented in the machine,is therefore on the order of em, 4310 which implies a fractional error in h oforder~em/h10-2!By equation (5.7.1) Recipes this immediately implies at least the same large fractional error in the derivative. We arrive at Lesson 1:Always choose h so that +h and z differ by an exactly representable number.This can usually be accomplished by the program steps North temp=x+h (5.7.4) h temp- Some optimizing compilers,and some computers whose floating-point chips have higher internal accuracy than is stored externally,can foil this trick;if so,it is usually enough to declare temp as volatile,or else to call a dummy function donothing(temp)between the two equations(5.7.4).This forces temp into and out of addressable memory
186 Chapter 5. Evaluation of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 5.7 Numerical Derivatives Imagine that you have a procedure which computes a function f(x), and now you want to compute its derivative f (x). Easy, right? The definition of the derivative, the limit as h → 0 of f (x) ≈ f(x + h) − f(x) h (5.7.1) practically suggests the program: Pick a small value h; evaluate f(x + h); you probably have f(x) already evaluated, but if not, do it too; finally apply equation (5.7.1). What more needs to be said? Quite a lot, actually. Applied uncritically, the above procedure is almost guaranteed to produce inaccurate results. Applied properly, it can be the right way to compute a derivative only when the function f is fiercely expensive to compute, when you already have invested in computing f(x), and when, therefore, you want to get the derivative in no more than a single additional function evaluation. In such a situation, the remaining issue is to choose h properly, an issue we now discuss: There are two sources of error in equation (5.7.1), truncation error and roundoff error. The truncation error comes from higher terms in the Taylor series expansion, f(x + h) = f(x) + hf (x) + 1 2 h2f(x) + 1 6 h3f(x) + ··· (5.7.2) whence f(x + h) − f(x) h = f + 1 2 hf + ··· (5.7.3) The roundoff error has various contributions. First there is roundoff error in h: Suppose, by way of an example, that you are at a point x = 10.3 and you blindly choose h = 0.0001. Neither x = 10.3 nor x + h = 10.30001 is a number with an exact representation in binary; each is therefore represented with some fractional error characteristic of the machine’s floating-point format, m, whose value in single precision may be ∼ 10−7. The error in the effective value of h, namely the difference between x + h and x as represented in the machine, is therefore on the order of mx, which implies a fractional error in h of order ∼ mx/h ∼ 10−2! By equation (5.7.1) this immediately implies at least the same large fractional error in the derivative. We arrive at Lesson 1: Always choose h so that x+ h and x differ by an exactly representable number. This can usually be accomplished by the program steps temp = x + h h = temp − x (5.7.4) Some optimizing compilers, and some computers whose floating-point chips have higher internal accuracy than is stored externally, can foil this trick; if so, it is usually enough to declare temp as volatile, or else to call a dummy function donothing(temp) between the two equations (5.7.4). This forces temp into and out of addressable memory.
