Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 38
Текст из файла (страница 38)
In Figure 5.6 ordered pairs (x,y)are plotted (as are pairs (x, 0) along the x-axis); x is a sampling point used by Adapt andy the corresponding value of the integrand. Notice how the points congregate aroundnthe singularity at zero as the theory predicts they should.EXERCISESUnless otherwise indicated, use as tolerances ABSERR = 10-8 and RELERR = 10-6 for the computing exercises.5.4 To test out Adapt try the following integrands.(a) 4/(1 + x2)(b) x 1/10(c) 1/[(3x - 2)2](d) 1 + sin2 38πx190CHAPTER 5NUMERICAL INTEGRATIONFigure 5.6 Sampling of the integrand of Example 5.8 by Adapt.FLAG = 0, ERREST = 0, yet ANSWER is totallywrong. Explain.Use ABSERR = 10-12, RELERR = 10-6, and A = 0,B = 1 for all parts.
Which problems require the mostfunction evaluations? Why? What are the exact answers? What are the actual errors? How good an estimate is ERREST? On which problems is it better thanothers? Why?5.9 A sphere of radius R floats half submerged in a liquid. If it is pushed down until the diametral plane isa distance p (0 < p < R) below the surface of theliquid and is then released, the period of the resultingvibration is5.5 The integrand in Exercise 5.4d has period l/38.where k2 = p2/(6R2 - p2) and g = 32.174 ft/sec2. ForR= 1 find T when p = 0.50, 0.75, 1.00.Rewrite this asand useAdapt on this with the same tolerances as in Exer- 5.10 A population is governed by a seasonally varying abilcise 5.4.
Is this approach faster? More accurate?ity of the environment to support it. A simple modelof this is the differential equation5.6 If a package is not available for evaluating the zeroth order Bessel function J0(x), then an alternativeis based on the formulaUse the formula to evaluate J0(x) for x = 1.0, 2.0, 3.0.Compare with tabulated results (e.g., [7]).5.7 The function y(x) =is called Dawson’sintegral. Tabulate this function for x = 0.0, 0.1, 0.2,0.3, 0.4, 0.5.
To avoid unnecessary function evaluations, split the integral into pieces.5.8 Derive a step function f(x) (a function that is constant on subintervals of [a,b]) for which Adapt returnswhere t measures time in months, P(t) is the population at time t, and the remaining parameters are knownconstants. This equation will be solved numerically inthe next chapter (Exercise 6.19); here we note that theproblem can be transformed intowhere F(t) = exp[kM( t - (6r) /πsin(πt /6))]. Assume that k = 0.001, M = 1000, r = 0.3, P(0) =5.4 SPECIAL DEVICES FOR INTEGRATION250 and use Adapt efficiently to table P(t) for t =0, 3, 6, 9,. .
. , 36.5.11 Examine the effect of noise on Adapt as follows. Forthe input function F(x) useT 1 := f(x) × 10n191y´(0). Observing that y´´ = exp(y) - 1 can be writtenin the formwe can integrate to obtainT 2 := T1 + f(x)F(x) := T2- T1where f(x) is the original integrand (use some of theexamples in Exercise 5.4). With ABSERR = RELERR(u the unit roundoff) tryand then nWhat is the behavior? Doesthe algorithm appear to be stable?5.12 To solve the nonlinear two-point boundary value problemusing standard initial value problem codes (e.g., Rke),it is necessary to find the missing initial conditionSince y(0) = 0, this says y´(0) =y´(x) (by separation of variables) yieldsSolving forwhich, when evaluated at x = 1, becomesUse Adapt and Zero to find c and then y´(0).5.4 SPECIAL DEVICES FOR INTEGRATIONQuadrature formulas approximate integrals using a finite number of samples.
If thesamples are not representative, the result may be inaccurate despite an estimated errorthat is acceptable. Put differently, the approximate integral and its error estimate arebased on an assumption that the integrand changes smoothly between samples. Adaptive codes generally do very well at recognizing the behavior of integrands, but f(x)with sharp peaks or many oscillations in the interval present special difficulties.
Sometimes it is necessary to assist a code by breaking up the interval of integration yourselfso as to ensure that the code will take samples in critical areas. A contrived examplewill help make the point.The family of integralscan be evaluated readily with the recurrenceWhen n is large, the integrand has a sharp peak at the midpoint of the interval. Withtolerances of ABSERR = 10-5 and RELERR = 10-11, for n = 200 Adapt returnedANSWER = 0.1252567600019366 and an estimated error of about -3.654 × 10-6 ata cost of 203 evaluations of the integrand. The true error is about -3.654 × 10-6. Thecode had no difficulty producing an accurate answer and an excellent error estimatebecause it samples at the midpoint and “sees” the peak.
However, if the interval is split192CHAPTER 5NUMERICAL INTEGRATIONinto [0,2.6] and [2.6,π], the results are quite different. Integrating over the two intervals and adding the results provides an approximate integral of 4.111202459491848 ×10-7 and an estimated error of about -1.896 × 10-7 at a cost of 14 evaluations of theintegrand. The code has been completely fooled! This happened because the initialsamples did not reveal the presence of the peak. The code took the minimum numberof samples from each interval, namely seven, showing that it believed the problem tobe very easy.
When this happens it is prudent to consider whether you agree that theproblem is very easy and if not, to break up the interval into pieces that will cause thecode to “see” the behavior of the integrand. Of course, one must be careful how to dothis breaking into pieces, as [0,2.6] and [2.6,π] won’t do for this problem.Adapt is based on approximation by polynomials over finite intervals.
As a consequence it may have to resort to a great many subintervals to integrate functions thatdo not behave like a polynomial near a critical point or to integrate functions over aninfinite interval. Gaussian quadrature formulas with a suitable weight function are agood way to handle such problems. Specific formulas can be found in sources like[13].
Substantial collections of quadrature codes such as those of QUADPACK, NAG,and IMSL contain specialized routines for a variety of difficulties. In this section wediscuss some techniques for preparing problems to make them more suitable for ageneral-purpose code like Adapt. With a little mathematical effort, a problem that thecode cannot solve might be put into a form that it can solve, and problems that it cansolve are solved more efficiently.OSCILLATORY INTEGRANDSIf the integral f(x) is periodic with period p, that is, f(x + p) = f(x) for all x, andb - a is some integer multiple of p, the integration should be performed over just oneperiod (or sometimes less) and symmetry used for the remainder.
For example, for theintegral in Exercise 5.3,If you have worked Exercise 5.3 you should have found that the composite trapezoidrule is very efficient; consequently, if you have many functions to be integrated over aperiod, it is advisable to use a special-purpose code based on the composite trapezoidalrule.Oscillatory nonperiodic functions are more difficult.
Generally the integral shouldbe split into subintervals so that there are few oscillations on each subinterval. Adaptmay do this automatically, but the computation can be made more reliable and possiblyless expensive by a little analysis. As an example, let us estimateThis integral could be rewritten many ways, one of which is5.4 SPECIAL DEVICES FOR INTEGRATION193Proceeding in this manner, Adapt is called 20 times, but on each interval the integrand has constant sign and varies smoothly between samples.
This is a reasonableway to deal with a single integral of this kind, but in contexts like Fourier analysis,where many functions having the form of a periodic function times a nonperiodicfunction need to be integrated, it is better to use special-purpose formulas such asproduct quadrature. This particular example is treated in Case Study 3 where the useof the general-purpose code Adapt is contrasted with Filon’s method for finite Fourierintegrals.INFINITE INTEGRATION INTERVALAdapt cannot be applied directly tobecause it deals only with finite intervals. One way to apply the code is to use the definitionThe idea is to determine an analytical bound for the tailWith it an endpoint b can be selected large enough thatequalsto the accuracyrequired. It does not matter if b is rather larger than necessary, so a crude bound forthe tail will suffice.Another way to get to a finite interval is to change the variable of integration.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
