c4-3 (779480)
Текст из файла
140Chapter 4.Integration of Functions4.3 Romberg Integration#include <math.h>#define EPS 1.0e-6#define JMAX 20#define JMAXP (JMAX+1)#define K 5Here EPS is the fractional accuracy desired, as determined by the extrapolation error estimate;JMAX limits the total number of steps; K is the number of points used in the extrapolation.float qromb(float (*func)(float), float a, float b)Returns the integral of the function func from a to b.
Integration is performed by Romberg’smethod of order 2K, where, e.g., K=2 is Simpson’s rule.{void polint(float xa[], float ya[], int n, float x, float *y, float *dy);float trapzd(float (*func)(float), float a, float b, int n);void nrerror(char error_text[]);float ss,dss;float s[JMAXP],h[JMAXP+1];These store the successive trapezoidal approxiint j;mations and their relative stepsizes.h[1]=1.0;for (j=1;j<=JMAX;j++) {s[j]=trapzd(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]=0.25*h[j];This is a key step: The factor is 0.25 even though the stepsize is decreased by only0.5.
This makes the extrapolation a polynomial in h2 as allowed by equation (4.2.1),not just a polynomial in h.}nrerror("Too many steps in routine qromb");return 0.0;Never get here.}The routine qromb, along with its required trapzd and polint, is quitepowerful for sufficiently smooth (e.g., analytic) integrands, integrated over intervalsSample 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).We can view Romberg’s method as the natural generalization of the routineqsimp in the last section to integration schemes that are of higher order thanSimpson’s rule. The basic idea is to use the results from k successive refinementsof the extended trapezoidal rule (implemented in trapzd) to remove all terms inthe error series up to but not including O(1/N 2k ). The routine qsimp is the caseof k = 2. This is one example of a very general idea that goes by the name ofRichardson’s deferred approach to the limit: Perform some numerical algorithm forvarious values of a parameter h, and then extrapolate the result to the continuumlimit h = 0.Equation (4.2.4), which subtracts off the leading error term, is a special case ofpolynomial extrapolation.
In the more general Romberg case, we can use Neville’salgorithm (see §3.1) to extrapolate the successive refinements to zero stepsize.Neville’s algorithm can in fact be coded very concisely within a Romberg integrationroutine. For clarity of the program, however, it seems better to do the extrapolationby function call to polint, already given in §3.1.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).0.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
