Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 37
Текст из файла (страница 37)
In 1964 Kronrod derived for each N-pointGaussian rule a corresponding rule of degree 3N + 1 or 3N + 2 (depending on whetherN is even or odd). The idea was to start with the N Gaussian nodes and add N + 1others chosen to obtain the highest possible degree of precision. Formulas for N < 40are tabulated in [10]. The N = 3 case iswherex1 =(as for three-point Gaussian quadrature)x2 = 0.9604912687080202x3 = 0.4342437493468026A1 = 0.2684880898683334A2 = 0.1046562260264672A3 = 0.4013974147759622A4 = 0.4509165386584744.The basic idea of adaptive quadrature is to approximate the integral over a subinterval (α,β).
If the approximation is not accurate enough, the interval is split into(usually) two parts and the integral approximated over each subinterval. Eventually,accurate approximations are computed on all the subintervals of [a,b] that are addedup to approximate the integral over [a,b] or the cost is deemed too high and the computation terminated. The error terms of the formulasquantify the benefit of splittingan interval; recall from (5.16) thatThe corresponding result for the Kronrod formula of seven points is186CHAPTER 5NUMERICAL INTEGRATIONFigure 5.5 Example where quadrature should be adaptive (for efficiency).Clearly, halving β - α will generally result in much more accurate approximations.This process will resort to short intervals only where f(x) changes rapidly and longones elsewhere.
For example, the function f(x) graphed in Figure 5.5 is likely torequire short intervals near 0 and [2, 3], but not elsewhere. A process like the onedescribed is called adaptive because where f(x) is evaluated in the approximation ofits integral depends on its behavior.It is important to understand that quadrature formulas of the kind taken up inthis chapter sample f at only a finite number of points, and if these points are notrepresentative of the behavior of the curve f(x), the result will not be accurate.
Whatis even worse is that the error estimate comes from comparing the result to that ofanother formula of the same kind, so it is possible that both are inaccurate. Becauseof this any quadrature formula and error estimate of the kind taken up here is doomedto be completely wrong on some problems.
As a concrete example, for the GaussKronrod (N = 3) case, letis a positive number. Yet both formulasClearly, f(x) > 0 on [-1,1], sosee only f(xi ) = 0, hence calculateThe result is terribly wrong and theerror estimate does not detect this. Applying a quadrature code blindly can get youinto trouble!The core of an adaptive quadrature code is a function that approximatesby Q and estimates the error of the approximation byI=to an accuracySuppose we wantTOL := max(ABSERR, RELERR × |ANSWER|).A state-of-the-art code like those of QUADPACK [12] proceeds as follows. At a givenstage of the process the interval [a,b] is partitioned into subintervals a = α1 < β 1 =5.2 ADAPTIVE QUADRATURE187α 2 < β 2 = α3 < ··· < β N = b; there is an estimate Qj available for the integral off(x) over each subinterval [αj ,β j] and an estimate Ej available of the error of Qj.
Byadding up the Qj and Ej, an approximation Q = ANSWER to I is available along withan estimate E = ERREST of its error. If the current approximation to I is not goodenough, the subinterval [αj, β j] with the largest error |Ej| is selected for improvement.It is split into two pieces [α j, (αj + β j) /2], [(αj + β j) /2, β j ], and approximate integralsover these pieces and estimates of their error are computed. The two subintervals andthe associated quantities replace the subinterval [αj , β j ] and its associated quantities.This global adaptive procedure is extremely efficient, but at a cost of keeping trackof all the subintervals and the information associated with them. In the code Adaptthe adaptation is more local and the implementation a little simpler.
The differencein Adapt is that when a subinterval is integrated with enough accuracy, it is no longera candidate for improvement. Thus a queue is kept of subintervals on which the estimated error of the integral over the subinterval is considered to be too large. At atypical stage of the process, the code tests whether the current approximation to I issufficiently accurate to satisfy the user’s error request. If it is not, the next subinterval[αj, β j ] is removed from the queue, split in half, and approximations to the integralover each half, along with error estimates, are computed. They are placed at the end ofthe queue. If an approximation to an integralis estimated to have an errorno more than |(β - α) / (b - a)] × TOL, it is accepted and [α, β] is not examined again.This is more stringent than necessary because adding up approximationsthat satisfy this condition will result in an error such thatThe global adaptation used in the collection QUADPACK [12] subdivides untilTOL.
The more local adaptation of Adapt subdivides untilTOL. Becauseand perhaps much larger, the local processis always at least as big asworks harder but should provide answers that are more accurate than the specifiedtolerance TOL.Let us now formalize the algorithm. To start the process, form the approximationQ to the integral over [a,b] and its estimated error ERREST. If |ERREST| < TOL, wetake ANSWER = Q as the result and return. If ERREST does not pass the error test,a, b, Q, and E = ERREST are placed in the queue and the following loop is entered:remove α, β, Q, E from top of queuecompute QLwith estimated error ELcompute QRwith estimated error ERANSWER := ANSWER + (( QL + QR) - Q)ERREST := ERREST + ((EL + ER) - E)if EL is too big, add α, (α + β)/2, QL, EL to end of queueif ER is too big, add (α + β)/2, β, QR, ER to end of queue.188CHAPTER 5NUMERICAL, INTEGRATIONThis procedure is repeated until one of the following events occurs:1.2.3.4.The queue becomes empty.|ERREST| < TOL.The queue becomes larger than the space allocated.Too many f evaluations are made.The first two cases represent a success.
The last two represent a failure in the sensethat the code has failed to achieve the desired accuracy in the work allotted. It maybe that the estimated error, although larger than specified by means of the tolerances,is acceptable and the answer reported will suffice. Even when it is not, an inaccuratevalue of the integral may be useful as a indication of the size of the integral whenselecting appropriate tolerances.Notice the way quantities have been grouped in the computation of ANSWER andERREST.
The quantity Q and the sum QL + QR both approximate the integral off over[α, β]. They normally agree in a few leading digits, so their difference involves cancellation and is computed without arithmetic error. Because the difference is normally asmall correction to ANSWER, it is possible to make a great many corrections withoutaccumulation of arithmetic error. If the quantities are not grouped, the correction of a“small” ANSWER may be inaccurate.5.3CODES FOR ADAPTIVE QUADRATUREThe code presented here uses the three-point Gaussian rule to estimate integrals alongwith the seven-point Kronrod rule to estimate errors.
In FORTRAN the routine Adapthas a typical callCALL ADAPT (F, A, B, ABSERR, RELERR, ANSWER, ERREST, FLAG)and in C++,flag = Adapt(f, a, b, abserr, relerr, answer, errest);and it isflag = Adapt(f, a, b, abserr, relerr, &answer, &errest);in the C version.The first argument, F or f, is the name of the function that provides values ofthe integral.
In FORTRAN F must be declared in an EXTERNAL statement. Thenext four variables are input parameters: A and B are the (finite) end points of theintegration interval, ABSERR and RELERR are the required absolute and relative error tolerances. The remaining variables are output quantities: ANSWER contains theapproximation to the integral and ERREST an estimate for the error of the approximation. The value of FLAG is 0 for a normal return with (5.18) satisfied. FLAG > 05.3 CODES FOR ADAPTIVE QUADRATURE189signals that there was not enough storage for the queue, or that too many function evaluations were needed.
Illegal input (ABSERR < 0 or RELERR < 10u) is indicated byFLAG = -1. In C and C++ the value of flag is the return value of the function Adapt.Example 5.6. Let us try the code on a problem with an analytical integral, for exfor which I = e - 1 = 1.71828182845904. Its output for a requestedample,= 10 - 5 and RELERR = 10 - 8 follows.accuracy of ABSERR0FLAG =Approximate value of the integral = 1.718281004372522Error in ANSWER is approximately8.240865232136876E-0077 evaluations of F were required.Routine Adapt evaluates f at least seven times, so the code found this task to be veryeasy.
It is seen that ERREST approximates the true error of 0.82408652 × 10-6 exntremely well.againExample 5.7. A more interesting problem is to estimate-5-8with ABSERR = 10 , RELERR = 10 . Although continuous at x = 0, the integrandhas infinite derivatives there, so it does not behave like a polynomial near this point.An adaptive code deals with this by using short intervals near the point.
Applicationof the Adapt code results in0FLAG =6.718072986891337E-001Approximate value of the integral =-6.448552584540975E-005Error in ANSWER is approximately119 evaluations of F were required.The techniques of the next section can be used to produce an accurate value of I =0.671800032402 for the integral, which tells us that the actual error is about -0.73 x10-5. In comparison to Example 5.6, this problem required many more evaluations off, and ERREST is not nearly as accurate. It should be appreciated that all we need is arough approximation to the error, so this one is perfectly acceptable. It is instructive tosee how Adapt samples the integrand in this example.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
