f18-2

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, and gwould 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. 18.2VolterraEquations 787Sample 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 voltra(n,m,t0,h,t,f,g,ak) INTEGER m,n,MMAX REAL h,t0,f(m,n),t(n),g,ak EXTERNAL ak,gPARAMETER (MMAX=5) C USES ak,g,lubksb,ludcmp Solves a set of mlinear Volterra equations of the second kind using the extended trapezoidal rule. On input, t0 is the starting point of the integration and n-1 is the number of steps of size hto be taken. g(k,t) is a user-supplied external function that returns gk(t), while ak(k,l,t,s) is another user-supplied external function that returns the (k, l )element of the matrix K(t, s ). The solution is returned in f(1:m,1:n) , with the corresponding abscissas in t(1:n) . INTEGER i,j,k,l,indx(MMAX) REAL d,sum,a(MMAX,MMAX),b(MMAX) t(1)=t0do 11k=1,m Initialize. f(k,1)=g(k,t(1)) enddo 11 do16i=2,n Take a step h. t(i)=t(i-1)+h do14k=1,m sum=g(k,t(i)) Accumulate right-hand side of linear equations in sum . do13l=1,m sum=sum+0.5*h*ak(k,l,t(i),t(1))*f(l,1) do12j=2,i-1 sum=sum+h*ak(k,l,t(i),t(j))*f(l,j) enddo 12 if(k.eq.l)then Left-hand side goes in matrix a. a(k,l)=1. else a(k,l)=0. endifa(k,l)=a(k,l)-0.5*h*ak(k,l,t(i),t(i)) enddo 13 b(k)=sum enddo 14 call ludcmp(a,m,MMAX,indx,d) Solve linear equations. call lubksb(a,m,MMAX,indx,b) do15k=1,m f(k,i)=b(k) enddo 15 enddo 16 returnEND FornonlinearVolterraequations,equation(18.2.4)holdswiththeproduct K iifi replaced by Kii(fi), and similarly for the other two products of K’s and f’s. Thus for each iwe solve a nonlinear equation for fiwith a known right-hand side. Newton’s method ( §9.4 or §9.6) with an initial guess of fi−1usually works very well provided the stepsize is not too big. Higher-ordermethodsfor solvingVolterraequationsare, in ouropinion,notas important as for Fredholm equations, since Volterra equations are relatively easy to solve. However, there is an extensive literature on the subject. Several difficulties arise. First,anymethodthatachieveshigherorderbyoperatingonseveralquadraturepoints simultaneously will need a special method to get started, when values at the first few points are not yet known. Second, stable quadrature rules can give rise to unexpected instabilities in integral equations. For example, suppose we try to replace the trapezoidal rule in 788 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).the algorithm above with Simpson’s rule. Simpson’s rule naturally integrates over an interval 2h, so we easily get the functionvalues at the even mesh points. For the oddmeshpoints,wecouldtryappendingonepaneloftrapezoidalrule. Buttowhich endoftheintegrationshouldweappendit? Wecoulddoonestepoftrapezoidalrule followed by all Simpson’s rule, or Simpson’s rule with one step of trapezoidal ruleat the end. Surprisingly,the formerscheme is unstable, while the latter is fine! A simple approach that can be used with the trapezoidal method given above is Richardson extrapolation: Compute the solution with stepsize hand h/2. Then, assuming the error scales with h 2, compute fE=4f(h/2)−f(h) 3(18.2.5 ) This procedure can be repeated as with Romberg integration. The general consensus is that the best of the higher order methods is the block-by-block method (see[1]). Another important topic is the use of variable stepsize methods, which are much more efficient if there are sharp features in Kor f. Variablestepsizemethodsarequiteabitmorecomplicatedthantheircounterparts for differentialequations; we refer you to the literature [1,2]for a discussion. Youshouldalso beon thelookoutforsingularitiesin theintegrand. If youfind them, then look to §18.3 for additional ideas. CITED REFERENCES AND FURTHER READING: Linz, P. 1985, Analytical and NumericalMethods for Volterra Equations (Philadelphia:S.I.A.M.). [1] Delves, L.M., and Mohamed, J.L. 1985, Computational Methods for Integral Equations (Cam- bridge, U.K.: Cambridge University Press). [2] 18.3 Integral Equations with Singular Kernels Many integral equations have singularities in either the kernel or the solution or both. A simple quadrature method will show poor convergence with Nif such singularities are ignored. There is sometimes art in how singularities are best handled. We start with a few straightforward suggestions:1. Integrablesingularitiescanoftenberemovedbyachangeofvariable. Forexample,the singular behavior K(t, s)∼s 1/2ors−1/2nears=0can be removed by the transformation z=s1/2. Note that we are assuming that the singular behavior is confined to K, whereas the quadrature actually involves the product K(t, s)f(s), and it is this product that must be “fixed.” Ideally,youmustdeducethesingularnatureoftheproductbeforeyoutryanumericalsolution, and take the appropriate action. Commonly, however, a singular kernel does not produce a singular solution f(t). (The highly singular kernel K(t, s)=δ(t−s)is simply the identity operator, for example.) 2. If K(t, s)can be factored as w(s) K(t, s), where w(s)is singular and K(t, s)is smooth, then aGaussian quadrature based on w(s)asa weightfunction willwork well. Even if the factorization is only approximate, the convergence is often improved dramatically. Allyou havetodoisreplace gaulegintheroutine fred2byanother quadrature routine. Section 4.5 explained how to construct such quadratures; or you can find tabulated abscissas andweights in the standard references [1,2]. You must of course supply Kinstead of K.