18.3 Integral Equations with Singular Kernels 797 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 K or f.Variable stepsize methods are quite a bit more complicated than their counterparts for differential equations:we refer you to the literature [1,2]for a discussion. You should also be on the lookout for singularities in the integrand.If you find them.then look to 818.3 for additional ideas. CITED REFERENCES AND FURTHER READING: 8 Linz,P.1985,Analytical and Numerical Methods for Volterra Equations (Philadelphia:S.I.A.M.). [1] Delves,L.M.,and Mohamed,J.L.1985,Computationa/Methods for Integral Equations (Cam- bridge,U.K.:Cambridge University Press).[2] Cam ICAL RECIPES 18.3 Integral Equations with Singular Kernels Many integral equations have singularities in either the kernel or the solution or both. 巴互2。 A simple quadrature method will show poor convergence with N if such singularities are ignored.There is sometimes art in how singularities are best handled. We start with a few straightforward suggestions: 1.Integrable singularities can often be removed by a change of variable.For example,the singular behavior K(t,s)s1/2 or s-1/2 near s=0 can be removed by the transformation =s1/2.Note that we are assuming that the singular behavior is confined to K,whereas SCIENTIFIC the quadrature actually involves the product K(t,s)f(s),and it is this product that must be 6 "fixed."Ideally,you must deduce the singular nature of the product before you try a numerical solution,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)=6(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 a Gaussian quadrature based on w(s)as a weight function will work well.Even if the factorization is only approximate,the convergence is often improved dramatically.All you have to do is replace gauleg in the routine fred2 by another quadrature routine.Section Numerica 10621 4.5 explained how to construct such quadratures;or you can find tabulated abscissas and 431 weights in the standard references [1,2].You must of course supply K instead of K. This method is a special case of the product Nystrom method [3,4],where one factors out Recipes a singular term p(t,s)depending on both t and s from K and constructs suitable weights for (outside its Gaussian quadrature.The calculations in the general case are quite cumbersome,because the weights depend on the chosent as well as the form of p(t,s). Software. 首 We prefer to implement the product Nystrom method on a uniform grid,with a quadrature scheme that generalizes the extended Simpson's 3/8 rule (equation 4.1.5)to arbitrary weight functions.We discuss this in the subsections below. 3.Special quadrature formulas are also useful when the kernel is not strictly singular but is "almost"so.One example is when the kernel is concentrated near t =s on a scale much 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 ti.In such a scheme the narrow width of the kernel becomes an asset,rather than a liability:The quadrature becomes exact as the width of the kernel goes to zero. 4.An infinite range of integration is also a form of singularity.Truncating the range at a large finite value should be used only as a last resort.If the kernel goes rapidly to zero,then
18.3 Integral Equations with Singular Kernels 797 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 K or f. Variable stepsize methods are quite a bit more complicated than their counterparts for differential equations; we refer you to the literature [1,2] for a discussion. You should also be on the lookout for singularities in the integrand. If you find them, then look to §18.3 for additional ideas. CITED REFERENCES AND FURTHER READING: Linz, P. 1985, Analytical and Numerical Methods for Volterra Equations (Philadelphia: S.I.A.M.). [1] Delves, L.M., and Mohamed, J.L. 1985, Computational Methods for Integral Equations (Cambridge, 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 N if such singularities are ignored. There is sometimes art in how singularities are best handled. We start with a few straightforward suggestions: 1. Integrable singularities can often be removed by a change of variable. For example, the singular behavior K(t, s) ∼ s1/2 or s−1/2 near s = 0 can 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, you must deduce the singular nature of the product before you try a numerical solution, 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 a Gaussian quadrature based on w(s) as a weight function will work well. Even if the factorization is only approximate, the convergence is often improved dramatically. All you have to do is replace gauleg in the routine fred2 by another quadrature routine. Section 4.5 explained how to construct such quadratures; or you can find tabulated abscissas and weights in the standard references [1,2]. You must of course supply K instead of K. This method is a special case of the product Nystrom method [3,4], where one factors out a singular term p(t, s) depending on both t and s from K and constructs suitable weights for its Gaussian quadrature. The calculations in the general case are quite cumbersome, because the weights depend on the chosen {ti} as well as the form of p(t, s). We prefer to implement the product Nystrom method on a uniform grid, with a quadrature scheme that generalizes the extended Simpson’s 3/8 rule (equation 4.1.5) to arbitrary weight functions. We discuss this in the subsections below. 3. Special quadrature formulas are also useful when the kernel is not strictly singular, but is “almost” so. One example is when the kernel is concentrated near t = s on a scale much 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 ti. In such a scheme the narrow width of the kernel becomes an asset, rather than a liability: The quadrature becomes exact as the width of the kernel goes to zero. 4. An infinite range of integration is also a form of singularity. Truncating the range at a large finite value should be used only as a last resort. If the kernel goes rapidly to zero, then
798 Chapter 18.Integral Equations and Inverse Theory a Gauss-Laguerre [w~exp(-as)]or Gauss-Hermite [w~exp(-s2)]quadrature should work well.Long-tailed functions often succumb to the transformation 2+1-a (18.3.1) which maps 0>-1 so that Gauss-Legendre integration can be used. Here o>0 is 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.Here the Nystrom method fails completely because the kernel gets evaluated at(,s). Subtraction of the singularity is one possible cure: K(t,s)f(s)ds= K(t,s)[f(s)-f(t)]ds+ K(t,s)f(t)ds (18.3.2) K(t,s)[f(s)-f(t)]ds+r(t)f(t) 是2%a nted for where r(t)=K(t,s)ds is 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),ve get wj Kislfs -fi]+Arifi+gi (18.3.3) (Nort server 2 America University Press. THE Sometimes the subtraction process must be repeated before the kernel is completely regularized. See [3]for details.(And read on for a different,we think better,way to handle diagonal ART singularities. 9 Programs Quadrature on a Uniform Mesh with Arbitrary Weight It is possible in general to find n-point linear quadrature rules that approximate the 是 integral of a function f(),times an arbitrary weight function ()over an arbitrary range 可 of integration (a,b),as the sum of weights times n evenly spaced values of the function f(), say at kh,(+1)h,...,(k+n-1)h.The general scheme for deriving such quadrature rules is to write down then linear equations that must be satisfied if the quadrature rule is to be exact for the n functions f()=const,x,,...,x" and then solve these for the OF SCIENTIFIC COMPUTING(ISBN 0-521- coefficients.This can be done analytically,once and for all,if the moments of the weight function over the same range of integration, 1 (18.3.4) 1988-1992 by Numerical Recipes -43198.5 are assumed to be known.Here the prefactor h"is chosen to make Wn scale as h if(as in the usual case)b-a is proportional to h. (outside Carrying out this prescription for the four-point case gives the result Software. 。w(ofei= Amer 若hk++20k+3wo-32+12k+1)w+3k+2w-Wm +fk+1h)-kk+2k+3)Wo+(3k2+10k+6)W1-(3+5)W2+W +2fk+2h)k+1k+3)W%-3k2+8k+3)w+(3k+4W-W +a0k+3[-kk+1k+2w%+(3k2+6k+2)W-3+1w+W (18.3.5)
798 Chapter 18. Integral Equations and Inverse Theory Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). a Gauss-Laguerre [w ∼ exp(−αs)] or Gauss-Hermite [w ∼ exp(−s2)] quadrature should work well. Long-tailed functions often succumb to the transformation s = 2α z + 1 − α (18.3.1) which maps 0 z> −1 so that Gauss-Legendre integration can be used. Here α > 0 is 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. Here the Nystrom method fails completely because the kernel gets evaluated at (ti, si). Subtraction of the singularity is one possible cure: b a K(t, s)f(s) ds = b a K(t, s)[f(s) − f(t)] ds + b a K(t, s)f(t) ds = b a K(t, s)[f(s) − f(t)] ds + r(t)f(t) (18.3.2) where r(t) = b a K(t, s) ds is 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 fi = λ N j=1 j=i wjKij [fj − fi] + λrifi + gi (18.3.3) Sometimes the subtraction process must be repeated before the kernel is completely regularized. See [3] for details. (And read on for a different, we think better, way to handle diagonal singularities.) Quadrature on a Uniform Mesh with Arbitrary Weight 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 of integration (a, b), as the sum of weights times n evenly spaced values of the function f(x), say at x = kh,(k + 1)h, . . . ,(k + n − 1)h. The general scheme for deriving such quadrature rules is to write down the n linear equations that must be satisfied if the quadrature rule is to be exact for the n functions 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 weight function over the same range of integration, Wn ≡ 1 hn b a xnw(x)dx (18.3.4) are assumed to be known. Here the prefactor h−n is chosen to make Wn scale as h if (as in the usual case) b − a is proportional to h. Carrying out this prescription for the four-point case gives the result b a w(x)f(x)dx = 1 6 f(kh) (k + 1)(k + 2)(k + 3)W0 − (3k2 + 12k + 11)W1 + 3(k + 2)W2 − W3 + 1 2 f([k + 1]h) − k(k + 2)(k + 3)W0 + (3k2 + 10k + 6)W1 − (3k + 5)W2 + W3 + 1 2 f([k + 2]h) k(k + 1)(k + 3)W0 − (3k2 + 8k + 3)W1 + (3k + 4)W2 − W3 + 1 6 f([k + 3]h) − k(k + 1)(k + 2)W0 + (3k2 + 6k + 2)W1 − 3(k + 1)W2 + W3 (18.3.5)
18.3 Integral Equations with Singular Kernels 799 While the terms in brackets superficially appear to scale ask2,there is typically cancellation at both O()and O(). Equation (18.3.5)can be specialized to various choices of (a,b).The obvious choice is a =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 setting w(x)=1,in which case (18.3.4)becomes n+lk+3)m+1-kn+] (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 r(k+3)h kh e)k=3距， +警1+++2物+警r+A(83刀 81 Back to the case of general w(x),some other choices for a and b are also useful.For example,we may want to choose (a,b)to be (k+1h,+3h)or (k+2h,k+3h), 速 allowing us to finish off an extended rule whose number of intervals is not a multiple of three,without loss of accuracy:The integral will be estimated using the four values f(kh),...,f([+3]h).Even more useful is to choose (a,b)to be (+1]h,[+2]h),thus using four points to integrate a centered single interval.These weights,when sewed together RECIPES I into an extended formula,give quadrature schemes that have smooth coefficients,i.e.,without 令 the Simpson-like 2.4,2.4.2 alternation.(In fact,this was the technique that we used to derive 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 for f()a cubic polynomial.Rules of lower order,if desired,are similarly obtained.The Press. three point formula is u(f)=f(k+1(k+2Wo-(2k+3w+W +f(k+)-k(k+2)W+2(+1)W1-W (18.3.8) SCIENTIFIC 6 +[k+2kk+1Wo-(2k+1)W+W Here the simple special case is to take,w()=1,so that w=车7k+2+1-k h (18.3.9) Then equation (18.3.8)becomes Simpson's rule, Numerical 色 e)证=专+警rk+)+专k+2 (k+2)h (18.3.10) 43126 For nonconstant weight functions w(),however,equation (18.3.8)gives rules of one order (outside less than Simpson,since they do not benefit from the extra symmetry of the constant case. The two point formula is simply North (k+1)h w(x)f(x)da=f(kh)[(k+1)Wo-W+f(k+1]h)[-Wo+W](18.3.11) /kh Here is a routine wwghts that uses the above formulas to return an extended N-point quadrature rule for the interval (a,b)=(0,[N-1]h).Input to wwghts is a user-supplied routine,kermom,that is called to get the first four indefinite-integral moments of w(),namely 「y Fm(y）= s"w(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,wwghts returns a rule of lower order than Simpson;you should structure your problem to avoid this
18.3 Integral Equations with Singular Kernels 799 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). While the terms in brackets superficially appear to scale as k2, there is typically cancellation at both O(k2) and O(k). Equation (18.3.5) can be specialized to various choices of (a, b). The obvious choice is a = 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 setting w(x)=1, in which case (18.3.4) becomes Wn = 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 (k+3)h kh f(x)dx = 3h 8 f(kh)+ 9h 8 f([k+1]h)+ 9h 8 f([k+2]h)+ 3h 8 f([k+3]h) (18.3.7) Back to the case of general w(x), some other choices for a and b are 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 multiple of three, without loss of accuracy: The integral will be estimated using the four values f(kh),...,f([k + 3]h). Even more useful is to choose (a, b) to be ([k + 1]h, [k + 2]h), thus using four points to integrate a centered single interval. These weights, when sewed together into an extended formula, give quadrature schemes that have smooth coefficients, i.e., without the Simpson-like 2, 4, 2, 4, 2 alternation. (In fact, this was the technique that we used to derive 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 for f(x) a cubic polynomial. Rules of lower order, if desired, are similarly obtained. The three point formula is b a w(x)f(x)dx = 1 2 f(kh) (k + 1)(k + 2)W0 − (2k + 3)W1 + W2 + f([k + 1]h) − k(k + 2)W0 + 2(k + 1)W1 − W2 + 1 2 f([k + 2]h) k(k + 1)W0 − (2k + 1)W1 + W2 (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, (k+2)h kh f(x)dx = h 3 f(kh) + 4h 3 f([k + 1]h) + h 3 f([k + 2]h) (18.3.10) 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 (k+1)h kh w(x)f(x)dx = f(kh)[(k + 1)W0 − W1] + f([k + 1]h)[−kW0 + W1] (18.3.11) Here is a routine wwghts that uses the above formulas to return an extended N-point quadrature rule for the interval (a, b) = (0, [N − 1]h). Input to wwghts is a user-supplied routine, kermom, that is called to get the first four indefinite-integral moments of w(x), namely Fm(y) ≡ y 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, wwghts returns a rule of lower order than Simpson; you should structure your problem to avoid this
800 Chapter 18.Integral Equations and Inverse Theory void wwghts(float wghts[],int n,float h, void (*kermom)(double double ,int)) Constructs in wghts[1..n]weights for the n-point equal-interval quadrature from 0 to(n-1)h of a function f(x)times an arbitrary (possibly singular)weight function w()whose indefinite- integral moments Fn(y)are provided by the user-supplied routine kermom. int j,k; double wold[5],wnew[5],w[5],hh,hi,c,fac,a,b; Double precision on internal calculations even though the interface is in single precision. hh=h; hi=1.0/hh: for (i=1;i=4)[ Use highest available order. 19881992 b=0.0; For another problem,you might change for(j=1;j<n-3;j+)( this lower limit. c=j-1; This is called k in equation (18.3.5) a=b; Set upper and lower limits for this step. b=a+hh; FO州fro 1f(灯=n-3)b=(n-1)*hh Last interval:go all the way to end. (*kermom)(wnew,b,4); for (fac=1.0,k=1;k<=4;k++,fac*=hi) Equation (18.3.4). w[k]=(wnew [k]-wold[k])*fac; (North wghts[j]+ Equation (18.3.5). ((c+1.0)*(c+2.0)*(c+3.0)*w[1] America server computer, -(11.0+c*(12.0+c*3.0))*w[2] one paper University Press. THE ART +3.0*(c+2.0)*w[3]-8[4])/6.0); wghts[j+1]+ 是 (-c*(c+2.0)*(c+3.0)*w[1] Programs +(6.0+c*(10.0+c*3.0))*w[2] -(3.0*c+5.0)*w[3]+w[4])*0.5); ghts[+2]+=( st st (c*(c+1.0)*(c+3.0)*w[1] -(3.0+c*(8.0+c*3.0))*w[2] to dir +(3.0*c+4.0)w[3]-w[4])*0.5): wghts[j+3]+ (-c*(c+1.0)*(c+2.0)*w[1] OF SCIENTIFIC COMPUTING(ISBN +(2.0+c*(6.0+c*3.0))*w[2] -3.0*(c+1.0)*w[3]+8[4])/6.0); for (k=1;k<=4;k++)wold[k]=wnew[k]; Reset lower limits for moments. v@cam else if (n ==3){ Lower-order cases;not recommended. (*kermom)(wnew,hh+hh,3); w[1]=wnew[1]-wo1d[1]; .Further reproduction, 1988-1992 by Numerical Recipes 10-621 43108 w[2]=hi*(wnew [2]-wold[2]); w[3]=hi*hi*(wnew [3]-wold[3]); wghts[1]=w[1]-1.5*w[2]+0.5*w[3]; (outside wghts[2]=2.0*w[2]-w[3]; wghts[3]=0.5*(w[3]-w[2]): Software. ]e1se1f(n=2){ (*kermom)(wnew,hh,2); ying of wghts [1]=wnev[1]-wold[1]-(wghts [2]=hi*(wnew [2]-wold[2])); We will now give an example of how to apply wwghts to a singular integral equation
800 Chapter 18. Integral Equations and Inverse Theory Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). void wwghts(float wghts[], int n, float h, void (*kermom)(double [], double ,int)) Constructs in wghts[1..n] weights for the n-point equal-interval quadrature from 0 to (n−1)h of a function f(x) times an arbitrary (possibly singular) weight function w(x) whose indefiniteintegral moments Fn(y) are provided by the user-supplied routine kermom. { int j,k; double wold[5],wnew[5],w[5],hh,hi,c,fac,a,b; Double precision on internal calculations even though the interface is in single precision. hh=h; hi=1.0/hh; for (j=1;j= 4) { Use highest available order. b=0.0; For another problem, you might change for (j=1;j<=n-3;j++) { this lower limit. c=j-1; This is called k in equation (18.3.5). a=b; Set upper and lower limits for this step. b=a+hh; if (j == n-3) b=(n-1)*hh; Last interval: go all the way to end. (*kermom)(wnew,b,4); for (fac=1.0,k=1;k<=4;k++,fac*=hi) Equation (18.3.4). w[k]=(wnew[k]-wold[k])*fac; wghts[j] += ( Equation (18.3.5). ((c+1.0)*(c+2.0)*(c+3.0)*w[1] -(11.0+c*(12.0+c*3.0))*w[2] +3.0*(c+2.0)*w[3]-w[4])/6.0); wghts[j+1] += ( (-c*(c+2.0)*(c+3.0)*w[1] +(6.0+c*(10.0+c*3.0))*w[2] -(3.0*c+5.0)*w[3]+w[4])*0.5); wghts[j+2] += ( (c*(c+1.0)*(c+3.0)*w[1] -(3.0+c*(8.0+c*3.0))*w[2] +(3.0*c+4.0)*w[3]-w[4])*0.5); wghts[j+3] += ( (-c*(c+1.0)*(c+2.0)*w[1] +(2.0+c*(6.0+c*3.0))*w[2] -3.0*(c+1.0)*w[3]+w[4])/6.0); for (k=1;k<=4;k++) wold[k]=wnew[k]; Reset lower limits for moments. } } else if (n == 3) { Lower-order cases; not recommended. (*kermom)(wnew,hh+hh,3); w[1]=wnew[1]-wold[1]; w[2]=hi*(wnew[2]-wold[2]); w[3]=hi*hi*(wnew[3]-wold[3]); wghts[1]=w[1]-1.5*w[2]+0.5*w[3]; wghts[2]=2.0*w[2]-w[3]; wghts[3]=0.5*(w[3]-w[2]); } else if (n == 2) { (*kermom)(wnew,hh,2); wghts[1]=wnew[1]-wold[1]-(wghts[2]=hi*(wnew[2]-wold[2])); } } We will now give an example of how to apply wwghts to a singular integral equation
18.3 Integral Equations with Singular Kernels 801 Worked Example:A Diagonally Singular Kernel As a particular example,consider the integral equation f(x)+ K(x,y)f(y)dy sin x (18.3.13) with the (arbitrarily chosen)nasty kernel K(x,）=cosxcosy× ∫-ln(e-)yx (18.3.15) from NUMERICAL RECIPES I 18881892 (Nort server 令 Fm(y;x）=- s"M In(x-s)ds (x-t)"Intdt ify extern double x; Defined in quadmx. 兰 void kermom(double w,double y,int m) Returns in w[1..m]the first m indefinite-integral moments of one row of the singular part of the kernel.(For this example,m is hard-wired to be 4.)The input variable y labels the column, while the global variable x is the row.We can take x as the lower limit of integration.Thus, we return the moment integrals either purely to the left or purely to the right of the diagonal. double d,df,clog,x2,x3,x4,y2; if (y >=x){ d=y-x; df=2.0*sgrt(d)*d; 和6 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING (ISBN 0-521 w[1]=df/3.0; Further reproduction,or -431085 w[2]=df*(x/3.0+d/5.0); w[3]=df*((x/3.0+0.4*d)*x+d*d/7.0): w[4]=df*(((x/3.0+0.6*d)*x+3.0*d*d/7.0)*x+d*d*d/9.0); else x3=(x2=x*x)*x; X4=X2*X2; y2=y*y; (outside North America). Software. d=x-yi w[1]=d*((c1og=1og(d))-1.0); w[2]=-0.25*(3.0*x+y-2.0*c1og*(x+y)*d; w[3]=(-11.0*x3+y*(6.0*x2+y*(3.0*x+2.0*y)) +6.0*c1og*(x3-y*y2))/18.0; w[4]=(-25.0*x4+y*(12.0*x3+y*(6.0*x2+y* (4.0*x+3.0*y))+12.0*c1og*(x4-(y2*y2))/48.0; 2
18.3 Integral Equations with Singular Kernels 801 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Worked Example: A Diagonally Singular Kernel As a particular example, consider the integral equation f(x) + π 0 K(x, y)f(y)dy = sin x (18.3.13) with the (arbitrarily chosen) nasty kernel K(x, y) = cos x cos y × − √ ln(x − y) yx (18.3.15) or Fm(y; x) = − y x sm ln(x − s)ds = x−y 0 (x − t) m ln t dt if y extern double x; Defined in quadmx. void kermom(double w[], double y, int m) Returns in w[1..m] the first m indefinite-integral moments of one row of the singular part of the kernel. (For this example, m is hard-wired to be 4.) The input variable y labels the column, while the global variable x is the row. We can take x as the lower limit of integration. Thus, we return the moment integrals either purely to the left or purely to the right of the diagonal. { double d,df,clog,x2,x3,x4,y2; if (y >= x) { d=y-x; df=2.0*sqrt(d)*d; w[1]=df/3.0; w[2]=df*(x/3.0+d/5.0); w[3]=df*((x/3.0 + 0.4*d)*x + d*d/7.0); w[4]=df*(((x/3.0 + 0.6*d)*x + 3.0*d*d/7.0)*x+d*d*d/9.0); } else { x3=(x2=x*x)*x; x4=x2*x2; y2=y*y; d=x-y; w[1]=d*((clog=log(d))-1.0); w[2] = -0.25*(3.0*x+y-2.0*clog*(x+y))*d; w[3]=(-11.0*x3+y*(6.0*x2+y*(3.0*x+2.0*y)) +6.0*clog*(x3-y*y2))/18.0; w[4]=(-25.0*x4+y*(12.0*x3+y*(6.0*x2+y* (4.0*x+3.0*y)))+12.0*clog*(x4-(y2*y2)))/48.0; } }
802 Chapter 18. Integral Equations and Inverse Theory Next,we write a routine that constructs the quadrature matrix. #include #include "nrutil.h" #def1nePI3.14159265 double x; Communicates with kermom void quadmx(float *a,int n) 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 http://www. the quadrature weights which integrate the singular part of the kernel are obtained via calls 常 to wwghts.An external routine kermom,which supplies indefinite-integral moments of the singular part of the kernel,is passed to wwghts void kermom(double w[],double y,int m); granted for 18881992 void wwghts(float wghts[],int n,float h, void (*kermom)(double []double ,int)); 100 (including this one) int j,ki float h,*wt,xx,cx; Cambridge from NUMERICAL RECIPES I wt=vector(1,n); h=PI/(n-1); for (j=1;j #include ectcustser OF SCIENTIFIC COMPUTING(ISBN #include "nrutil.h" 18881920 #define PI 3.14159265 #define N 40 Here the size of the grid is specified. int main(void)/*Program fredex * 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, Numerical Recipes 10-621 -43108 kernel. void lubksb(float **a,int n,int *indx,float b[]); (outside void ludcmp(float **a,int n,int *indx,float *d); void quadmx(float **a,int n); North Software. float **a,d,*g,x; int *indx,ji Amer indx=ivector(1,N); a=matrix(1,N,1,N); g=vector(1,N); quadmx(a,N) Make the quadrature matrix;all the action is here. ludcmp(a,N,indx,&d); Decompose the matrix. for(j=1;j<=N;j++)g[j]=s1n((j-1)*PI/(W-1)); Construct the right hand side,here sinz. lubksb(a,N,indx,g); Backsubstitute. for (j=1;j<=N;j++) Write out the solution. x=(j-1)*PI/(W-1);
802 Chapter 18. Integral Equations and Inverse Theory Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Next, we write a routine that constructs the quadrature matrix. #include #include "nrutil.h" #define PI 3.14159265 double x; Communicates with kermom. void quadmx(float **a, int n) 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 calls to wwghts. An external routine kermom, which supplies indefinite-integral moments of the singular part of the kernel, is passed to wwghts. { void kermom(double w[], double y, int m); void wwghts(float wghts[], int n, float h, void (*kermom)(double [], double ,int)); int j,k; float h,*wt,xx,cx; wt=vector(1,n); h=PI/(n-1); for (j=1;j #include #include "nrutil.h" #define PI 3.14159265 #define N 40 Here the size of the grid is specified. int main(void) /* Program fredex */ 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. { void lubksb(float **a, int n, int *indx, float b[]); void ludcmp(float **a, int n, int *indx, float *d); void quadmx(float **a, int n); float **a,d,*g,x; int *indx,j; indx=ivector(1,N); a=matrix(1,N,1,N); g=vector(1,N); quadmx(a,N); Make the quadrature matrix; all the action is here. ludcmp(a,N,indx,&d); Decompose the matrix. for (j=1;j<=N;j++) g[j]=sin((j-1)*PI/(N-1)); Construct the right hand side, here sin x. lubksb(a,N,indx,g); Backsubstitute. for (j=1;j<=N;j++) { Write out the solution. x=(j-1)*PI/(N-1);
18.3 Integral Equations with Singular Kernels 803 TTTT http://www.nr Permission is read able files …n=10 -----n=20 -n=40 .com or call (including this one) granted fori 19881992 11-800 -.5 LLLL1LLLLLLLLLLLLLLLLLLLLLLLI 0 .5 1 1.5 2 2.5 3 Figure 18.3.1.Solution of the example integral equation(18.3.14)with grid sizes N=10,20,and 40. The tabulated solution values have been connected by straight lines;in practice one would interpolate to any server from NUMERICAL RECIPES IN C: a small N solution more smoothly. (North America Cambridge University Press. tusers to make one paper computer, THE printf("%6.2d%12.6f%12.6f\n",j,x,g[j]); 是 ART free_vector(g,1,N); free_matrix(a,1,N,1,N); st st Programs free_ivector(indx,1,N); return 0; With N=40,this program gives accuracy at about the 10-5 level.The accuracy to dir increases as N(as it should for our Simpson-order quadrature scheme)despite the highly singular kernel.Figure 18.3.1 shows the solution obtained,also plotting the solution for smaller values of N,which are themselves seen to be remarkably faithful.Notice that the 1788-1982 OF SCIENTIFIC COMPUTING(ISBN solution is smooth,even though the kernel is singular,a common occurrence. Further reproduction Numerical Recipes 10-621 CITED REFERENCES AND FURTHER READING: 43108 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 (outside Dover Publications,New York).[1] Stroud,A.H.,and Secrest,D.1966,Gaussian Quadrature Formulas (Englewood Cliffs,NJ: North Software. Prentice-Hall).[2] 6 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.3 Integral Equations with Singular Kernels 803 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 0 .5 1 1.5 2 2.5 3 0 .5 1 −.5 f(x) x n = 10 n = 20 n = 40 Figure 18.3.1. Solution of the example integral equation (18.3.14) with grid sizes N = 10, 20, and 40. The tabulated solution values have been connected by straight lines; in practice one would interpolate a small N solution more smoothly. printf("%6.2d %12.6f %12.6f\n",j,x,g[j]); } free_vector(g,1,N); free_matrix(a,1,N,1,N); free_ivector(indx,1,N); return 0; } With N = 40, this program gives accuracy at about the 10−5 level. The accuracy increases as N4 (as it should for our Simpson-order quadrature scheme) despite the highly singular kernel. Figure 18.3.1 shows the solution obtained, also plotting the solution for smaller 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 Mathematics 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 (Cambridge, 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]
804 Chapter 18.Integral Equations and Inverse Theory 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 u is an "unknown"vector that we plan to determine by some minimization principle.Let A[u]>0 and Blu]>0 be two positive functionals of u,so that we can try to determine u by either minimize:Au or minimize:B[u] (18.4.1) 8 (Of course these will generally give different answers for u.)As another possibility, now suppose that we want to minimize Au]subject to the constraint that Bu]have some particular value,say b.The method of Lagrange multipliers gives the variation 0 {4[u+X1(Bu-b)}=。(4[u+1B[▣)=0 (18.4.2) RECIPES where A is a Lagrange multiplier.Notice that b is absent in the second equality, 9 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 Alu]have a particular value,a.Instead of equation(18.4.2) we have 、豆之w 9 9 6 0 {B网+(A回-a}=a(B[回+4)=0 (18.4.3) with,this time,A2 the Lagrange multiplier.Multiplying equation(18.4.3)by the 61 constant 1/A2,and identifying 1/A2 with A1,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).As A1 varies from 0 to oo,the solution u(1) varies along a so-called trade-off curve between the problem of minimizing A and the problem of minimizing B.Any solution along this curve can equally well be thought of as either (i)a minimization ofA for some constrained value of B, or(ii)a minimization of B for some constrained value ofA,or(iii)a weighted Numerica 10621 minimization of the sum A+AB. 43106 The second preliminary point has to do with degenerate minimization principles. In the example above,now suppose that Au]has the particular form A[回=A·u-c2 (18.4.4) for some matrix A and vector c.If A has fewer rows than columns,or if A is square but degenerate (has a nontrivial nullspace,see $2.6,especially Figure 2.6.1),then minimizing A[u]will not give a unique solution for u.(To see why,review 815.4, and note that for a"design matrix"A with fewer rows than columns,the matrix A.A in the normal equations 15.4.10 is degenerate.)However,if we add any multiple A times a nondegenerate quadratic form B[u],for example u.H.u with H a positive definite matrix,then minimization of Au+ABu]will lead to a unique solution for u.(The sum of two quadratic forms is itself a quadratic form,with the second piece guaranteeing nondegeneracy.)
804 Chapter 18. Integral Equations and Inverse Theory Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 u is an “unknown” vector that we plan to determine by some minimization principle. Let A[u] > 0 and B[u] > 0 be two positive functionals of u, so that we can try to determine u by either minimize: A[u] or minimize: B[u] (18.4.1) (Of course these will generally give different answers for u.) As another possibility, now suppose that we want to minimize A[u] subject to the constraint that B[u] have some particular value, say b. The method of Lagrange multipliers gives the variation δ δu {A[u] + λ1(B[u] − b)} = δ δu (A[u] + λ1B[u]) = 0 (18.4.2) where λ1 is a Lagrange multiplier. Notice that b is 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, λ2 the Lagrange multiplier. Multiplying equation (18.4.3) by the constant 1/λ2, and identifying 1/λ2 with λ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). As λ1 varies from 0 to ∞, the solution u(λ1) varies along a so-called trade-off curve between the problem of minimizing A and the problem of minimizing B. Any solution along this curve can equally well be thought of as either (i) a minimization of A for some constrained value of B, or (ii) a minimization of B for some constrained value of A, or (iii) a weighted minimization of the sum A + λ1B. The second preliminary point has to do with degenerateminimization principles. In the example above, now suppose that A[u] has the particular form A[u] = |A · u − c| 2 (18.4.4) for some matrix A and vector c. If A has fewer rows than columns, or if A is square but degenerate (has a nontrivial nullspace, see §2.6, especially Figure 2.6.1), then minimizing A[u] will not give a unique solution for u. (To see why, review §15.4, and note that for a “design matrix” A with fewer rows than columns, the matrix AT · A in the normal equations 15.4.10 is degenerate.) However, if we add any multiple λ times a nondegenerate quadratic form B[u], for example u · H · u with H a positive definite matrix, then minimization of A[u] + λB[u] will lead to a unique solution for u. (The sum of two quadratic forms is itself a quadratic form, with the second piece guaranteeing nondegeneracy.)