c4-2 (779479)
Текст из файла
136Chapter 4.Integration of FunctionsN=1234Figure 4.2.1. Sequential calls to the routine trapzd incorporate the information from previous calls andevaluate the integrand only at those new points necessary to refine the grid. The bottom line shows thetotality of function evaluations after the fourth call. The routine qsimp, by weighting the intermediateresults, transforms the trapezoid rule into Simpson’s rule with essentially no additional overhead.There are also formulas of higher order for this situation, but we will refrain fromgiving them.The semi-open formulas are just the obvious combinations of equations (4.1.11)–(4.1.14) with (4.1.15)–(4.1.18), respectively.
At the closed end of the integration,use the weights from the former equations; at the open end use the weights fromthe latter equations. One example should give the idea, the formula with error termdecreasing as 1/N 3 which is closed on the right and open on the left:Z xN723f2 + f3 + f4 + f5 +f(x)dx = h1212x1 (4.1.20)1351+O· · · + fN−2 + fN−1 + fN1212N3CITED 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 byDover Publications, New York), §25.4. [1]Isaacson, E., and Keller, H.B. 1966, Analysis of Numerical Methods (New York: Wiley), §7.1.4.2 Elementary AlgorithmsOur starting point is equation (4.1.11), the extended trapezoidal rule. There aretwo facts about the trapezoidal rule which make it the starting point for a variety ofalgorithms. One fact is rather obvious, while the second is rather “deep.”The obvious fact is that, for a fixed function f(x) to be integrated between fixedlimits a and b, one can double the number of intervals in the extended trapezoidalrule without losing the benefit of previous work. The coarsest implementation ofthe trapezoidal rule is to average the function at its endpoints a and b.
The firststage of refinement is to add to this average the value of the function at the halfwaypoint. The second stage of refinement is to add the values at the 1/4 and 3/4 points.And so on (see Figure 4.2.1).Without further ado we can write a routine with this kind of logic to it: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).(total after N = 4)4.2 Elementary Algorithms137#define FUNC(x) ((*func)(x))float trapzd(float (*func)(float), float a, float b, int n)This routine computes the nth stage of refinement of an extended trapezoidal 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,...if (n == 1) {return (s=0.5*(b-a)*(FUNC(a)+FUNC(b)));} else {for (it=1,j=1;j<n-1;j++) it <<= 1;tnm=it;del=(b-a)/tnm;This is the spacing of the points to be added.x=a+0.5*del;for (sum=0.0,j=1;j<=it;j++,x+=del) sum += FUNC(x);s=0.5*(s+(b-a)*sum/tnm);This replaces s by its refined value.return s;}}The above routine (trapzd) is a workhorse that can be harnessed in severalways. The simplest and crudest is to integrate a function by the extended trapezoidalrule where you know in advance (we can’t imagine how!) the number of steps youwant.
If you want 2M + 1, you can accomplish this by the fragmentfor(j=1;j<=m+1;j++) s=trapzd(func,a,b,j);with the answer returned as s.Much better, of course, is to refine the trapezoidal rule until some specifieddegree of accuracy has been achieved:#include <math.h>#define EPS 1.0e-5#define JMAX 20float qtrap(float (*func)(float), float a, float b)Returns the integral of the function func from a to b. The parameters EPS can be set to thedesired fractional accuracy and JMAX so that 2 to the power JMAX-1 is the maximum allowednumber of steps. Integration is performed by the trapezoidal rule.{float trapzd(float (*func)(float), float a, float b, int n);void nrerror(char error_text[]);int j;float s,olds;olds = -1.0e30;Any number that is unlikely to be the average of thefor (j=1;j<=JMAX;j++) {function at its endpoints will do here.s=trapzd(func,a,b,j);if (j > 5)Avoid spurious early convergence.if (fabs(s-olds) < EPS*fabs(olds) ||(s == 0.0 && olds == 0.0)) return s;olds=s;}nrerror("Too many steps in routine qtrap");return 0.0;Never get here.}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).(in that sequential order) will improve the accuracy by adding 2n-2 additional interior points.{float x,tnm,sum,del;static float s;int it,j;138Chapter 4.Integration of FunctionsWe come now to the “deep” fact about the extended trapezoidal rule, equation(4.1.11). It is this: The error of the approximation, which begins with a term oforder 1/N 2 , is in fact entirely even when expressed in powers of 1/N .
This followsdirectly from the Euler-Maclaurin Summation Formula,ZxNf(x)dx = hx111f1 + f2 + f3 + · · · + fN−1 + fN22B2 h2 0B2k h2k (2k−1)(2k−1)(fN − f10 ) − · · · −(f−− f1) −···2!(2k)! N(4.2.1)Here B2k is a Bernoulli number, defined by the generating function∞Xtnt=Bnte − 1 n=0n!(4.2.2)with the first few even values (odd values vanish except for B1 = −1/2)B0 = 1B8 = −B2 =13016B10 =B4 = −566130B6 =B12 = −1426912730(4.2.3)Equation (4.2.1) is not a convergent expansion, but rather only an asymptoticexpansion whose error when truncated at any point is always less than twice themagnitude of the first neglected term. The reason that it is not convergent is thatthe Bernoulli numbers become very large, e.g.,B50 =49505720524107964821247752566The key point is that only even powers of h occur in the error series of (4.2.1).This fact is not, in general, shared by the higher-order quadrature rules in §4.1.For example, equation (4.1.12) has an error series beginning with O(1/N 3 ), butcontinuing with all subsequent powers of N : 1/N 4 , 1/N 5 , etc.Suppose we evaluate (4.1.11) with N steps, getting a result SN , and then againwith 2N steps, getting a result S2N .
(This is done by any two consecutive calls ofSample 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).Unsophisticated as it is, routine qtrap is in fact a fairly robust way of doingintegrals of functions that are not very smooth.
Increased sophistication will usuallytranslate into a higher-order method whose efficiency will be greater only forsufficiently smooth integrands. qtrap is the method of choice, e.g., for an integrandwhich is a function of a variable that is linearly interpolated between measured datapoints. Be sure that you do not require too stringent an EPS, however: If qtrap takestoo many steps in trying to achieve your required accuracy, accumulated roundofferrors may start increasing, and the routine may never converge. A value 10−6is just on the edge of trouble for most 32-bit machines; it is achievable when theconvergence is moderately rapid, but not otherwise.4.2 Elementary Algorithms139trapzd.) The leading error term in the second evaluation will be 1/4 the size of theerror in the first evaluation.
Therefore the combinationS=14S2N − SN33(4.2.4)#include <math.h>#define EPS 1.0e-6#define JMAX 20float qsimp(float (*func)(float), float a, float b)Returns the integral of the function func from a to b. The parameters EPS can be set to thedesired fractional accuracy and JMAX so that 2 to the power JMAX-1 is the maximum allowednumber of steps. Integration is performed by Simpson’s rule.{float trapzd(float (*func)(float), float a, float b, int n);void nrerror(char error_text[]);int j;float s,st,ost,os;ost = os = -1.0e30;for (j=1;j<=JMAX;j++) {st=trapzd(func,a,b,j);s=(4.0*st-ost)/3.0;Compare equation (4.2.4), above.if (j > 5)Avoid spurious early convergence.if (fabs(s-os) < EPS*fabs(os) ||(s == 0.0 && os == 0.0)) return s;os=s;ost=st;}nrerror("Too many steps in routine qsimp");return 0.0;Never get here.}The routine qsimp will in general be more efficient than qtrap (i.e., requirefewer function evaluations) when the function to be integrated has a finite 4thderivative (i.e., a continuous 3rd derivative).
The combination of qsimp and itsnecessary workhorse trapzd is a good one for light-duty work.CITED REFERENCES AND FURTHER READING:Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§3.3.Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),§§7.4.1–7.4.2.Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for MathematicalComputations (Englewood Cliffs, NJ: Prentice-Hall), §5.3.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.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
