Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 41
Текст из файла (страница 41)
There is the additional and very serious complicationof splitting up the region of integration and transforming properly each piece to astandard region. The whole area of integration of functions of more than one variableis still the subject of research.EXERCISES5.24 Given the triangular domain T with vertices (0,0),(1,0), and (01), we would like to approximateover T.(a) Derive a quadrature approximation of the formA 1 f (0, 0) + A2f(0, 1) + A3f(1, 0)where the coefficients are chosen to make the approximation exact for f = 1, x, and y.(b) Derive a corresponding composite quadrature formula for the subdivision obtained by cutting T intofour subtriangles by connecting the midpoints of theedges of T.5.7 CASE STUDY 5This case study has two parts, one devoted to a problem involving singular integrandsand the other to problems with integrands that have sharp peaks.
As usual, our aim isto understand the difficulties and what might be done to make it possible to solve theproblems effectively with standard codes. The first part is an example of a problemrequiring more than one code for its solution. The second part develops a technique ofclassical applied mathematics and uses it for insight.Reference [5] discusses the motion of a particle of mass m at position q in a potential field V(q). In suitable circumstances the motion is a libration (or oscillation or204CHAPTER 5NUMERICAL INTEGRATIONvibration).
For an energy level E, it is found that if the equation E - V(q) = 0 hassimple roots q1 < q2, the period T of the motion is given byThe integrand is singular at both end points, so the first thing we must do is sort out itsbehavior there. Near a simple root qi , a Taylor series expansion tells us that ashence the singularity is integrable. The argument shows why the theory requires thatthe roots be simple. If qi were a multiple root, say double, the integrand would behavelikeand the integral would not exist. Percival and Richards [5] write, “Gaussian numericalintegration is a very effective method of obtaining the periods of libration in practice.”Gaussian quadrature is generally quite effective, but a crucial point in the present situation is that it can be applied directly because it does not evaluate the integrand at theends of the interval.
The formulas of Adapt have the same property.As a simple example let us take V(q) = (q + l)(q - 0.8)7 and E = -4. As canbe deduced easily from the derivative V´(q) = (8q + 6.2)(q - 0.8)6 or by inspectionof Figure 5.8, V(q) has a minimum at q = -0.775, and for V(-0.775) < E < 0, theequation E - V(q) = 0 has two simple roots, one in [-1, -0.775] and the other in[-0.775,0.8]. To evaluate the integral, we must first compute these roots numerically.The roots should be located very accurately because they define the integral that wewish to approximate. Moreover, if they are not determined accurately, a quadraturecode might evaluate the integrand at an argument q for which the integrand is notdefined because of taking the square root of a negative number or dividing by zero.A problem like this is very easy for Zero, so we take ABSERR = 10-l4 and RELERR = 10-14.
Even with these stringent tolerances the code requires only 21 functionevaluations to find that-0.9041816 with a residual of about 1.3 × 10-l5 and-0.5797068 with a residual of about 1.9 × 10-15. Because E provides a naturalmeasure of scale for this equation, the residuals and a backward error analysis make itclear that the roots are very accurate, although of course we already know that becauseZero reported that it had obtained the specified accuracy.Using the computed values for the end points, it is possible to approximate theintegral by a simple call to Adapt.
With ABSERR = 10-6 and RELERR = 10-6 thisresults in an answer of about 0.444687 with an estimated error of about 9.0 × 10-7 at acost of 3171 function evaluations. This is relatively expensive for Adapt with its (arbitrary) limit of 3577 evaluations. If we wanted to do a number of integrals of this kind,some preparation of the problem would be worthwhile. We have already worked outthe dominant behavior of the integrand at the end points and the value of the derivativeappearing in the expression is readily available, so the method of subtracting out thesingularity would be easy to apply.
However, because the singularity is algebraic, it5.7 CASE STUDY205Figure 5.8 The potential V(q) = (q+ l)(q - 0.8)7.is easy to change variables and this technique is generally more effective. We need tobreak the interval of integration so that there is at most one singularity in each pieceand when present, it is an end point. A natural choice is [q1, -0.775] and [-0.775, q2 ].In the notation used earlier in the chapter, γ = -1/2, so we might choose β = 2, thatis, introduce the variable t2 = q - q1 in the portion to the left and t2 = q2 - q in theportion to the right.
Adapt is used twice with the same tolerances as before and theintegrals are added to obtain an answer of about 0.444688. The estimated errors canalso be added since they are estimates of the error rather than the magnitude of theerror, and this gives an estimated error of about 8.2 × 10-7. The value for the integralis consistent with that of the direct approach, but it costs a total of only 56 functionevaluations.Let us turn now to a different kind of difficulty, integrands with sharp peaks.
Manyimportant problems of applied mathematics are solved by transform methods that leadto the task of evaluating integrals. A method developed by Laplace illustrates the use ofasymptotic methods for this purpose. A family of integrals depending on a parameters >> 1 of the formis considered. If the function f(x) has a unique maximum at x0 with a < x0 <theintegrand has a sharp peak at x0, and the greater s is, the sharper the peak. The idea isfirst to approximate f(x) near x0 byand then to approximate the integral by integrating this function. Because the integranddecays so rapidly for large s, the approximation to the integral is scarcely affected byextending the interval of integration toThis is done in order to integrate the206CHAPTER 5NUMERICAL INTEGRATIONapproximating function analytically.
It amounts to approximating the integrand by thefamiliar bell-shaped curve of a normal distribution with mean x0 and standard deviationThe result of these approximations is Laplace’s formula,A classic application of this formula is to the Stirling approximation of the gammafunction seen in Example 1.10. Some manipulation of the integral definition of thegamma function puts it into the form required:(Here t = xs.) Laplace’s formula with f(x) = lnx - x and x0 = 1 then givesLaplace’s formula and the class of integrals it approximates illustrate how the approximations of classical applied mathematics differ from those of numerical analysis.A general-purpose code like Adapt that accepts any smooth integrand is likely to failwhen presented an integrand with a sharp peak because it is not likely to “see” thepeak.
By this we mean that the code is not likely to take samples from the small subinterval where the peak occurs, so it, in effect, approximates the integrand by a functionthat does not have a peak. The approach of applied mathematics is much less generalbecause it requires that the location of the peak be supplied along with informationabout its width.
An advantage of the approach is that it provides important qualitativeinformation about how the integral depends on a parameter. On the other hand, theaccuracy of the formula depends on the parameter, and for a given value of s it mightnot be enough. As is often the case, when used to obtain a numerical approximationto an integral, the approaches are complementary with the one working better as theparameter increases and the other working worse, and vice versa. Let us take advantage of the insight provided by Laplace’s method to compute an accurate value foran integral with a sharp peak using a general-purpose code.
Clearly we ought first tolocate the peak, then get some information about its width, and finally break up the interval appropriately. D. Amos [l] does exactly this to evaluate integral representationsof statistical distributions. In each case an evaluation of the function amounts to thenumerical integration of a bell-shaped integrand with a single maximum at x0. FirstNewton’s method is used to locate x0. A scale σ is then determined using the Laplacemethod. Finally quadratures over intervals of length σ to the left and right of x0 aresummed until a limit of integration is reached or the truncation error is small enough.The basic idea is simple, but the generality of the distributions allowed and the refinements needed for fast evaluation of the functions make the details too complex todescribe here.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
