f18-1

782 Chapter18. IntegralEquationsandInverseTheorySample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).specialquadraturerules,buttheyarealsosometimesblessingsindisguise,sincethey can spoil a kernel’s smoothing and make problems well-conditioned. In§§18.4–18.7 we face up to the issues of inverse problems. §18.4 is an introduction to this large subject. We should note here that wavelet transforms, already discussed in §13.10, are applicable not only to data compressionand signal processing, but can also be used totransformsomeclassesofintegralequationsintosparselinearproblemsthatallow fast solution. You may wish to review §13.10 as part of reading this chapter. Some subjects, such as integro-differential equations , we must simply declare to be beyondour scope. For a review of methods for integro-differentialequations,see Brunner [4]. Itshouldgowithoutsayingthatthis oneshortchaptercanonlybarelytouchon a few of the most basic methods involved in this complicated subject. CITED REFERENCES AND FURTHER READING: Delves, L.M., and Mohamed, J.L. 1985, Computational Methods for Integral Equations (Cam- bridge, U.K.: Cambridge University Press). [1] Linz, P. 1985, Analytical and NumericalMethods for Volterra Equations (Philadelphia:S.I.A.M.). [2] Atkinson, K.E. 1976, A Survey of Numerical Methods for the Solution of Fredholm Integral Equations of the Second Kind (Philadelphia: S.I.A.M.). [3] Brunner,H.1988,in NumericalAnalysis1987 ,PitmanResearchNotesinMathematicsvol.170, D.F. Griffiths and G.A. Watson, eds. (Harlow, Essex, U.K.: Longman Scientific and Tech-nical), pp. 18–38. [4] Smithies, F. 1958, Integral Equations (Cambridge, U.K.: Cambridge University Press). Kanwal, R.P. 1971, Linear Integral Equations (New York: Academic Press). Green, C.D. 1969, Integral Equation Methods (New York: Barnes & Noble). 18.1 Fredholm Equations of the Second Kind We desire a numerical solution for f(t)in the equation f(t)=λ/integraldisplayb aK(t, s)f(s)ds+g(t)( 18.1.1 ) Themethodwe describe,a verybasic one,is calledthe Nystrom method . It requires the choice of some approximate quadrature rule : /integraldisplayb ay(s)ds=N/summationdisplay j=1wjy(sj)( 18.1.2 ) Here the set {wj}are the weights of the quadrature rule, while the Npoints {sj} are the abscissas. What quadrature rule should we use? It is certainly possible to solve integral equationswithlow-orderquadraturerulesliketherepeatedtrapezoidalorSimpson’s 18.1FredholmEquationsoftheSecondKind 783Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).rules. We will see, however, that the solution method involves O(N3)operations, and so the most efficient methods tend to use high-order quadrature rules to keepNas small as possible. For smooth, nonsingularproblems, nothing beats Gaussian quadrature (e.g., Gauss-Legendre quadrature, §4.5). (For non-smooth or singular kernels, see §18.3.) Delves and Mohamed [1]investigated methods more complicated than the Nystrom method. For straightforward Fredholm equations of the second kind, they concluded“ ...theclearwinnerofthiscontesthasbeentheNystromroutine ...with theN-pointGauss-Legendrerule. Thisroutineisextremelysimple .... Suchresults are enough to make a numerical analyst weep.” If we apply the quadrature rule (18.1.2) to equation (18.1.1), we get f(t)=λN/summationdisplay j=1wjK(t, s j)f(sj)+g(t)( 18.1.3 ) Evaluate equation (18.1.3) at the quadrature points: f(ti)=λN/summationdisplay j=1wjK(ti,sj)f(sj)+g(ti)( 18.1.4 ) Letfibe the vector f(ti),githe vector g(ti),Kijthematrix K(ti,sj), anddefine /tildewideKij=Kijwj (18.1.5 ) Then in matrix notation equation (18.1.4) becomes (1−λ/tildewideK)·f=g (18.1.6 ) This is a set of Nlinear algebraic equations in Nunknowns that can be solved by standard triangular decompositiontechniques ( §2.3) — that is where the O(N3) operations count comes in. The solution is usually well-conditioned, unless λis very close to an eigenvalue. Having obtained the solution at the quadraturepoints {ti}, how do you get the solution at some other point t? You do notsimply use polynomial interpolation. This destroys all the accuracy you have worked so hard to achieve. Nystrom’s key observation was that you should use equation (18.1.3) as an interpolatory formula, maintaining the accuracy of the solution. We here give two subroutines for use with linear Fredholm equations of the second kind. The routine fred2sets up equation (18.1.6)and then solves it by LU decompositionwith calls to the routines ludcmpandlubksb. The Gauss-Legendre quadrature is implemented by first getting the weights and abscissas with a call to gauleg. Routine fred2requires that you provide an external function that returns g(t)and another that returns λK ij. It then returns the solution fat the quadrature points. It also returns the quadrature points and weights. These are used by the second routine fredinto carry out the Nystrom interpolation of equation (18.1.3) and return the value of fat any point in the interval [a, b ]. 784 Chapter18. IntegralEquationsandInverseTheorySample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE fred2(n,a,b,t,f,w,g,ak) INTEGER n,NMAX REAL a,b,f(n),t(n),w(n),g,ak EXTERNAL ak,gPARAMETER (NMAX=200) C USES ak,g,gauleg,lubksb,ludcmp Solves a linear Fredholm equation of the second kind. On input, aandbare the limits of integration, and nis the number of points to use in the Gaussian quadrature. gandak are user-supplied external functions that respectively return g(t)andλK (t, s ). The routine returns arrays t(1:n) andf(1:n) containing the abscissas tiof the Gaussian quadrature and the solution fat these abscissas. Also returned is the array w(1:n) of Gaussian weights for use with the Nystrom interpolation routine fredin . INTEGER i,j,indx(NMAX) REAL d,omk(NMAX,NMAX) if(n.gt.NMAX) pause ’increase NMAX in fred2’call gauleg(a,b,t,w,n) Replace gauleg with another routine if not using Gauss-Legendre quadrature. do 12i=1,n do11j=1,n Form1−λ/tildewideK. if(i.eq.j)then omk(i,j)=1. else omk(i,j)=0. endif omk(i,j)=omk(i,j)-ak(t(i),t(j))*w(j) enddo 11 f(i)=g(t(i)) enddo 12 call ludcmp(omk,n,NMAX,indx,d) Solve linear equations. call lubksb(omk,n,NMAX,indx,f)return END FUNCTION fredin(x,n,a,b,t,f,w,g,ak) INTEGER n REAL fredin,a,b,x,f(n),t(n),w(n),g,akEXTERNAL ak,g C USES ak,g Given arrays t(1:n) andw(1:n) containing the abscissas and weights of the Gaussian quadrature, and given the solution array f(1:n) from fred2 , this function returns the value of fatxusing the Nystrom interpolation formula. On input, aandbare the limits of integration, and nis the number of points used in the Gaussian quadrature. gandak are user-supplied external functions that respectively return g(t)andλK (t, s ). INTEGER iREAL sum sum=0. do 11i=1,n sum=sum+ak(x,t(i))*w(i)*f(i) enddo 11 fredin=g(x)+sum returnEND One disadvantageof a methodbasedon Gaussian quadratureis that thereis no simplewaytoobtainanestimateoftheerrorintheresult. Thebestpracticalmethod istoincrease Nby50%,say,andtreatthedifferencebetweenthetwoestimatesasa conservativeestimate of theerrorin theresult obtainedwith thelargervalueof N. 18.1FredholmEquationsoftheSecondKind 785Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Turn now to solutions of the homogeneous equation. If we set λ=1/σand g=0, then equation (18.1.6) becomes a standard eigenvalue equation /tildewideK·f=σf (18.1.7 ) which we can solve with any convenient matrix eigenvalue routine (see Chapter 11). Note that if our original problem had a symmetric kernel, then the matrix K is symmetric. However, since the weights wjare not equal for most quadrature rules, the matrix /tildewideK(equation 18.1.5) is not symmetric. The matrix eigenvalue problem is much easier for symmetric matrices, and so we should restore the symmetryifpossible. Providedtheweightsarepositive(whichtheyareforGaussianquadrature), we can define the diagonal matrix D=diag (w j)and its square root, D1/2=diag (√wj). Then equation (18.1.7) becomes K·D·f=σf Multiplying by D1/2, we get /parenleftBig D1/2·K·D1/2/parenrightBig ·h=σh (18.1.8 ) whereh=D1/2·f. Equation(18.1.8)is nowin theformofa symmetriceigenvalue problem. Solution of equations (18.1.7) or (18.1.8) will in general give Neigenvalues, where Nis the number of quadrature points used. For square-integrable kernels, these will provide good approximations to the lowest Neigenvalues of the integral equation. Kernels of finite rank (also called degenerate orseparable kernels) have only a finite number of nonzero eigenvalues (possibly none). You can diagnosethis situation by a cluster of eigenvalues σthat are zero to machine precision. The number of nonzero eigenvalues will stay constant as you increase Nto improve their accuracy. Some care is required here: A nondegenerate kernel can have aninfinite number of eigenvalues that have an accumulation point at σ=0.Y o u distinguish the two cases by the behavior of the solution as you increase N. If you suspectadegeneratekernel,youwillusuallybeabletosolvetheproblembyanalytic techniques described in all the textbooks. CITED REFERENCES AND FURTHER READING: Delves, L.M., and Mohamed, J.L. 1985, Computational Methods for Integral Equations (Cam- bridge, U.K.: Cambridge University Press). [1] Atkinson, K.E. 1976, A Survey of Numerical Methods for the Solution of Fredholm Integral Equations of the Second Kind (Philadelphia: S.I.A.M.). 786 Chapter18. IntegralEquationsandInverseTheorySample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).18.2 Volterra Equations Let us now turn to Volterra equations, of which our prototype is the Volterra equation of the second kind, f(t)=/integraldisplayt aK(t, s)f(s)ds+g(t)( 18.2.1 ) MostalgorithmsforVolterraequationsmarchoutfrom t=a,buildingupthesolution as they go. In this sense they resemble not only forward substitution (as discussed in§18.0),but also initial-valueproblemsfor ordinarydifferentialequations. In fact, many algorithms for ODEs have counterparts for Volterra equations. The simplest way to proceed is to solve the equation on a mesh with uniform spacing: ti=a+ih, i =0,1,...,N, h ≡b−a N(18.2.2 ) To do so, we must choose a quadrature rule. For a uniform mesh, the simplest scheme is the trapezoidal rule, equation (4.1.11): /integraldisplayti aK(ti,s)f(s)ds=h 1 2Ki0f0+i−1/summationdisplay j=1Kijfj+1 2Kiifi  (18.2.3 ) Thus the trapezoidal method for equation (18.2.1) is: f0=g0 (1−1 2hK ii)fi=h 1 2Ki0f0+i−1/summationdisplay j=1Kijfj +gi,i =1,...,N(18.2.4 ) (For a Volterra equation of the first kind, the leading 1on the left would be absent, andgwould have opposite sign, with correspondingstraightforwardchanges in the rest of the discussion.) Equation (18.2.4) is an explicit prescription that gives the solution in O(N2) operations. UnlikeFredholmequations,itisnotnecessarytosolveasystemoflinear equations. Volterra equationsthus usually involveless work than the correspondingFredholmequationswhich, as we haveseen, do involvethe inversionof,sometimes large, linear systems. The efficiency of solving Volterra equations is somewhat counterbalanced by the fact that systemsof these equations occur more frequently in practice. If we interpret equation (18.2.1) as a vectorequation for the vector of mfunctions f(t), then the kernel K(t, s)is an m×mmatrix. Equation (18.2.4) must now also be understood as a vector equation. For each i, we have to solve the m×mset of linear algebraic equations by Gaussian elimination. The routine voltrabelow implements this algorithm. You must supply an external function that returns the kth function of the vector g(t)at the point t, and another that returns the (k,l)element of the matrix K(t, s)at(t, s). The routine voltrathen returns the vector f(t)at the regularly spaced points t i.