Heath - Scientific Computing (523150), страница 70
Текст из файла (страница 70)
In evaluating a two-dimensional integral, we wishto compute the volume under a surface over a region in the plane. For a rectangular region,258CHAPTER 8. NUMERICAL INTEGRATION AND DIFFERENTIATIONa double integral has the formZ bZadf (x, y) dx dy.cFor a more general two-dimensional domain Ω, the integral takes the formZZf (x, y) dA.ΩBy analogy with numerical quadrature for one-dimensional integrals, the numerical approximation of two-dimensional integrals is sometimes called numerical cubature.To evaluate a double integral, a number of approaches are available, including the following:• Use an automatic one-dimensional quadrature routine for each dimension, one for theouter integral and the other for the inner integral. Each time the outer routine callsits integrand function, the latter will call the inner quadrature routine.
This approachrequires some care in setting the error tolerances for the respective quadrature routines.• Use a product quadrature rule, which results from applying a one-dimensional rule tosuccessive dimensions. This approach is limited to standard domains, such as rectangles.• Use a nonproduct quadrature rule. In recent years, such rules, including error estimates,have become available. The most important case for automatic adaptive use is for triangles, since many two-dimensional regions can be efficiently triangulated to any desireddegree of refinement.8.5.4Multiple IntegralsTo evaluate a multiple integral in dimensions higher than two, the only generally viableapproach is the Monte Carlo method.
The function is sampled at n points distributedrandomly in the domain of integration, and then the mean of these function values ismultiplied by the area (or volume, etc.) of the domain to obtain an estimate for theintegral. The error in this estimate goes to zero as n−1/2 , which means, for example, thatto gain an extra decimal place of accuracy the number of sample points must be increasedby a factor of 100. For this reason, it is not unusual for Monte Carlo calculations of integralsto require millions of evaluations of the integrand.The Monte Carlo method is not competitive for integrals in one or two dimensions,but the beauty of the method is that its convergence rate is independent of the number ofdimensions.
Thus, for example, one million points in six dimensions amounts to only tenpoints per dimension, which is vastly better than any type of conventional quadrature rulewould require for the same level of accuracy. The efficiency of Monte Carlo integration canbe enhanced by various methods for biasing the sampling, either to achieve more uniformcoverage of the sampled volume (e.g., by avoiding undesirable random clumping of thesample points; see Section 13.4) or to concentrate sampling in regions where the integrandis largest in magnitude (importance sampling) or in variability (stratified sampling), in aspirit similar to adaptive quadrature.
See Chapter 13 for further information on the use ofrandom sampling for numerical integration as well as other types of problems.8.6. INTEGRAL EQUATIONS8.6259Integral EquationsAn integral equation is an equation in which the unknown to be determined is a functioninside an integral sign. An integral equation can be thought of as a continuous analogue,or limiting case, of a system of algebraic equations. For example, the analogue of a linearsystem Ax = y is a Fredholm integral equation of the first kind, which has the formZbK(s, t)u(t) dt = f (s),awhere the functions K, called the kernel , and f are known, and the function u is to bedetermined. Integral equations arise naturally in many fields of science and engineering,particularly observational sciences (e.g., astronomy, seismology, spectrometry), where thekernel K represents the response function of an instrument (determined by calibration withknown signals), f represents measured data, and u represents the underlying signal thatis sought.
In effect, we are trying to resolve the measured data f as a (continuous) linearcombination of standard signals. Integral equations can also result from Green’s functionmethods [214] or boundary element methods [154] for solving differential equations (topicsbeyond the scope of this book).Establishing the existence and uniqueness of solutions to integral equations is much moreproblematic than with algebraic equations. Moreover, when a solution does exist, it may beextremely sensitive to perturbations in the input data f , which are often subject to randomexperimental or measurement errors. The reason for this sensitivity is that integration is asmoothing process, so its inverse (i.e., determining the integrand from the integral) is justthe opposite.
Integrating an arbitrary function u against a smooth kernel K dampens anyhigh-frequency oscillation, so solving for u tends to introduce high-frequency oscillation inthe result. For example, Riemann showed that for any integrable kernel K,limZn→∞ abK(s, t) sin(nt) dt = 0,which implies that an arbitrarily high-frequency component of u has an arbitrarily smalleffect on f . Thus, integral equations of the first kind with smooth kernels are alwaysill-conditioned.A standard technique for solving integral equations numerically is to use a quadratureformula to replace the integral by an approximating finite sum. Denote the nodes andweights of the quadrature rule by tj and wj , j = 1, .
. . n. We also choose n points si for thevariable s, often the same as the tj , but not necessarily so. Then the approximation to theintegral equation becomesnXK(si , tj )wj u(tj ) = f (si ),i = 1, . . . n.j=1This system of linear algebraic equations Ax = y, where aij = K(si , tj )wj , yi = f (si ), andxj = u(tj ), can then be solved for x to obtain a discrete sample of approximate values ofthe function u.260CHAPTER 8. NUMERICAL INTEGRATION AND DIFFERENTIATIONExample 8.6 Integral Equation. Consider the integral equationZ1(1 + αst)u(t) dt = 1−1[i.e., K(s, t) = 1 + αst and f (s) = 1], where α is a known positive constant whose value isunspecified for now. Using the composite midpoint quadrature rule with two panels, takingt1 = − 12 , t2 = 12 , and w1 = w2 = 1, and also taking s1 = − 12 and s2 = 12 , we obtain thelinear system 1 + α/4 1 − α/4x11=.1 − α/4 1 + α/4x21It is easily verified that the solution to this linear system is x = [ 12 12 ]T , independent ofthe value of α.Now suppose that the measured values of y1 = f (s1 ) and y2 = f (s2 ) are in error by 1and 2 , respectively.
Then by linearity, the change in the solution x is given by the samelinear system, but with a right-hand side of [ 1 2 ]T . The resulting change in x is thereforegiven by ∆x1(1 − 2 )/α + (1 + 2 )/4=.∆x2(2 − 1 )/α + (1 + 2 )/4Thus, if α is sufficiently small, the relative error in the computed value for x can be arbitrarily large. A very small value for α in this particular kernel corresponds to a very insensitiveinstrument with a very flat response. This is reflected in the conditioning of the matrix A,whose columns become more nearly linearly dependent as α decreases in magnitude. Thissimple example is typical of integral equations with smooth kernels.Note that the sensitivity in the previous example is inherent in the problem and is notdue to the method of solving it. In general, such an integral operator with a smooth kernelhas zero as an eigenvalue (i.e., there are nonzero functions that it annihilates), and henceusing a more accurate quadrature rule makes the conditioning of the linear system worseand the resulting solution more erratic.
Because of this behavior, additional informationmay be required to obtain a physically meaningful solution. Such techniques include:• Truncated singular value decomposition. The solution to the system Ax = y is computedusing the SVD of A; but the small singular values of A, which reflect the ill-conditioning,are omitted from the solution (see Section 4.5.2).• Regularization. A damped solution is obtained by solving the minimization problemmin(ky − Axk22 + µkxk22 ),xwhere the parameter µ determines the relative weight given to the norm of the residualand the norm of the solution. This minimization problem is equivalent to the linear leastsquares problem Ay√x≈,µIo8.7.
NUMERICAL DIFFERENTIATION261which can be solved by the methods discussed in Chapter 3. More generally, other norms,usually based on first or second differences between its components, can also be used toweight the smoothness of the solution. The Levenberg-Marquardt method for nonlinearleast squares problems (see Section 6.4.2) is another example of regularization.• Constrained optimization. Some norm of the residual ky − Axk is minimized subjectto constraints on x that disallow nonphysical solutions. In many applications, for example, the components of the solution x are required to be nonnegative. The resultingconstrained optimization problem can then be solved by one of the methods discussed inSection 6.5.A variety of such methods are implemented in the MATLAB toolbox documented in [121].We have considered only Fredholm integral equations of the first kind. Many othertypes arise in practice, including integral equations of the second kind (eigenvalue problems), Volterra integral equations (in which the upper limit of integration is s instead ofb), singular integral equations (in which one or both of the limits of integration are infinite), and nonlinear integral equations.
All types of integral equations can be discretizedby means of numerical quadrature, yielding a system of algebraic equations.PAlternatively,the unknown function u can be expressed as a linear combination u(t) = nj=1 cj φj (t) ofsuitably chosen basis functions φj , which leads to a system of algebraic equations for thecoefficients cj . This type of approach will be examined in more detail in Section 10.5, whenwe consider finite element methods for boundary value problems in differential equations.8.7Numerical DifferentiationWe now turn briefly to numerical differentiation. It is important to realize that differentiation is an inherently sensitive problem, as small perturbations in the data can cause largechanges in the result.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
