f5-6

178 Chapter5. EvaluationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).5.6 Quadratic and Cubic Equations Therootsofsimplealgebraicequationscanbeviewedasbeingfunctionsofthe 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 ) 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 aorc(or both) are small, then one of the roots will involvethe subtraction ofbfroma verynearlyequal quantity(the discriminant);you will get that rootvery inaccurately. The correct way to compute the roots is q≡−1 2/bracketleftBig b+sgn(b)/radicalbig b2−4ac/bracketrightBig (5.6.4 ) Then the two roots are x1=q aand x2=c q(5.6.5 ) If the coefficients a, b, c, are complexrather than real, then the aboveformulas still hold, except that in equation (5.6.4) the sign of the square root should be chosen so as to make Re(b*/radicalbig 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−1and cosh−1are in fact just logarithms of solutions to such equations, sinh−1(x)= l n/parenleftbig x+/radicalbig x2+1/parenrightbig (5.6.7 ) cosh−1(x)=±ln/parenleftbig x+/radicalbig x2−1/parenrightbig (5.6.8 ) Equation (5.6.7)is numericallyrobust 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. Since FORTRAN mysteriously omits the inverse hyperbolic functions from its list of intrinsic functions, equations (5.6.7)–(5.6.8)are sometimes quite essential. 5.6QuadraticandCubicEquations 179Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).For thecubic equation x3+ax2+bx+c=0 ( 5.6.9 ) with real or complex coefficients a, b, c, first compute Q≡a2−3b 9and R≡2a3−9ab+2 7c 54(5.6.10 ) IfQandRare real (always true when a, b, care real)andR2<Q3, then the cubic equation has three real roots. Find them by computing θ=arccos (R//radicalbig Q3)( 5.6.11 ) in terms of which the three roots are x1=−2/radicalbig Qcos/parenleftbiggθ 3/parenrightbigg −a 3 x2=−2/radicalbig Qcos/parenleftbiggθ+2π 3/parenrightbigg −a 3 x3=−2/radicalbig Qcos/parenleftbiggθ−2π 3/parenrightbigg −a 3(5.6.12 ) (This equation first appears in Chapter VI of Fran¸ cois Vi`ete’s treatise “De emen- datione,” published in 1615!) Otherwise, compute A=−/bracketleftBig R+/radicalbig R2−Q3/bracketrightBig1/3 (5.6.13 ) where the sign of the square root is chosen to make Re(R*/radicalbig R2−Q3)≥0( 5.6.14 ) (asterisk again denoting complex conjugation). If QandRare both real, equations (5.6.13)–(5.6.14) are equivalent to A=−sgn(R)/bracketleftBig |R|+/radicalbig R2−Q3/bracketrightBig1/3 (5.6.15 ) where the positive square root is assumed. Next compute B=/braceleftbigg Q/A (A/negationslash=0 ) 0( A=0 )(5.6.16 ) in terms of which the three roots are x1=(A+B)−a 3(5.6.17 ) 180 Chapter5. EvaluationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).(the single real root when a, b, care 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 arrangedbothtominimizeroundofferror,andalso(aspointedoutbyA.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, HandbookofTablesforMathematics , 3rded.(Cleveland:TheChemical Rubber Co.), pp. 130–133. Pachner, J. 1983, Handbookof NumericalAnalysis Applications (NewYork: McGraw-Hill), §6.1. McKelvey,J.P.1984, AmericanJournalofPhysics ,vol.52,pp.269–270;seealsovol.53,p.775, and vol. 55, pp. 374–375. 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/prime(x). Easy, right? The definition of the derivative, the limit as h→ 0of f/prime(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 fisfiercelyexpensive to compute, when you already have invested in computing f(x), and when, therefore, you want to get thederivativein nomorethan a singleadditionalfunctionevaluation. Insuch a situation, the remainingissue is to choose hproperly,an issue we nowdiscuss: Therearetwosourcesoferrorinequation(5.7.1),truncationerrorandroundoff error. Thetruncationerrorcomesfromhighertermsin theTaylorseries expansion, f(x+h)=f(x)+hf/prime(x)+1 2h2f/prime/prime(x)+1 6h3f/prime/prime/prime(x)+··· (5.7.2 ) whence f(x+h)−f(x) h=f/prime+1 2hf/prime/prime+··· (5.7.3 )