c18-3 (779610), страница 2
Текст из файла (страница 2)
Input to wwghts is a user-suppliedroutine, kermom, that is called to get the first four indefinite-integral moments of w(x), namelyZ yFm (y) ≡sm w(s)dsm = 0, 1, 2, 3(18.3.12)(The lower limit is arbitrary and can be chosen for convenience.) Cautionary note: Whencalled with N < 4, wwghts returns a rule of lower order than Simpson; you should structureyour problem to avoid this.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).h(18.3.6)[(k + 3)n+1 − kn+1 ]n+1The four terms in square brackets equation (18.3.5) each become independent of k, and(18.3.5) in fact reduces toZ (k+3)h3h9h9h3hf (x)dx =f (kh)+ f ([k +1]h)+ f ([k +2]h)+ f ([k +3]h) (18.3.7)8888khWn =800Chapter 18.Integral Equations and Inverse Theoryhh=h;hi=1.0/hh;for (j=1;j<=n;j++) wghts[j]=0.0;Zero all the weights so we can sum into them.(*kermom)(wold,0.0,4);Evaluate indefinite integrals at lower end.if (n >= 4) {Use highest available order.b=0.0;For another problem, you might changefor (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.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (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)hof 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.80118.3 Integral Equations with Singular KernelsWorked Example: A Diagonally Singular KernelAs a particular example, consider the integral equationZ πf (x) +K(x, y)f (y)dy = sin x(18.3.13)0K(x, y) = cos x cos y ×ln(x√ − y)y−xy<xy≥x(18.3.14)which has a logarithmic singularity on the left of the diagonal, combined with a square-rootdiscontinuity on the right.The first step is to do (analytically, in this case) the required moment integrals overthe singular part of the kernel, equation (18.3.12).
Since these integrals are done at a fixedvalue of x, we can use x as the lower limit. For any specified value of y, the requiredindefinite integral is then eitherZ yZ y−xFm (y; x) =sm (s − x)1/2 ds =(x + t)m t1/2 dtif y > x(18.3.15)xorZyFm (y; x) =xZ0x−ysm ln(x − s)ds =(x − t)m ln t dtif y < x(18.3.16)0(where a change of variable has been made in the second equality in each case). Doing theseintegrals analytically (actually, we used a symbolic integration package!), we package theresulting formulas in the following routine.
Note that w(j + 1) returns Fj (y; x).#include <math.h>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 ofthe 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;}}Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).with the (arbitrarily chosen) nasty kernel802Chapter 18.Integral Equations and Inverse TheoryNext, we write a routine that constructs the quadrature matrix.#include <math.h>#include "nrutil.h"#define PI 3.14159265double x;Communicates with kermom.wt=vector(1,n);h=PI/(n-1);for (j=1;j<=n;j++) {x=xx=(j-1)*h;Put x in global variable for use by kermom.wwghts(wt,n,h,kermom);cx=cos(xx);Part of nonsingular kernel.for (k=1;k<=n;k++) a[j][k]=wt[k]*cx*cos((k-1)*h);Put together all the pieces of the kernel.++a[j][j];Since equation of the second kind, there is diagonal piece}independent of h.free_vector(wt,1,n);}Finally, we solve the linear system for any particular right-hand side, here sin x.#include <stdio.h>#include <math.h>#include "nrutil.h"#define PI 3.14159265#define N 40Here 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 theproduct 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);Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).void quadmx(float **a, int n)Constructs in a[1..n][1..n] the quadrature matrix for an example Fredholm equation ofthe second kind. The nonsingular part of the kernel is computed within this routine, whilethe 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 thesingular 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;80318.3 Integral Equations with Singular Kernels10− .5n = 10n = 20n = 400.511.522.53xFigure 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 interpolatea 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 accuracyincreases as N 4 (as it should for our Simpson-order quadrature scheme) despite the highlysingular 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 thesolution 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 byDover Publications, New York). [1]Stroud, A.H., and Secrest, D. 1966, Gaussian Quadrature Formulas (Englewood Cliffs, NJ:Prentice-Hall).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
