Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 36
Текст из файла (страница 36)
. . , such that θ N+1 (x) is of degree N and(5.17)When w(x) = 1, a = -1, and b = 1, these polynomials are the Legendre polynomials(see [8, p. 202]). It is also known that the N distinct roots of θ N+1 (x) are real and lie in(a, b). Suppose that an interpolatory quadrature formula (5.2) is based on interpolationat the roots of θN+1 (x).
If f(x) is any polynomial of degree 2N - 1, it can be writtenasf(x) = q(x)θθN+l (x) + r(x),where the quotient q(x) and remainder r(x) polynomials are of degree at most N - 1.Thenwhere the first integral vanishes because of (5.17). For any choice of the nodes {x i } ,the formula of (5.2) integrates a polynomial of degree N exactly, soThe formula applied to f(x) has5.1 BASIC QUADRATURE RULES181Now we use the fact that the xi are roots of θ n(x) to see thatSince any polynomial f(x) of degree 2N - 1 is integrated exactly, this formula has adegree of precision that is at least 2N - 1.There are computationally convenient ways to derive Gaussian quadrature rules,and formulas may be found in specialized books.
Gaussian formulas are valuable because they provide the highest degree of precision for the number of values of f(x).An important fact about Gaussian formulas is that the Ai are all positive. As discussedin connection with the error bounds, this means that we can use formulas of a very highdegree of precision, even for integrands that are not smooth. Gaussian formulas incorporating weight functions are especially valuable tools for dealing with integrands thatare singular and intervals that are infinite. Whether or not there is a weight function,the nodes of a Gaussian formula all lie in the interior of [a, b]. This means that theformula does not use f(a) or f(b). We shall see that this is quite helpful in dealingwith singular integrands.So far we have been considering procedures based on approximating f(x) over thewhole interval [a,b]. Just as with polynomial interpolation, the error depends stronglyon the length of the interval. This suggests that we break up the interval and so approximate the function by a piecewise polynomial function rather than a single polynomial.The simplest approach is to split the interval into pieces specified in advance.
If wepartition [a,b] into a = x1 < x2 < ··· < xn+1 = b, thenand we can apply standard quadrature rules to each of the n integrals. The resultingformula is known as a composite or compound rule. Traditionally the {x i} have beenchosen to be equally spaced in [a,b] and the same formula used on each piece, but thisis not necessary.Example 5.5. Composite Trapezoid Rule. The composite trapezoid rule approximates I =by splitting [a,b] into n pieces of length h = (b - a)/n and applying the trapezoid rule to each piece.
With the definition xi = a + ih, this iswhich simplifies toFigure 5.4 illustrates the composite trapezoid rule.nAn ingenious use of repeated integration by parts establishes the Euler-Maclaurinsum formula. It states that if f(2v)(x) is continuous on [a,b], then182CHAPTER 5NUMERICAL INTEGRATIONFigure 5.4 Composite trapezoid rule.for some x in [a,b]. The coefficients B2k appearing here are known as the Bernoullinumbers. The first few terms of the error expansion areThe basic trapezoid rule applied to an interval of length h has an error that goes to zerolike h3. When the n = 1/h terms are combined, the error of the approximation to theintegral goes to zero like h2.
However, notice that if it should happen that f(1)(b) =f(1)(a), the formula is more accurate than usual. If in addition other derivatives havethe same values at the ends of the interval, the formula is still more accurate. Whenintegrating a periodic function over a multiple of a period, all the derivatives at theends of the interval are equal and this formula is extraordinarily accurate.
In fact, ifthe periodic function is analytic, so that it has derivatives of all orders, TnI fasterthan any power of h! Although rather special, this is extremely important in the contextof Fourier analysis.The error of Tn can be estimated by comparing it to the more accurate result T2nobtained by halving each subinterval. A convenient way to evaluate the formula iswhereIt is important to note that all the evaluations of f made in forming Tn are reused inT2n.There is a way of exploiting the error expansion of the composite trapezoid ruledue to Romberg that is popular for general integrands. The idea is to combine Tn and5.1 BASIC QUADRATURE RULES183T2n in such a way as to obtain a higher order result. According to the error expansion,A little manipulation then shows thatThe formulais of higher order than each of the individual formulas.
As it turns out, this formulais the composite Simpson’s rule. Romberg developed a computationally convenientway of successively combining results so as to increase the order by two with eachcomputation of a composite trapezoid rule. The process is called extrapolation.Romberg integration can be very effective. It adapts the order of the method tothe problem. It does, however, depend on the integrand being smooth throughout theinterval. Also, it evaluates at the ends of the interval, which is sometimes inconvenient.MATHCAD uses Romberg integration for quadrature.
If there is a singularity at an endof an interval or if the process does not converge, the code switches to a variant basedon the midpoint rule that does not evaluate at the ends of the intervals and dividesintervals by 3 rather than 2.nEXERCISES5.1 Use the method of undetermined coefficients to deriveNewton’sand calculate A1 and x1 in the usual manner. Assuming E(f) = cf(d+1)(ξ), find d and c.
What is the corresponding formula and associated error on the general interval [a,b]?5.3 Implement the composite trapezoid rule and apply ittoCalculate A1, A2, A3, A4, d, and c in the usual manner.5.2 Use the method of undetermined coefficients to findthe two-point Gaussian quadrature formula with its associated error. Begin withOf course, you must choose h small enough that samples are taken in each period. Approximate the integral for a number of values of h that tend to 0. According to the theory of Example 5.5, the approximationsTn ought to converge extremely fast.
Is that what youfind?184CHAPTER 5 NUMERICAL INTEGRATION5.2 ADAPTIVE QUADRATUREIn this section a code is developed that approximates I =to an accuracyspecified by the user. This is done by splitting the interval [a,b] into pieces and applying a basic formula to each piece. The interval is split in a way adapted to the behaviorof f(x). A fundamental tool is an error estimate. With the capability of estimating theerror of integrating over a subinterval, we ask if the error is acceptable, and if it is not,the subinterval is split again. As we have seen, for a formula of even modest degree ofprecision, reducing the length of the interval increases substantially the accuracy of theapproximation.
Proceeding in this manner, the formula is applied over long subintervals where f(x) is easy to approximate and over short ones where it is difficult. Codeslike the one developed here are in very wide use, being found, for example, in thecollection of state-of-the-art codes QUADPACK [12], libraries like NAG and IMSL,and computing environments like MATLAB.When the code is supplied absolute and relative error tolerances ABSERR andRELERR, it attempts to calculate a value ANSWER such that|I - ANSWER| < max(ABSERR, RELERR × |I|).The computational form of this uses an error estimateand replaces I by ANSWER:|ERREST| < max(ABSERR, RELERR × |ANSWER|).(5.18)We cannot hope to get a more accurate approximate integral than the correctly roundedtrue value, so it makes no sense to take RELERR < u, the unit roundoff of the computerused.
Indeed, we require RELERR > 10u so that we do not work with error estimatesthat are nothing but roundoff. We also require ABSERR > 0 so as to deal with the raresituation that I = 0.The method employed in the code breaks the interval [a,b] into pieces [α,β] onwhich the basic quadrature rule is sufficiently accurate.
To decide this we must beable to estimate the error of the rule. This is done with a basic principle of numericalanalysis, namely to estimate the error of a result by comparing it to a more accurateresult. Besides the approximationanother approximationsays that whenis formed that is believed to be more accurate. Thenis more accurate thanthe error inis approximately equal to5.2 ADAPTIVE QUADRATURE185To keep down the cost of estimating the error, we use evaluations of f(x) in bothformulas.
As a simple example, we might take Q to be the trapezoid rule and to beSimpson’s rule. The trapezoid rule is based on the values f (α) and f(β). The errorestimate is computed with Simpson’s rule, which needs only the one additional valuef((α+β) /2 Simpson’s rule is considerably more accurate than the trapezoid rule,giving us reason to hope that the estimate will be a good one.The code Adapt uses for Q the three-point Gaussian quadrature formula of Example 5.4 that has degree of precision d1 = 5. A formula of much higher degree ofprecision is used forThe error analysis based on best possible polynomial approximation gives us confidence thatwill be more accurate than Q. It would be possibleto use another Gaussian rule forbut the N-point Gaussian rule for N3 uses Ncompletely different {x i} (except possibly x = 0). To keep the number off evaluationsto a minimum, another approach is taken.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
