c4-4 (779481)
Текст из файла
4.4 Improper Integrals141which contain no singularities, and where the endpoints are also nonsingular. qromb,in such circumstances, takes many, many fewer function evaluations than either ofthe routines in §4.2. For example, the integralZ2x4 log(x +px2 + 1)dxconverges (with parameters as shown above) on the very first extrapolation, afterjust 5 calls to trapzd, while qsimp requires 8 calls (8 times as many evaluations ofthe integrand) and qtrap requires 13 calls (making 256 times as many evaluationsof the integrand).CITED REFERENCES AND FURTHER READING:Stoer, J., and Bulirsch, R.
1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§§3.4–3.5.Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),§§7.4.1–7.4.2.Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (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 limitsR∞If an integral is infinite (e.g., 1 x−1 dx), or does not exist in a limiting senseR∞(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” propertySample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).0142Chapter 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,ZxNf(x)dx = h[f3/2 + f5/2 + f7/2 + · · · + fN−3/2 + fN−1/2 ]x1(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 withRn=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.}}Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).B2 h2 0(fN − f10 ) + · · ·4B2k h2k(2k−1)(2k−1)(1 − 2−2k+1 )(fN− f1) +···+(2k)!+4.4 Improper Integrals143s=(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);Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).The 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 be144Chapter 4.Integration of FunctionsZZb1/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.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
