Đang chuẩn bị liên kết để tải về tài liệu:
Integration of Functions part 4

We can view Romberg’s method as the natural generalization of the routine qsimp in the last section to integration schemes that are of higher order than Simpson’s rule. The basic idea is to use the results from k successive refinements of the extended trapezoidal rule (implemented in trapzd) | 140 Chapter4. Integration of Functions 4.3 Romberg Integration We can view Romberg s method as the natural generalization of the routine qsimp in the last section to integration schemes that are of higher order than Simpson s rule. The basic idea is to use the results from k successive refinements of the extended trapezoidal rule implemented in trapzd to remove all terms in the error series up to but not including O l N2k . The routine qsimp is the case of k 2. This is one example of a very general idea that goes by the name of Richardson s deferred approach to the limit Perform some numerical algorithm for various values of a parameter h and then extrapolate the result to the continuum limit h 0. Equation 4.2.4 which subtracts off the leading error term is a special case of polynomial extrapolation. In the more general Romberg case we can use Neville s algorithm see 3.1 to extrapolate the successive refinements to zero stepsize. Neville s algorithm can in fact be coded very concisely withina Romberg integration routine. For clarity of the program however it seems better to do the extrapolation by function call to polint already given in 3.1. include math.h define EPS 1.0e-6 define JMAX 20 define JMAXP JMAX 1 define K 5 Here 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 s method 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 approxi- int 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 .