Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 38
Текст из файла (страница 38)
(New York:McGraw-Hill), §4.10–2.4.4 Improper IntegralsFor our present purposes, an integral will be “improper” if it has any of thefollowing problems:• its integrand goes to a finite limiting value at finite upper and lower limits,but cannot be evaluated right on one of those limits (e.g., sin x/x at x = 0)• its upper limit is ∞ , or its lower limit is −∞• it has an integrable singularity at either limit (e.g., x−1/2 at x = 0)• it has an integrable singularity at a known place between its upper andlower limits• it has an integrable singularity at an unknown place between its upperand lower limits'∞If an integral is infinite (e.g., 1 x−1 dx), or does not exist in a limiting sense'∞(e.g., −∞ cos xdx), we do not call it improper; we call it impossible.
No amount ofclever algorithmics will return a meaningful answer to an ill-posed problem.In this section we will generalize the techniques of the preceding two sectionsto cover the first four problems on the above list. A more advanced discussion ofquadrature with integrable singularities occurs in Chapter 18, notably §18.3. Thefifth problem, singularity at unknown location, can really only be handled by theuse of a variable stepsize differential equation integration routine, as will be givenin Chapter 16.We need a workhorse like the extended trapezoidal rule (equation 4.1.11), butone which is an open formula in the sense of §4.1, i.e., does not require the integrandto be evaluated at the endpoints.
Equation (4.1.19), the extended midpoint rule, isthe best choice. The reason is that (4.1.19) shares with (4.1.11) the “deep” property142Chapter 4.Integration of Functionsof having an error series that is entirely even in h. Indeed there is a formula, not aswell known as it ought to be, called the Second Euler-Maclaurin summation formula,&xNx1f(x)dx = h[f3/2 + f5/2 + f7/2 + · · · + fN−3/2 + fN−1/2 ]B2 h2 (fN − f1 ) + · · ·4B2k h2k(2k−1)(2k−1)(1 − 2−2k+1 )(fN− f1) +···+(2k)!+(4.4.1)This equation can be derived by writing out (4.2.1) with stepsize h, then writing itout again with stepsize h/2, then subtracting the first from twice the second.It is not possible to double the number of steps in the extended midpoint ruleand still have the benefit of previous function evaluations (try it!).
However, it ispossible to triple the number of steps and do so. Shall we√ do this, or double andaccept the loss? On the average, tripling does a factor 3 of unnecessary work,since the “right” number of steps for a desired accuracy criterion may in fact fallanywhere√ in the logarithmic interval implied by tripling. For doubling, the factoris only 2, but we lose an extra factor of 2 in being unable to use all the previousevaluations. Since 1.732 < 2 × 1.414, it is better to triple.Here is the resulting routine, which is directly comparable to trapzd.#define FUNC(x) ((*func)(x))float midpnt(float (*func)(float), float a, float b, int n)This routine computes the nth stage of refinement of an extended midpoint rule.
func is inputas a pointer to the function to be integrated between limits a and b, also input. When called with'n=1, the routine returns the crudest estimate of ab f (x)dx. Subsequent calls with n=2,3,...(in that sequential order) will improve the accuracy of s by adding (2/3) × 3n-1 additionalinterior points. s should not be modified between sequential calls.{float x,tnm,sum,del,ddel;static float s;int it,j;if (n == 1) {return (s=(b-a)*FUNC(0.5*(a+b)));} else {for(it=1,j=1;j<n-1;j++) it *= 3;tnm=it;del=(b-a)/(3.0*tnm);ddel=del+del;The added points alternate in spacing betweenx=a+0.5*del;del and ddel.sum=0.0;for (j=1;j<=it;j++) {sum += FUNC(x);x += ddel;sum += FUNC(x);x += del;}s=(s+(b-a)*sum/tnm)/3.0;The new sum is combined with the old integralreturn s;to give a refined integral.}}4.4 Improper Integrals143The routine midpnt can exactly replace trapzd in a driver routine like qtrap(§4.2); one simply changes trapzd(func,a,b,j) to midpnt(func,a,b, j), andperhaps also decreases the parameter JMAX since 3JMAX−1 (from step tripling) is amuch larger number than 2JMAX−1 (step doubling).The open formula implementation analogous to Simpson’s rule (qsimp in §4.2)substitutes midpnt for trapzd and decreases JMAX as above, but now also changesthe extrapolation step to bes=(9.0*st-ost)/8.0;since, when the number of steps is tripled, the error decreases to 1/9th its size, not1/4th as with step doubling.Either the modified qtrap or the modified qsimp will fix the first problemon the list at the beginning of this section.
Yet more sophisticated is to generalizeRomberg integration in like manner:#include <math.h>#define EPS 1.0e-6#define JMAX 14#define JMAXP (JMAX+1)#define K 5float qromo(float (*func)(float), float a, float b,float (*choose)(float(*)(float), float, float, int))Romberg integration on an open interval. Returns the integral of the function func from a to b,using any specified integrating function choose and Romberg’s method. Normally choose willbe an open formula, not evaluating the function at the endpoints. It is assumed that choosetriples the number of steps on each call, and that its error series contains only even powers ofthe number of steps. The routines midpnt, midinf, midsql, midsqu, midexp, are possiblechoices for choose.
The parameters have the same meaning as in qromb.{void polint(float xa[], float ya[], int n, float x, float *y, float *dy);void nrerror(char error_text[]);int j;float ss,dss,h[JMAXP+1],s[JMAXP];h[1]=1.0;for (j=1;j<=JMAX;j++) {s[j]=(*choose)(func,a,b,j);if (j >= K) {polint(&h[j-K],&s[j-K],K,0.0,&ss,&dss);if (fabs(dss) <= EPS*fabs(ss)) return ss;}h[j+1]=h[j]/9.0;This is where the assumption of step tripling and an even}error series is used.nrerror("Too many steps in routing qromo");return 0.0;Never get here.}Don’t be put off by qromo’s complicated ANSI declaration. A typical invocation(integrating the Bessel function Y0 (x) from 0 to 2) is simply#include "nr.h"float answer;...answer=qromo(bessy0,0.0,2.0,midpnt);144Chapter 4.Integration of FunctionsThe differences between qromo and qromb (§4.3) are so slight that it is perhapsgratuitous to list qromo in full.
It, however, is an excellent driver routine for solvingall the other problems of improper integrals in our first list (except the intractablefifth), as we shall now see.The basic trick for improper integrals is to make a change of variables toeliminate the singularity, or to map an infinite range of integration to a finite one.For example, the identity&&b1/af(x)dx =a1/b 11dtft2tab > 0(4.4.2)can be used with either b → ∞ and a positive, or with a → −∞ and b negative, andworks for any function which decreases towards infinity faster than 1/x2 .You can make the change of variable implied by (4.4.2) either analytically andthen use (e.g.) qromo and midpnt to do the numerical evaluation, or you can letthe numerical algorithm make the change of variable for you. We prefer the lattermethod as being more transparent to the user.
To implement equation (4.4.2) wesimply write a modified version of midpnt, called midinf, which allows b to beinfinite (or, more precisely, a very large number on your particular machine, suchas 1 × 1030), or a to be negative and infinite.#define FUNC(x) ((*funk)(1.0/(x))/((x)*(x)))Effects the change of variable.float midinf(float (*funk)(float), float aa, float bb, int n)This routine is an exact replacement for midpnt, i.e., returns the nth stage of refinement ofthe integral of funk from aa to bb, except that the function is evaluated at evenly spacedpoints in 1/x rather than in x.
This allows the upper limit bb to be as large and positive asthe computer allows, or the lower limit aa to be as large and negative, but not both. aa andbb must have the same sign.{float x,tnm,sum,del,ddel,b,a;static float s;int it,j;b=1.0/aa;These two statements change the limits of integration.a=1.0/bb;if (n == 1) {From this point on, the routine is identical to midpnt.return (s=(b-a)*FUNC(0.5*(a+b)));} else {for(it=1,j=1;j<n-1;j++) it *= 3;tnm=it;del=(b-a)/(3.0*tnm);ddel=del+del;x=a+0.5*del;sum=0.0;for (j=1;j<=it;j++) {sum += FUNC(x);x += ddel;sum += FUNC(x);x += del;}return (s=(s+(b-a)*sum/tnm)/3.0);}}1454.4 Improper IntegralsIf you need to integrate from a negative lower limit to positive infinity, you dothis by breaking the integral into two pieces at some positive value, for example,answer=qromo(funk,-5.0,2.0,midpnt)+qromo(funk,2.0,1.0e30,midinf);Where should you choose the breakpoint? At a sufficiently large positive value sothat the function funk is at least beginning to approach its asymptotic decrease tozero value at infinity.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
