Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 39
Текст из файла (страница 39)
Forexample, to estimatethe new variable s = 1/x yields the equivalentintegralon a finite interval. Generally this trades one difficulty(an infinite interval) for another (a singularity at an end point). For this particularso the integrand is continuous at s = 0 andexample,Adapt can be applied directly. If the original problem werethens = x -1 would not help because it produces the same infinite interval.
The choices = e-x leads tofor which the integrand approaches 0 asOn the integration intervalthe transformation s = arctan x is often useful.SINGULAR INTEGRANDSBy a singularity we mean a point where the function or a low order derivative is discontinuous. There are two cases: (1) finite jumps and (2) infinite values. In either case,if there is a known singularity at a point c in (a,b), a basic tactic is to split the intervalinto (a,c) and (c,b).
In this way we can arrange that we approximate integrals withintegrands that are finite except possibly at one end.If f(x) is finite at the singularity (as for step functions), Adapt can be used on eachpiece and the results added to produce an approximation toThis is clearenough, but it must be kept in mind that a function that is smooth to the eye may not besmooth to the code. Surprisingly often someone integrates a spline with a quadratureroutine. Of course, such an integration should be done analytically, but it may beconvenient to resort to an adaptive quadrature routine.
Just remember that splines havediscontinuities in low order derivatives. If a routine like Adapt is applied to each pieceof a piecewise polynomial function, it will have no difficulty. Indeed, it will be exact194CHAPTER 5NUMERICAL INTEGRATIONif the degree of precision of the basic formula is sufficiently high. However, if noattention is paid to the lack of smoothness at the knots, the code will have to deal withit automatically. This can represent quite a substantial, and completely unnecessary,expense as the code locates the knots and resorts to short intervals to deal with them.Infinite discontinuities require some study. First we reduce the problem to comwith a singularity at one end, let us say a for definiteness.puting an integralThen we have to sort out the behavior of the integrand at a to convince ourselves thatthe integral even exists.
We have already seen problems for which the singularity islogarithmic, f(x) ~ c ln (x) asand (in Example 5.7) algebraic, f(x) ~ x1/7 asIn the case of an algebraic singularity, f(x) ~ c (x - a)γ asit is necessarythat γ > -1 for the integral to exist. The behavior at the singular point tells us whatkind of weight function to introduce if we wish to deal with the difficulty by means ofa special formula.One way to deal with singularities using a general-purpose code like Adapt is tointroduce a new variable with the aim of removing, or at least weakening, the singularity.
This can be done quite generally for singularities of the form f(x) ~ c ( x - a) γ .Let us try a new variable t with x = a + t β. ThenThe new integrand G(t) ~It will be nonsingular if β(γ + 1) - 1 > 0. Bychoosing β = l/(γ + 1) (recall that γ > -l), the function G(t) ~ cβ asand wecan apply our favorite code to the new problem. If the code evaluates at the ends ofthe interval of integration, we have to code the subprogram for G(t) carefully so thatthe proper limit value is returned when the argument is t = 0. For such codes wemight prefer to take, say, β = 2/(γ + 1) so that G(t) ~ cβt and the limit value is simplyG(0) = 0.To illustrate the technique, suppose we want to computeOften series expansions are an easy way to understand the behavior of an integrand.Here the difficulty is at x = 0.
Using the series for exp(x),Thus f(x) ~ x-1/3 and γ = -1/3. The integral exists because γ > -1. If we takeβ = 2/(γ + 1) = 3, the integral becomeswhich presents no difficulty at all.5.4 SPECIAL DEVICES FOR INTEGRATIONExample 5.8.195As a more substantial example, suppose we wish to approximateUsing the series expansionwe haveThe integrand f(x) ~ x-1/4 asso the integral exists and a suitable change ofvariables would be x = t4/3.
Since Adapt does not evaluate at endpoints, it can beapplied in a straightforward way to the transformed integralApplying Adapt directly to the original problem with ABSERR = 10-10 and RELERR = 10-8, results in ANSWER = 1.913146656196971 and an estimated error ofabout 1.485 × 10-8 at a cost of 1841 evaluations of f. On the other hand, applying the code to the problem after the change of variables results in ANSWER =1.913146663755082 and an estimated error of about 1.003 × 10-8 at a cost of 147evaluations off.nA technique called subtracting out the singularity is a valuable alternative to achange of variables.
The integrand is split into two terms, one containing the singularity that is integrated analytically and another that is smoother and is approximated numerically. The technique is an alternative to using special formulas based on a weightfunction suggested earlier for the integrals arising in plane potential computations:Before it was assumed implicitly that f(x) is well behaved at x = 0. Thenf(x) ln(x) ~ f(0) ln(x)196CHAPTER 5NUMERICAL INTEGRATIONFigure 5.7 Integrands for Examples 5.8 and 5.9 near x = 0.and we can writeThe first integral is done analytically and the second numerically.
The integral tobe computed numerically is easier than the original one. This is seen by expandingf(x) = f(0) + f´(0)x + ··· and observing that the integrand behaves like f´(0)x ln(x)asIt is easier to apply this technique than to change variables, but a good change ofvariables will be more effective because it deals with all the singular behavior ratherthan just the leading term(s) that are subtracted out.Example 5.9. Let us compare the techniques for the problem in Example 5.8. Insubtracting the singularity, we writeWith the same tolerances as before, the output from Adapt isANSWER = 1.913146679435873, ERREST = -5.476 × 10-9at a cost of 287 calls to f.
See Figure 5.7 for a plot near x = 0 of the original integrandand those from Examples 5.8 and 5.9.nWe end this section with a more realistic problem.5.4 SPECIAL DEVICES FOR INTEGRATION197Example 5.10. A conducting ellipsoidal column projecting above a flat conductingplane has been used as a model of a lightning rod [2]. When the ellipsoid is given bythe equationthe potential function is known to beThe constant A depends only on the shape of the rod and is given byThe quantity λ is a function of the place where we want to evaluate the potential. It isthe unique nonnegative root of the equationSuppose that for a tall, slender rod described by a = 1, b = 2, and c = 100, we seekthe potential V at x = y = z = 50.
As we saw in the preceding chapter, the root λ of φis approximately 5928.359. To compute(5.19)we note that the integrand tends to zero like u-5/2 asso the integral is welldefined and the change of variables u = 1/t is satisfactory. This substitution producesThis is acceptable mathematically but is in poor computational form.
Here it is easy torearrange to obtainUsing Adapt with RELERR = 10-5 and ABSERR = 10-l4 producesANSWER = 0.5705983 × 10-6, and ERREST = -0.52 × 10-11,which requires 105 f evaluations. The integrand is approximatelynear t = 0, whichsuggests that a better change of variables could have been made (e.g., u = 1/w 2 ). Thenthe integral is198CHAPTER 5NUMERICAL INTEGRATIONThis is a somewhat easier integral to evaluate. Routine Adapt producesANSWER = 0.5705883 × 10-6, and ERREST = 0.46 × 10-11,which requires 35 f evaluations. Note that the two results agree to the requested fivefigures.These transformations cannot be applied directly to computing A because the interval remains infinite.
However, since (5.19) has to be computed anyway, we couldsplitintoThe first integral is computed to be 0.7218051 × 10-5 with an errorestimate of 0.31 × 10-10 and requires 357 f evaluations, while the sum using thesecond value from above is 0.7788639 × 10-5. This yields A = 0.1283921 × 106 andfinally V = -46.3370.nEXERCISESUnless otherwise indicated, use as tolerances ABSERR = 10-8 and RELERR = 10-6 for the computing exercises.5.13 Evaluateusing Adapt(a) on the problem as is,Try Adapt with RELERR = 10-2 and ABSERR =½ × 10-3, ½ × 10-6, and ½ × 10-9.(c) Compare the answers in (a) and (b). Discuss efficiency and accuracy.5.16 The integral(b) using one of the series techniques, and(c) using a change of variable.Compare the results.5.14 Repeat Exercise 5.13 with5.15 The integralhas a rather nasty singularity at x = 0.
Analyze thenature of the singularity and argue informally that theintegral exists.(a) Use Adapt for this problem as it is. How is thiseven possible?can be difficult to estimate because of the large number of oscillations as(although the amplitudesare approaching zero, too).(b) Treat the singularity using a technique from Section 5.4, then use Adapt.(c) Compare the answers in (a) and (b).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
