208 Chapter 5.Evaluation of Functions Figure 5.13.I shows the discrepancies for the first five iterations of ratlsq when it is applied to find the m=k=4 rational fit to the function f(x)=cos/(1+e")in the interval (0,).One sees that after the first iteration,the results are virtually as good as the minimax solution.The iterations do not converge in the order that the figure suggests:In fact,it is the second iteration that is best (has smallest maximum deviation).The routine ratlsq accordingly returns the best of its iterations,not necessarily the last one;there is no advantage in doing more than five iterations. CITED REFERENCES AND FURTHER READING: Ralston,A.and Wilf,H.S.1960,Mathematical Methods for Digital Computers(New York:Wiley). Chapter 13.[1] 5.14 Evaluation of Functions by Path Cam Integration In computer programming,the technique of choice is not necessarily the most 9 efficient,or elegant,or fastest executing one.Instead,it may be the one that is quick to implement,general,and easy to check. Q>%,入 One sometimes needs only a few,or a few thousand,evaluations of a special function,perhaps a complex valued function of a complex variable,that has many 9 different parameters,or asymptotic regimes,or both.Use of the usual tricks(series, continued fractions,rational function approximations,recurrence relations,and so forth)may result in a patchwork program with tests and branches to different formulas.While such a program may be highly efficient in execution,it is often not the shortest way to the answer from a standing start. A different technique of considerable generality is direct integration of a function's defining differential equation-an ab initio integration for each desired function value-along a path in the complex plane if necessary.While this may at first seem like swatting a fly with a golden brick,it turns out that when you already have the brick,and the fly is asleep right under it,all you have to do is let it fall! 10621 As a specific example,let us consider the complex hypergeometric func- Numerica tion 2F1(a,b,c;z),which is defined as the analytic continuation of the so-called 431 hypergeometric series, Recipes 2(a,6,c2)=1+驰名+aa+b+2 c(c+1)2z万+… North +a(a+1)..(a+j-1b6+1)..6+j-1)z c(c+1)..(c+j-1) (5.14.1) The series converges only within the unit circle z<1 (see [1]).but one's interest in the function is often not confined to this region. The hypergeometric function 2F is a solution(in fact the solution that is regular at the origin)of the hypergeometric differential equation,which we can write as z(1-2)F"=abF-[c-(a+b+1)zF' (5.14.2)
208 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). Figure 5.13.1 shows the discrepancies for the first five iterations of ratlsq when it is applied to find the m = k = 4 rational fit to the function f(x) = cos x/(1 + ex) in the interval (0, π). One sees that after the first iteration, the results are virtually as good as the minimax solution. The iterations do not converge in the order that the figure suggests: In fact, it is the second iteration that is best (has smallest maximum deviation). The routine ratlsq accordingly returns the best of its iterations, not necessarily the last one; there is no advantage in doing more than five iterations. CITED REFERENCES AND FURTHER READING: Ralston, A. and Wilf, H.S. 1960, Mathematical Methods for Digital Computers (New York: Wiley), Chapter 13. [1] 5.14 Evaluation of Functions by Path Integration In computer programming, the technique of choice is not necessarily the most efficient, or elegant, or fastest executing one. Instead, it may be the one that is quick to implement, general, and easy to check. One sometimes needs only a few, or a few thousand, evaluations of a special function, perhaps a complex valued function of a complex variable, that has many different parameters, or asymptotic regimes, or both. Use of the usual tricks (series, continued fractions, rational function approximations, recurrence relations, and so forth) may result in a patchwork program with tests and branches to different formulas. While such a program may be highly efficient in execution, it is often not the shortest way to the answer from a standing start. A different technique of considerable generality is direct integration of a function’s defining differential equation – an ab initio integration for each desired function value — along a path in the complex plane if necessary. While this may at first seem like swatting a fly with a golden brick, it turns out that when you already have the brick, and the fly is asleep right under it, all you have to do is let it fall! As a specific example, let us consider the complex hypergeometric function 2F1(a, b, c; z), which is defined as the analytic continuation of the so-called hypergeometric series, 2F1(a, b, c; z)=1+ ab c z 1! + a(a + 1)b(b + 1) c(c + 1) z2 2! + ··· + a(a + 1)...(a + j − 1)b(b + 1)...(b + j − 1) c(c + 1)...(c + j − 1) zj j! + ··· (5.14.1) The series converges only within the unit circle |z| < 1 (see [1]), but one’s interest in the function is often not confined to this region. The hypergeometric function 2F1 is a solution (in fact the solution that is regular at the origin) of the hypergeometric differential equation, which we can write as z(1 − z)F = abF − [c − (a + b + 1)z]F (5.14.2)
5.14 Evaluation of Functions by Path Integration 209 Here prime denotes d/dz.One can see that the equation has regular singular points at z=0,1,and oo.Since the desired solution is regular at z=0,the values 1 and oo will in general be branch points.If we want 2F to be a single valued function. we must have a branch cut connecting these two points.A conventional position for this cut is along the positive real axis from 1 to oo,though we may wish to keep open the possibility of altering this choice for some applications. Our golden brick consists of a collection of routines for the integration of sets of ordinary differential equations,which we will develop in detail later,in Chapter 16.For now,we need only a high-level,"black-box"routine that integrates such a set from initial conditions at one value of a(real)independent variable to final conditions at some other value of the independent variable,while automatically adjusting its internal stepsize to maintain some specified accuracy.That routine is called odeint and,in one particular invocation,calculates its individual steps with a sophisticated Bulirsch-Stoer technique. Suppose that we know values for F and its derivative Fat some value zo,and that we want to find F at some other point z in the complex plane.The straight-line path connecting these two points is parametrized by z(s)=20+s(21-20) (5.14.3) 9 with s a real parameter.The differential equation(5.14.2)can now be written as a set of two first-order equations, 9 dF 豆三 ds =(21-20)F dFr abF-[c-(a+b+1)z]F (5.14.4) ds =(21-20） 2(1-2) 6 to be integrated from s =0 to s 1.Here F and F are to be viewed as two independent complex variables.The fact that prime means d/dz can be ignored:it will emerge as a consequence of the first equation in(5.14.4).Moreover,the real and imaginary parts of equation(5.14.4)define a set of four real differential equations, 、p2 Numerica 10621 with independent variable s.The complex arithmetic on the right-hand side can be viewed as mere shorthand for how the four components are to be coupled.It is precisely this point of view that gets passed to the routine odeint,since it knows nothing of either complex functions or complex independent variables. It remains only to decide where to start,and what path to take in the complex plane,to get to an arbitrary point z.This is where consideration of the function's North singularities,and the adopted branch cut,enter.Figure 5.14.1 shows the strategy that we adopt.For 1,above and below the branch cut;the purpose of starting away from the real axis in these cases is to avoid passing too close to the singularity at =1(see Figure 5.14.1).The location of the branch cut is defined by the fact that our adopted strategy never integrates across the real axis for Re (z)>1. An implementation of this algorithm is given in $6.12 as the routine hypgeo
5.14 Evaluation of Functions by Path Integration 209 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). Here prime denotes d/dz. One can see that the equation has regular singular points at z = 0, 1, and ∞. Since the desired solution is regular at z = 0, the values 1 and ∞ will in general be branch points. If we want 2F1 to be a single valued function, we must have a branch cut connecting these two points. A conventional position for this cut is along the positive real axis from 1 to ∞, though we may wish to keep open the possibility of altering this choice for some applications. Our golden brick consists of a collection of routines for the integration of sets of ordinary differential equations, which we will develop in detail later, in Chapter 16. For now, we need only a high-level, “black-box” routine that integrates such a set from initial conditions at one value of a (real) independent variable to final conditions at some other value of the independent variable, while automatically adjusting its internal stepsize to maintain some specified accuracy. That routine is called odeint and, in one particular invocation, calculates its individual steps with a sophisticated Bulirsch-Stoer technique. Suppose that we know values for F and its derivative F at some value z0, and that we want to find F at some other point z1 in the complex plane. The straight-line path connecting these two points is parametrized by z(s) = z0 + s(z1 − z0) (5.14.3) with s a real parameter. The differential equation (5.14.2) can now be written as a set of two first-order equations, dF ds = (z1 − z0)F dF ds = (z1 − z0) abF − [c − (a + b + 1)z]F z(1 − z) (5.14.4) to be integrated from s = 0 to s = 1. Here F and F are to be viewed as two independent complex variables. The fact that prime means d/dz can be ignored; it will emerge as a consequence of the first equation in (5.14.4). Moreover, the real and imaginary parts of equation (5.14.4) define a set of four real differential equations, with independent variable s. The complex arithmetic on the right-hand side can be viewed as mere shorthand for how the four components are to be coupled. It is precisely this point of view that gets passed to the routine odeint, since it knows nothing of either complex functions or complex independent variables. It remains only to decide where to start, and what path to take in the complex plane, to get to an arbitrary point z. This is where consideration of the function’s singularities, and the adopted branch cut, enter. Figure 5.14.1 shows the strategy that we adopt. For |z| ≤ 1/2, the series in equation (5.14.1) will in general converge rapidly, and it makes sense to use it directly. Otherwise, we integrate along a straight line path from one of the starting points (±1/2, 0) or (0, ±1/2). The former choices are natural for 0 1, above and below the branch cut; the purpose of starting away from the real axis in these cases is to avoid passing too close to the singularity at z = 1 (see Figure 5.14.1). The location of the branch cut is defined by the fact that our adopted strategy never integrates across the real axis for Re (z) > 1. An implementation of this algorithm is given in §6.12 as the routine hypgeo.
210 Chapter 5. Evaluation of Functions Im : http://www.nr. use power series branch cut 83 Re granted for (including this one) 19881992 1-.200 ● 872 Cambridge from NUMERICAL RECIPES IN C: (North America server computer, users to make one paper University Press. THE ART Programs Figure 5.14.1.Complex plane showing the singular points of the hypergeometric function,its branch ictly proh cut,and some integration paths from the circle z=1/2 (where the power series converges rapidly) to other points in the plane. to dir A number of variants on the procedure described thus far are possible,and easy to program.Ifsuccessively called values of z are close together (with identical values OF SCIENTIFIC COMPUTING(ISBN of a,b,and c),then you can save the state vector(F,F)and the corresponding value 1988-19920 of z on each call,and use these as starting values for the next call.The incremental integration may then take only one or two steps.Avoid integrating across the branch 10-521 cut unintentionally:the function value will be "correct,"but not the one you want. Alternatively,you may wish to integrate to some position z by a dog-leg path that does cross the real axis Re z>1,as a means of moving the branch cut.For Numerical Recipes 43106 example,in some cases you might want to integrate from (0,1/2)to (3/2,1/2), and go from there to any point with Re z>1-with either sign of Im z.(If (outside you are,for example,finding roots of a function by an iterative method,you do North Software. not want the integration for nearby values to take different paths around a branch point.If it does,your root-finder will see discontinuous function values,and will likely not converge correctly!) visit website In any case,be aware that a loss of numerical accuracy can result if you integrate machine through a region of large function value on your way to a final answer where the function value is small.(For the hypergeometric function,a particular case of this is when a and b are both large and positive,with c and1.)In such cases,you'll need to find a better dog-leg path. The general technique of evaluating a function by integrating its differential equation in the complex plane can also be applied to other special functions.For
210 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). use power series branch cut Im 0 1 Re Figure 5.14.1. Complex plane showing the singular points of the hypergeometric function, its branch cut, and some integration paths from the circle |z| = 1/2 (where the power series converges rapidly) to other points in the plane. A number of variants on the procedure described thus far are possible, and easy to program. If successively called values of z are close together (with identical values of a, b, and c), then you can save the state vector(F, F ) and the corresponding value of z on each call, and use these as starting values for the next call. The incremental integration may then take only one or two steps. Avoid integrating across the branch cut unintentionally: the function value will be “correct,” but not the one you want. Alternatively, you may wish to integrate to some position z by a dog-leg path that does cross the real axis Re z > 1, as a means of moving the branch cut. For example, in some cases you might want to integrate from (0, 1/2) to (3/2, 1/2), and go from there to any point with Re z > 1 — with either sign of Im z. (If you are, for example, finding roots of a function by an iterative method, you do not want the integration for nearby values to take different paths around a branch point. If it does, your root-finder will see discontinuous function values, and will likely not converge correctly!) In any case, be aware that a loss of numerical accuracy can result if you integrate through a region of large function value on your way to a final answer where the function value is small. (For the hypergeometric function, a particular case of this is when a and b are both large and positive, with c and x >∼ 1.) In such cases, you’ll need to find a better dog-leg path. The general technique of evaluating a function by integrating its differential equation in the complex plane can also be applied to other special functions. For
5.14 Evaluation of Functions by Path Integration 211 example,the complex Bessel function.Airy function.Coulomb wave function,and Weber function are all special cases of the confluent hypergeometric fumnction,with a differential equation similar to the one used above (see,e.g.,[1]813.6,for a table of special cases).The confluent hypergeometric function has no singularities at finite z: That makes it easy to integrate.However,its essential singularity at infinity means that it can have,along some paths and for some parameters,highly oscillatory or exponentially decreasing behavior:That makes it hard to integrate. Some case by case judgment (or experimentation)is therefore required. http://www.nr Permission is read able files Sample page opyright(C) CITED REFERENCES AND FURTHER READING: Abramowitz,M.,and Stegun,I.A.1964.Handbook of Mathematical Functions.Applied Mathe- 83 matics Series,Volume 55(Washington:National Bureau of Standards:reprinted 1968 by granted for Dover Publications,New York).[1] call 1-800-872-7423(North America only).orsend email to directcustserv@cambridge.org(outside North America). 1988-1992 by Cambridge University Press.Programs Copyright(C)1988-1992 by Numerical Recipes from NUMERICAL RECIPES IN C:THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) Software
5.14 Evaluation of Functions by Path Integration 211 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). example, the complex Bessel function, Airy function, Coulomb wave function, and Weber function are all special cases of the confluent hypergeometric function, with a differential equation similar to the one used above (see, e.g., [1] §13.6, for a table of special cases). The confluent hypergeometric function has no singularities at finite z: That makes it easy to integrate. However, its essential singularity at infinity means that it can have, along some paths and for some parameters, highly oscillatory or exponentially decreasing behavior: That makes it hard to integrate. Some case by case judgment (or experimentation) is therefore required. 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). [1]