f18-3

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 handh/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. 18.3IntegralEquationswithSingularKernels 789Sample 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).Thismethodisaspecialcaseofthe productNystrommethod [3,4],whereonefactorsout a singular term p(t, s)depending on both tandsfromKand constructs suitable weights for its Gaussian quadrature. The calculations in the general case are quite cumbersome, becausethe weights depend on the chosen {t i}as well as the form of p(t, s). WeprefertoimplementtheproductNystrommethodonauniformgrid,withaquadrature scheme that generalizes the extended Simpson’s 3/8 rule (equation 4.1.5) to arbitrary weightfunctions. We discuss this in the subsections below. 3. Special quadrature formulas are also useful when the kernel is not strictly singular, butis“almost”so. Oneexampleiswhenthekernelisconcentratednear t=sonascalemuch smaller than the scale on which the solution f(t)varies. In that case, a quadrature formula can be based on locally approximating f(s)by a polynomial or spline, while calculating the first few moments of the kernel K(t, s)at the tabulation points t i. In such a scheme the narrow width of the kernel becomes an asset, rather than a liability: The quadrature becomesexact as the width of the kernel goes to zero. 4. Aninfiniterangeofintegration isalsoaformofsingularity. Truncatingtherange ata large finite value should be used only as a last resort. If the kernel goes rapidly to zero, then a Gauss-Laguerre [ w∼exp(−αs)] or Gauss-Hermite [ w∼exp(−s 2)] quadrature should work well. Long-tailed functions often succumb to the transformation s=2α z+1−α (18.3.1 ) which maps 0<s< ∞to1>z> −1so that Gauss-Legendre integration can be used. Here α> 0is a constant that you adjust to improve the convergence. 5. A common situation in practice is that K(t, s)is singular along the diagonal line t=s. HeretheNystrommethodfailscompletelybecausethekernelgetsevaluatedat (ti,si). Subtraction of the singularity is one possible cure: /integraldisplayb aK(t, s)f(s)ds=/integraldisplayb aK(t, s)[f(s)−f(t)]ds+/integraldisplayb aK(t, s)f(t)ds =/integraldisplayb aK(t, s)[f(s)−f(t)]ds+r(t)f(t)(18.3.2 ) where r(t)=/integraltextb aK(t, s)dsis computed analytically or numerically. If the first term on the right-hand side is now regular, we can use the Nystrom method. Instead of equation(18.1.4), we get f i=λN/summationdisplay j=1 j/negationslash=iwjKij[fj−fi]+λrifi+gi (18.3.3 ) Sometimesthesubtractionprocessmustberepeatedbeforethekerneliscompletelyregularized. See[3]for details. (And read on for a different, we think better, way to handle diagonal singularities.) Quadrature ona UniformMesh withArbitraryWeight It is possible in general to find n-point linear quadrature rules that approximate the integral of a function f(x), times an arbitrary weight function w(x), over an arbitrary range ofintegration (a, b),asthesumofweightstimes nevenlyspaced valuesofthefunction f(x), sayat x=kh, (k+1 )h ,..., (k+n−1)h. Thegeneralschemeforderivingsuchquadrature rules is to write down the nlinear equations that must be satisfied if the quadrature rule is to be exact for the nfunctions f(x)=const,x ,x2,...,xn−1, and then solve these for the coefficients. This can be done analytically, once and for all, if the moments of the weightfunction over the same range of integration, W n≡1 hn/integraldisplayb axnw(x)dx (18.3.4 ) 790 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).are assumed to be known. Here the prefactor h−nis chosen to make Wnscale as hif (as in the usual case) b−ais proportional to h. Carrying out this prescription for the four-point case gives the result /integraldisplayb aw(x)f(x)dx= 1 6f(kh)/bracketleftbigg (k+1 ) ( k+2 ) ( k+3 ) W0−(3k2+1 2 k+ 11) W1+3 ( k+2 ) W2−W3/bracketrightbigg +1 2f([k+1 ] h)/bracketleftbigg −k(k+2 ) ( k+3 ) W0+( 3 k2+1 0 k+6 ) W1−(3k+5 ) W2+W3/bracketrightbigg +1 2f([k+2 ] h)/bracketleftbigg k(k+1 ) ( k+3 ) W0−(3k2+8 k+3 ) W1+( 3 k+4 ) W2−W3/bracketrightbigg +1 6f([k+3 ] h)/bracketleftbigg −k(k+1 ) ( k+2 ) W0+( 3 k2+6 k+2 ) W1−3(k+1 ) W2+W3/bracketrightbigg (18.3.5 ) While the terms in brackets superficially appear to scale as k2, there is typically cancellation at both O(k2)andO(k). Equation (18.3.5) can be specialized to various choices of (a, b). The obvious choice isa=kh,b=(k+3 )h, in which case we get a four-point quadrature rule that generalizes Simpson’s 3/8 rule (equation 4.1.5). In fact, we can recover this special case by settingw(x)=1, in which case (18.3.4) becomes W n=h n+1[(k+3 )n+1−kn+1]( 18.3.6 ) The four terms in square brackets equation (18.3.5) each become independent of k, and (18.3.5) in fact reduces to /integraldisplay(k+3)h khf(x)dx=3h 8f(kh)+9h 8f([k+1]h)+9h 8f([k+2]h)+3h 8f([k+3]h)(18.3.7 ) Back to the case of general w(x), some other choices for aandbare also useful. For example, we may want to choose (a, b)to be ([k+1 ]h,[k+3 ]h)or([k+2 ]h,[k+3 ]h), allowing us to finish off an extended rule whose number of intervals is not a multipleof three, without loss of accuracy: The integral will be estimated using the four valuesf(kh),...,f ([k+3 ]h). Evenmoreusefulistochoose (a, b )tobe ([k+1 ]h,[k+2 ]h),thus using four points to integrate a centered single interval. These weights, when sewed together intoanextended formula,givequadrature schemes thathavesmooth coefficients,i.e.,withouttheSimpson-like 2,4,2,4,2alternation. (Infact,thiswasthetechniquethatweusedtoderive equation 4.1.14, which you may now wish to reexamine.) All these rules are of the same order as the extended Simpson’s rule, that is, exact forf(x)a cubic polynomial. Rules of lower order, if desired, are similarly obtained. The three point formula is /integraldisplay b aw(x)f(x)dx=1 2f(kh)/bracketleftbigg (k+1 ) ( k+2 )W0−(2k+3 )W1+W2/bracketrightbigg +f([k+1 ]h)/bracketleftbigg −k(k+2 )W0+2 (k+1 )W1−W2/bracketrightbigg +1 2f([k+2 ]h)/bracketleftbigg k(k+1 )W0−(2k+1 )W1+W2/bracketrightbigg(18.3.8 ) Here the simple special case is to take, w(x)=1, so that Wn=h n+1[(k+2 )n+1−kn+1]( 18.3.9 ) Then equation (18.3.8) becomes Simpson’s rule, /integraldisplay(k+2)h khf(x)dx=h 3f(kh)+4h 3f([k+1 ]h)+h 3f([k+2 ]h)( 18.3.10 ) 18.3IntegralEquationswithSingularKernels 791Sample 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).For nonconstant weight functions w(x), however, equation (18.3.8) gives rules of one order less than Simpson, since they do not benefit from the extra symmetry of the constant case. The two point formula is simply /integraldisplay(k+1)h khw(x)f(x)dx=f(kh)[(k+1 )W0−W1]+f([k+1 ]h)[−kW 0+W1](18.3.11 ) Here is a routine wwghtsthat uses the above formulas to return an extended N-point quadrature rule for the interval (a, b)=( 0 ,[N−1]h). Input to wwghtsis a user-supplied routine, kermom,thatiscalledtogetthefirstfour indefinite-integral momentsof w(x),namely Fm(y)≡/integraldisplayy smw(s)ds m =0,1,2,3( 18.3.12 ) (The lower limit is arbitrary and can be chosen for convenience.) Cautionary note: When called with N< 4,wwghtsreturns a rule of lower order than Simpson; you should structure your problem to avoid this. SUBROUTINE wwghts(wghts,n,h,kermom) INTEGER nREAL wghts(n),h EXTERNAL kermom C USES kermom Constructs in wghts(1:n) weights for the n-point equal-interval quadrature from 0to (n−1)hof a function f(x)times an arbitrary (possibly singular) weight function w(x)whose indefinite-integral moments Fn(y)are provided by the user-supplied subroutine kermom . INTEGER j,kDOUBLE PRECISION wold(4),wnew(4),w(4),hh,hi,c,fac,a,b hh=h Double precision on internal calculations even though the interface is in single precision. hi=1.d0/hh do 11j=1,n Zero all the weights so we can sum into them. wghts(j)=0. enddo 11 call kermom(wold,0.d0,4) Evaluate indefinite integrals at lower end. if (n.ge.4) then Use highest available order. b=0.d0 For another problem, you might change this lower limit. do14j=1,n-3 c=j-1 This is called kin equation (18.3.5). a=b Set upper and lower limits for this step. b=a+hhif (j.eq.n-3) b=(n-1)*hh Last interval: go all the way to end. call kermom(wnew,b,4) fac=1.d0 do 12k=1,4 Equation (18.3.4). w(k)=(wnew(k)-wold(k))*facfac=fac*hi enddo 12 wghts(j)=wghts(j)+ Equation (18.3.5). * ((c+1.d0)*(c+2.d0)*(c+3.d0)*w(1) * -(11.d0+c*(12.d0+c*3.d0))*w(2) * +3.d0*(c+2.d0)*w(3)-w(4))/6.d0 wghts(j+1)=wghts(j+1)+ * (-c*(c+2.d0)*(c+3.d0)*w(1) * +(6.d0+c*(10.d0+c*3.d0))*w(2) * -(3.d0*c+5.d0)*w(3)+w(4))*.5d0 wghts(j+2)=wghts(j+2)+ * (c*(c+1.d0)*(c+3.d0)*w(1) * -(3.d0+c*(8.d0+c*3.d0))*w(2)* +(3.d0*c+4.d0)*w(3)-w(4))*.5d0 wghts(j+3)=wghts(j+3)+ * (-c*(c+1.d0)*(c+2.d0)*w(1) * +(2.d0+c*(6.d0+c*3.d0))*w(2)* -3.d0*(c+1.d0)*w(3)+w(4))/6.d0 do 13k=1,4 Reset lower limits for moments. 792 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).wold(k)=wnew(k) enddo 13 enddo 14 else if (n.eq.3) then Lower-order cases; not recommended. call kermom(wnew,hh+hh,3) w(1)=wnew(1)-wold(1) w(2)=hi*(wnew(2)-wold(2))w(3)=hi**2*(wnew(3)-wold(3))wghts(1)=w(1)-1.5d0*w(2)+0.5d0*w(3) wghts(2)=2.d0*w(2)-w(3) wghts(3)=0.5d0*(w(3)-w(2)) else if (n.eq.2) then call kermom(wnew,hh,2) wghts(2)=hi*(wnew(2)-wold(2)) wghts(1)=wnew(1)-wold(1)-wghts(2) endif END We willnow give an example of how to apply wwghtsto a singular integral equation. Worked Example: ADiagonallySingular Kernel As a particular example, consider the integral equation f(x)+/integraldisplayπ 0K(x, y )f(y)dy=s i n x (18.3.13 ) with the (arbitrarily chosen) nasty kernel K(x, y )=c o s xcosy×/braceleftbigg −ln(x−y)y<x√y−xy ≥x(18.3.14 ) which has a logarithmic singularity on the left of the diagonal, combined with a square-root discontinuity on the right. The first step is to do (analytically, in this case) the required moment integrals over the singular part of the kernel, equation (18.3.12). Since these integrals are done at a fixedvalue of x, we can use xas the lower limit. For any specified value of y, the required indefinite integral is then either F m(y;x)=/integraldisplayy xsm(s−x)1/2ds=/integraldisplayy−x 0(x+t)mt1/2dtify>x (18.3.15 ) or Fm(y;x)=−/integraldisplayy xsmln(x−s)ds=/integraldisplayx−y 0(x−t)mlntd tify<x (18.3.16 ) (where a change of variable has been made in the second equality in each case). Doing these integrals analytically (actually, we used a symbolic integration package!), we package theresulting formulas in the following routine. Note that w(j+1 )returns F j(y;x). SUBROUTINE kermom(w,y,m) Returns in w(1:m) the first mindefinite-integral moments of one row of the singular part of the kernel. (For this example, mis hard-wired to be 4.) The input variable ylabels the column, while x(inCOMMON )i st h er o w . INTEGER m DOUBLE PRECISION w(m),y,x,d,df,clog,x2,x3,x4 COMMON /momcom/ x We can take xas the lower limit of integration. Thus, we return the moment integrals either purely to the left or purely to the right of the diagonal. if (y.ge.x) then d=y-xdf=2.d0*sqrt(d)*d w(1)=df/3.d0 18.3IntegralEquationswithSingularKernels 793Sample 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).w(2)=df*(x/3.d0+d/5.d0) w(3)=df*((x/3.d0 + 0.4d0*d)*x + d**2/7.d0) w(4)=df*(((x/3.d0 + 0.6d0*d)*x + 3.d0*d**2/7.d0)*x * + d**3/9.d0) else x2=x**2 x3=x2*xx4=x2*x2d=x-y clog=log(d) w(1)=d*(clog-1.d0)w(2)=-0.25d0*(3.d0*x+y-2.d0*clog*(x+y))*dw(3)=(-11.d0*x3+y*(6.d0*x2+y*(3.d0*x+2.d0*y)) * +6.d0*clog*(x3-y**3))/18.d0 w(4)=(-25.d0*x4+y*(12.d0*x3+y*(6.d0*x2+y* * (4.d0*x+3.d0*y)))+12.d0*clog*(x4-y**4))/48.d0 endif returnEND Next, we write a routine that constructs the quadrature matrix. SUBROUTINE quadmx(a,n,np)INTEGER n,np,NMAX REAL a(np,np),PI DOUBLE PRECISION xxPARAMETER (PI=3.14159265,NMAX=257)COMMON /momcom/ xx EXTERNAL kermom C USES wwghts,kermom Constructs in a(1:n,1:n) the quadrature matrix for an example Fredholm equation of the second kind. The nonsingular part of the kernel is computed within this routine, while the quadrature weights which integrate the singular part of the kernel are obtained via callsto wwghts . An external routine kermom , which supplies indefinite-integral moments of the singular part of the kernel, is passed to wwghts . INTEGER j,k REAL h,wt(NMAX),x,cx,yh=PI/(n-1) do 12j=1,n x=(j-1)*hxx=x Put xinCOMMON for use by kermom . call wwghts(wt,n,h,kermom) cx=cos(x) Part of nonsingular kernel. do 11k=1,n y=(k-1)*h a(j,k)=wt(k)*cx*cos(y) Put together all the pieces of the kernel. enddo 11 a(j,j)=a(j,j)+1. Since equation of the second kind, there is diagonal piece independent of h. enddo 12 return END Finally, we solve the linear system for any particular right-hand side, here sinx. PROGRAM fredex INTEGER NMAXREAL PI PARAMETER (NMAX=100,PI=3.14159265) INTEGER indx(NMAX),j,nREAL a(NMAX,NMAX),g(NMAX),x,d C USES quadmx,ludcmp,lubksb 794 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).3 2.52 1.51 .5 00.51 −.5f(x) xn = 10 n = 20 n = 40 Figure 18.3.1. Solution of the example integral equation (18.3.14) with grid sizes N=1 0,20, and 40. The tabulated solution values have been connected by straight lines; in practice one would interpolate a small Nsolution more smoothly. This sample program shows how to solve a Fredholm equation of the second kind using the product Nystrom method and a quadrature rule especially constructed for a particular, singular, kernel. n=40 Here the size of the grid is specified. call quadmx(a,n,NMAX) Make the quadrature matrix; all the action is here. call ludcmp(a,n,NMAX,indx,d) Decompose the matrix. do11j=1,n Construct the right hand side, here sin x. x=(j-1)*PI/(n-1)g(j)=sin(x) enddo 11 call lubksb(a,n,NMAX,indx,g) Backsubstitute. do12j=1,n Write out the solution. x=(j-1)*PI/(n-1) write (*,*) j,x,g(j) enddo 12 write (*,*) ’normal completion’ END With N=4 0, this program gives accuracy at about the 10−5level. The accuracy increases as N4(as it should for our Simpson-order quadrature scheme) despitethe highly singular kernel. Figure 18.3.1 shows the solution obtained, also plotting the solution forsmaller values of N, which are themselves seen to be remarkably faithful. Notice that the solution is smooth, even though the kernel is singular, a common occurrence. CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York). [1] Stroud, A.H., and Secrest, D. 1966, Gaussian Quadrature Formulas (Englewood Cliffs, NJ: Prentice-Hall). [2] Delves, L.M., and Mohamed, J.L. 1985, Computational Methods for Integral Equations (Cam- bridge, U.K.: Cambridge University Press). [3] 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.). [4] 18.4InverseProblemsandtheUse ofA PrioriInformation 795Sample 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.4 Inverse Problems and the Use of A Priori Information Later discussion will be facilitated by some preliminary mention of a couple of mathematical points. Suppose that uis an“unknown”vector that we plan to determine by some minimization principle. Let A[u]>0andB[u]>0be two positive functionals of u, so that we can try to determine uby either minimize: A[u]or minimize: B[u]( 18.4.1 ) (Ofcoursethese will generallygivedifferentanswers for u.) As anotherpossibility, nowsupposethat we want to minimize A[u]subject to the constraint thatB[u]have someparticularvalue,say b. ThemethodofLagrangemultipliersgivesthevariation δ δu{A[u]+λ1(B[u]−b)}=δ δu(A[u]+λ1B[u]) = 0 ( 18.4.2 ) where λ1is a Lagrange multiplier. Notice that bis absent in the second equality, since it doesn ’t depend on u. Next, suppose that we change our minds and decide to minimize B[u]subject to the constraint that A[u]have a particular value, a. Instead of equation (18.4.2) we have δ δu{B[u]+λ2(A[u]−a)}=δ δu(B[u]+λ2A[u]) = 0 ( 18.4.3 ) with, this time, λ2the Lagrange multiplier. Multiplying equation (18.4.3) by the constant 1/λ 2, and identifying 1/λ 2withλ1, we see that the actual variations are exactly the same in the two cases. Both cases will yield the same one-parameter family of solutions, say, u(λ1).A s λ1varies from 0to∞, the solution u(λ1) varies along a so-called trade-off curve between the problem of minimizing Aand the problem of minimizing B. Any solution along this curve can equally well be thought of as either (i) a minimization of Afor some constrained value of B, or (ii) a minimization of Bfor some constrained value of A, or (iii) a weighted minimization of the sum A+λ1B. Thesecondpreliminarypointhastodowith degenerate minimizationprinciples. In the example above, now suppose that A[u]has the particular form A[u]=|A·u−c|2(18.4.4 ) forsomematrix Aandvector c.I fAhasfewerrowsthancolumns,orif Ais square but degenerate (has a nontrivial nullspace, see §2.6, especially Figure 2.6.1), then minimizing A[u]willnotgive a unique solution for u. (To see why, review §15.4, and note that for a “design matrix ”Awith fewer rows than columns, the matrix AT·Ain the normal equations 15.4.10 is degenerate.) However, if we add any multiple λtimes a nondegeneratequadraticform B[u], forexample u·H·uwithH a positive de finite matrix, then minimization of A[u]+λB[u]willlead to a unique solution for u. (The sum of two quadratic forms is itself a quadratic form, with the second piece guaranteeing nondegeneracy.)