Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 40
Текст из файла (страница 40)
Discuss efficiency and accuracy.(a) Use Adapt for this problem with RELERR =10-l2 and the three absolute tolerances ABSERR = 5.17 The exponential integral10 -3, 10 -6, and 10-9.(b) The change of variable t = 1/x produces theequivalent integralApproximate thisbywhere b is chosen sufficiently large.How large should b be in order to guarantee thatt > 0, arises in the study of radiative transfer and transport theory [9]. Some manipulation shows that5.4 SPECIAL DEVICES FOR INTEGRATION1995.20 The Gibb’s effect describes the behavior of a Fourierseries near a discontinuity of the function approximated.
The magnitude of the jump can be related tothe integralThe expression in braces is known to have the value0.5772157, the Euler-Mascheroni constant. Thesecond term integrates analytically to - lnt. Hence,Evaluate E1 (t) for t = 1.0, 2.0, 3.0. Does there appearto be any difficulty caused by the behavior of the integrand at x = 0?5.18 In studying the conduction properties of a V-groovedheat pipe, a heat flux constant C(θ) satisfieswhich is about 0.281. Routine Adapt cannot be applied directly because it requires a finite interval. Oneway of dealing with this is to drop the tail of the integral, that is, approximate the integral byThis does not work well in this instance because theintegrand decays slowly asTo see this, workout an analytical bound for the error made in droppingthe tail, that is, bound( θ in radians).
Compute values of C for θ = π/ 6, π/4,π /3, 5π/12, 17π/36. It can be shown mathematicallythat C(θ) is strictly decreasing and 0 < C(θ) < 1. Thedenominator in the integrand vanishes at π/2. Use series expansions to sort out the behavior of the integrand at this point, and if it causes difficulty, fix it.Integration by parts leads to5.19 Debye’s specific heat law gives the molar specific heatof a solid, Cv, as a function of its temperature, T:where R = 1.986 calories/mole is the molar gas constant, UD = θD/T, and θ D is the Debye temperaturethat is characteristic of each substance. For diamondθD = 1900K, evaluate Cv(T) at the temperatures indicated in the accompanying table. Compare with theexperimental values given.
Does the behavior of theintegrand at x = 0 appear to cause problems?The integral arising here is easier to approximate because the integrand decays faster. Integration by partscan be repeated to get integrals that are even easierto approximate. Approximate the original integral accurately enough to verify that its magnitude is about0.281. Do this by applying Adapt to the integral thatarises after integrating by parts a few times. To apply the code you will have to drop the tail. There aretwo sources of error in the computation. One is theerror made by Adapt, which you can control by thetolerances that you specify. The other error is due todropping the tail.
You can control this error by working out an analytical bound for the error and choosinga value b large enough to guarantee the accuracy thatyou need.200CHAPTER 5NUMERICAL INTEGRATION5.5 INTEGRATION OF TABULAR DATAThe problem addressed here is the approximation ofgiven only (xn,yn) for1 < n < N, where yn = f(xn). Adaptive quadrature routines cannot be used becausethey automatically choose the points where f is to be evaluated and it is unlikely thatsuch points will correspond to the data points {xn}. The basic approach is the same asin Section 5.1: approximate f(x) by a piecewise polynomial function F(x), which isthen integrated exactly.Since the complete cubic interpolating spline worked so well as an approximatingfunction in Chapter 3, it is a natural choice for F(x). For the sake of simplicity, let usassume that a = x1 and b = xN.
Then, using the notation of Chapter 3 for the spline,Substituting the expressions for an, bn, and dn in terms of the data and cn leads to(5.20)An algorithm based on this technique first uses SPCOEF/Spline_coeff from Chapter 3to getand then evaluates (5.20). In terms of h = maxn(xn - xn-1), the complete cubic spline S(x) provides an approximation to a smooth function f that is accurate to O(h4). Accordingly, if the sample points xi are sufficiently dense in [x1 ,xN] andthe function is sufficiently smooth, we can expectto be an O(h4) estimateofThe cubic spline is familiar and it is easy enough to manipulate once the linearsystem forhas been solved, but there are other interesting possibilities. Thelinear system arises because S(x) is required to have two continuous derivatives on[x 1 ,x N ]. This smoothness is unnecessary to the approximation ofTheshape-preserving spline of Section 3.5.2 is less smooth but is attractive here for severalreasons. The reasons and the following formula are left as an exercise:(5.21)A widely used scheme is based on local quadratic interpolation.
To approximatef(n) by a quadratic over [x n ,x n+1 ], we must interpolate it at three points. One possi-5.5 INTEGRATION OF TABULAR DATA201bility is to interpolate at {x n - 1 ,x n ,x n + l }. Another is to interpolate at {x n ,x n + 1 ,x n + 2 } .There is no obvious reason to prefer one to the other and the computation is inexpensive, so a reasonable way to proceed is to compute both and average the two results.This provides a symmetrical formula that smooths out mild measurement errors in thetabulated values off. Of course, at the ends n = 1 and n = N - 1, only one quadraticis used.
The formula for a typical (interior) interval [x n ,x n + l ] is(5.22)Reference [4] contains further discussion and a FORTRAN code AVINT.There is no obvious way to obtain a good estimate of the error of these quadraturerules, much less to control it. The difficulty is inherent to functions defined solely bytables.EXERCISES5.21 In performing an arginine tolerance test, a doctor measures glucose, insulin, glucagon, and growth hormonelevels in the blood over a l-hour time period at 10minute intervals to obtain the following data:time0102030405060glucose102114122132115107100insulin11263647392715timeglucagon01020304050601881300230026001800840460growthhormone1.701.701.202.507.258.108.00The doctor is interested in the integrated effect of eachresponse.
For example, if the glucose curve is represented by g(t), thenis desired. Computeone of the integrals by(a) the method (5.20) based on splines,(b) the method (5.21) based on splines, and(c) the method based on averaging quadratics (5.22).Compare the results.5.22 Supply the details of the derivation of (5.22). Workout the special forms for the end points.5.23 Derive the formula (5.21) for the integration of tabulardata using the shape-preserving spline. When mightyou prefer this to formula (5.20) based on the complete cubic spline? Consider the cost of computingthe spline and how well the spline might fit the data.202CHAPTER 5NUMERICAL INTEGRATION5.6 INTEGRATION IN TWO VARIABLESDefinite integrals in two or more variables are generally much more difficult to approximate than those in one variable because the geometry of the region causes trouble.
Inthis section we make a few observations about the common case of two variables,especially as it relates to finite elements.Integration over a rectangle,can be handled easily with the formulas for one variable by treating I as an iteratedintegral. Thus we first approximatewith one quadrature rule using N1 points {x i } and thenusing another rule of N2 points {y j}. This approach generalizes toIt is even possible to use an adaptive code for the integrals in one variable, but thematter is a little complicated.
In Fritsch, Kahaner, and Lyness [6] it is explained howto go about this; for pitfalls to be avoided, see [11, Section 9].Degree of precision now refers to polynomials in two variables, so a formula of,for example, degree 2 must integrate all polynomials of the formexactly on the region of interest. Just as in the case of one variable, we can derivequadrature formulas by interpolating f(x,y) and integrating the interpolant.
This isquite practical for integration over a square or a triangle and is often done. The schemefor rectangles based on iterated integrals can be quite efficient when the formulas forintegration in one variable are Gaussian. They are not necessarily the best that can bedone. As in the one-dimensional case, formulas can be derived that use the smallestnumber of evaluations of f(x,y) to obtain a given degree of precision. Nevertheless,the most efficient formulas may not be the most attractive in practice.
The approachbased on interpolation can be very convenient when the interpolation is done at pointsinteresting for other reasons, as in finite elements. The iterated integral approach canbe very convenient because of its simplicity and generality.In one dimension the transformation of any finite interval [a,b] to a standard onesuch as [-1,l] is trivial. In two dimensions the matter is far more important and5.7 CASE STUDY 5203difficult.
Now an integral over a general region R,must be broken up into pieces that can be transformed into a standard square or triangle. Discretizing a region R in this way is an important part of any finite element code.If the region is decomposed into triangles (with straight sides), the easy transformationwas stated in Chapter 3. An integral over a general triangle T becomes an integral overthe standard, reference triangle T*,Here |D*| is the determinant of the Jacobian of the transformation. It relates the infinitesimal area dydx in the one set of variables to that in the other set. For the affinetransformation from one triangle (with straight sides) to another, this matter is easy.In the general case it is necessary to investigate whether the transformation is a properone, meaning that the image covers all the triangle T* and has no overlaps.The main point of this section is that the basic ideas of the one variable casecarry over to several variables.