Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 26
Текст из файла (страница 26)
Namely, first approximate f(x) with a convenientfunction S(x) and then compute analyticallyas an approximation to the desired integral. In detail there are important differencesbecause here the weight function w(x) does not have one sign and oscillates rapidlyfor large frequencies ω. Also, the approach of Chapter 5 would apply to particularω, and we would like a method that can be applied conveniently for any ω. Insightis provided by a classic technique of applied mathematics for approximating Fourierintegrals when w is “large.” If derivatives of f(x) are available, asymptotic approximations can be obtained by means of integration by parts. For example, integrating byparts twice giveswhereIf M1 is a bound on the magnitude of f´(x), thenthat is, R(ω) is 0(ω−2 ) asAccordingly, the asymptotic approximation3.7 CASE STUDY 3129is accurate to 0(ω-2).
However, the integral itself ordinarily goes to zero like ω - 1 ,so the relative error is ordinarily only 0( ω-1). The situation is typical of classicalasymptotic approximation of integrals. The bigger w is, the better the asymptotic approximation and the more difficult the integral is for conventional numerical methods.On the other hand, for a given w, the asymptotic approximation may not be sufficientlyaccurate and there is no easy way to improve it. When w is “small,” finite Fourier integrals are easy for conventional numerical methods, and when w is “large,” asymptoticapproximations are satisfactory.
Filon’s method provides a way to compute accurateintegrals when w is of moderate size.Filon divides the interval [a,b] into 2N subintervals of equal length h. Let us definexj = a + jh for j = 0, . . . ,2N. The function f(x) is approximated by a continuous splineS(x) that is a quadratic polynomial on each [x 2 m , x m + 2 ] defined there by interpolationto f(xj) for j = 2m, 2m + 1,2m + 2. Each of the integrals in the approximationcan be evaluated analytically by integration by parts.
A lengthy calculation results inFilon’s method:Here θ = ωh andAlso,There is a similar approximation when cosine is replaced by sine. The formula isThe coefficients α, β, γ are the same, and Se and So are like Ce and Co with the cosinesreplaced by sines.Using the results developed in this chapter for the error of interpolation, it is easyto bound the error of Filon’s method.
On each subinterval [x 2 m , x 2 m + 2 ], S(x) is aquadratic interpolant to f(x). If M3 is a bound on the third derivative of f(x) onall of [a,b], we found that130CHAPTER 3INTERPOLATIONuniformly for a < x < b. ThenThis is a bound on the absolute error that depends on f(x) and the sampling intervalh, but not on ω. If we want a meaningful result for “large” ω, we have to take intoaccount that the integral is 0(ω -l). This leads to a bound on the relative error that is-1 40( θ h ).In a subsection of Chapter 5 about applying a general-purpose code to problemswith oscillatory integrands, the exampleis discussed. The function f(x) = l/(1 + x2 ) ought to be approximated well by aquadratic interpolatory spline with a relatively crude mesh, so Filon’s method oughtto be quite effective.
A matter not usually discussed with Filon’s method is how toestimate the accuracy of the result. One way to proceed is to compute a result that webelieve to be more accurate and estimate the error of the less accurate result by comparison. If the accuracy is acceptable, often the more accurate result is the one taken asthe answer.
Inspection of the formulas shows that if h, or equivalently θ, is halved, wecan reuse all the evaluations of f, sin(x), and cos(x) made in the first approximationto keep down the cost. According to the bounds, the error should be reduced enoughby halving h to get a good estimate of the error by comparison. Using Filon’s methodwith θ = 0.4, we obtained an approximate integral of 0.04566373122996.
Halving θresulted in the approximation 0.04566373690838. Estimating the accuracy by comparison suggests that we have an answer with an error smaller than 6 × 10-9. Reuseof the function evaluations by virtue of halving θ holds the cost to 315 evaluationsof the integrand. The quadrature code Adapt developed in Chapter 5 asks the userto specify the accuracy desired. Using the code in the manner outlined in that chapter and experimenting some with the error tolerances, we obtained an approximation0.04566373866324 with an estimated error of about -1.6 × 10-9.
This cost 1022evaluations of the integrand.Let us now change the subject from applying a continuous spline to applying asmooth spline. In this chapter we have been looking at the approximation of functionsy = f(x), but sometimes we want to approximate curves. This will be discussed inthe plane for it will then be clear how to deal with a curve in three dimensions. Thebasic idea is to use a parametric representation (x(s),y(s)) for the curve and approximate independently the coordinate functions x(s) and y(s).
The parameter s can beanything, but often in the theory it is taken to be arc length. Having chosen somehownodes si , i = l, . . . , N, we can interpolate the data xi = x(si ) by a spline Sx(s) and likewise the data yi = y(si ) by a spline Sy(s). The curve (x(s),y(s)) is then approximatedby (Sx(s),Sy(s)). This yields a curve in the plane that passes through all the points(x(si ),y(si )) in order. It is natural to use the smooth cubic spline of SPCOEF becauseit leads to a curve with continuous curvature, but if the data are sparse, we might haveto resort to the shape-preserving spline to get a curve of the expected shape.
All thesecomputations are familiar except for the selection of the nodes si . One way to proceedREFERENCES131Figure 3.21 Curve fit for a sewing machine pattern.is to choose them so that the parameter s approximates the arc length of the curve.This is done by taking s1 = 0 and defining the difference between si and si+1 as thedistance between the points (xi ,yi ) and (x i+1 ,y i+l ), namelyExercise 3.39 suggests an alternative. See Farin [6] for many other approaches.An interesting example of the technique is furnished by a need to approximatecurves for automatic control of sewing machines. Arc length is the natural parameterbecause a constant increment in arc length corresponds to a constant stitch length.An example taken from [18] fits the data (2.5, -2.5), (3.5, -0.5), (5, 2), (7.5, 4),(9.5, 4.5), (11.8, 3.5), (13, 0.5), (11.5, -2), (9, -3), (6, -3.3), (2.5, -2.5), (0, 0),(-1.5, 2), (-3, 5), (-3.5, 9), (-2, 11), (0, 11.5), (2, 11), (3.5, 9), (3, 5), (1.5, 2),(0, 0), (-2.5, -2.5), (-6, -3.3), (-9, -3), (-11.5, -2), (-13, 0.5), (-11.8, 3.5),(-9.5, 4.5), (-7.5, 4), (-5, 2), (-3.5, -0.5), (-2.5, -2.5).
The resulting spline curveis seen in Figure 3.21.nREFERENCES1. E. Becker, G. Carey and J. T. Oden, Finite Elements: An Introduction, Vol. I, Prentice Hall,Englewood Cliffs, N.J., 1981.2. G. Birkhoff and A. Priver, “Hermite interpolation errors for derivatives,” J. Math. and Physics,46 (1967), pp. 440-447.132CHAPTER 3 INTERPOLATION3.
C. de Boor, A Practical Guide to Splines, Springer-Verlag, New York, 1978.4. E. W. Cheney, Introduction to Approximation Theory, McGraw-Hill, New York, 1966.5. J. Crank and G. Park, “Evaluation of the diffusion coefficient for CHCl3 in polystyrene fromsimple absorption experiments,” Trans. Faraday Soc., 45 (1949), pp. 240-249.6. G. Farin, Curves and Surfaces for Computer Aided Geometric Design, Academic Press, SanDiego, 1988.7. J. Ferguson and K.
Miller, “Characterization of shape in a class of third degree algebraic curves,”TRW Report 5322-3-5, 1969.8. L. N. G. Filon, “On a quadrature formula for trigonometric integrals,” Proc. Roy. Soc. Edinburgh,49 (1928-1929), pp. 38-47.9. F. Fritsch and J.
Butland, “A method for constructing local monotone piecewise cubic interpolants,” SIAM J. Sci. Stat. Comp., 5 (1984), pp. 300-304.10. F. Fritsch and R. Carlson, “Monotone piecewise cubic interpolation,” SIAM J. Numer Anal., 17(1980), pp. 238-246.11. C. Hall and W. Meyer, “Optimal error bounds for cubic spline interpolation,” J. Approx. Theory,16 (1976), pp. 105-122.12. Handbook of Chemistry and Physics, 63rd ed., CRC Press, Cleveland, 1982-1983.13.
Handbook of Mathematical Functions, M. Abramowitz and I. Stegun, eds., Dover, Mineola,N.Y., 1964.14. H. Huynh, “Accurate monotone cubic interpolation,” SIAM J. Nume,: Anal,, 30 (1993), pp. 57100.15. E. Isaacson and H. Keller, Analysis of Numerical Methods, Dover, Mineola, N.Y., 1994.16. W. J. Kammerer and G. W. Reddien, “Local convergence of smooth cubic spline interpolates,”SIAM J. Numer Anal., 9 (1972), pp. 687-694.17.
M. J. D. Powell, “On the maximum errors of polynomial approximation defined by interpolationand by least squares criteria,” Comp. J., 9 (1967), pp. 404-407.18. P Rentrop and W. Wever, “Interpolation algorithms for the control of a sewing machine,” Proceedings of ECM II in Oberwolfach, H. Neunzert, ed., Teubner-Kluwer, Dortdrecht, Netherlands, 1988, pp. 251-268.19. B. Savage, and J. Mathis, “Observed properties of interstellar dust,” in Annual Review of Astronomy and Astrophysics, 17 (1979).MISCELLANEOUS EXERCISES FOR CHAPTER 33.39 As mentioned in the case study, there are other ways to select the nodes whenfitting data with curves.
A simple one ist1 = 0ti+l = ti + |xi+l - xi | + |yi+l - yi | 1 < i < N.Using this scheme and SPCOEF or Spline_coeff, fit the data133Using SVALUE or Spline_value, evaluate at sufficiently many points to sketch asmooth curve in the xy plane.3.40 Find a technique that gives a good fit to the model rocket data from Example 3.11.Interpolate the indicated (>) data and avoid undesirable oscillations.3.41 Implement the shape-preserving cubic spline described in Section 3.5. Test it outon some of the data sets from this chapter, for example, Exercises 3.22, 3.25,and 3.27.
Sketch the graph of H(x). Also try it on the nonmonotone data inExample 3.11, Exercise 3.26, and Exercise 3.28. How well does it do?3.42 Implement the bilinear interpolating function Q(x,y) given in (3.48). Test it onseveral different functions and several different grids.3.43 Implement the linear interpolating function Q(n,y) given in (3.50). Test it onseveral different functions and several different triangulations.Previous Home NextCHAPTER 4ROOTS OF NONLINEAR EQUATIONSFinding solutions of a system of nonlinear equationsf(x) = 0(4.1)is a computational task that occurs frequently both on its own and as a part of a morecomplicated problem. Most of this chapter is devoted to the case of a continuous realfunction f(x) of a single real variable x because it is important and can be discussedin elementary terms. The general case of n nonlinear equations in n unknowns ismuch more difficult both in theory and practice.
Although the theory is too involvedto be developed here, some of the basic methods are discussed briefly at the end of thechapter.A root of (4.l), or a zero of f(x), is a number a such that f(α) = 0. A root isdescribed more fully by its multiplicity m.
This means that for x near α, f(x) can bewritten in the form(4.2)f(x) = (x - α) mg(x)where g(x) is continuous near α and g(α)0. If m = 1, the root is said to be simpleand otherwise, multiple. The basic definition permits m to be a fraction. For example,with the functionequation (4.1) has α = 1 as a root of multiplicity l/2 (and α = 0 as a simple root).However, if f(x) is sufficiently smooth, then m must be a positive integer. Indeed, iff(x) has its first m derivatives continuous on an interval that includes a and(4.3)then a is a root of multiplicity m.