802 Chapter 18. Integral Equations and Inverse Theory Next,we write a routine that constructs the quadrature matrix. #include <math.h> #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<=n;j++){ server computer, to make x=xx=(-1)*h; Put x in global variable for use by kermom wwghts(wt,n,h,kermom); UnN电.t THE CX=cOs(xx); Part of nonsingular kernel. (North America one paper for (k=1;k<=n;k++)a[j][k]=wt[k]*cx*cos((k-1)*h); ART Put together all the pieces of the kernel. ++a[j][j]; Since equation of the second kind,there is diagonal piece independent of h. st st Programs free_vector(wt,1,n); Finally,we solve the linear system for any particular right-hand side,here sin. #include <stdio.h> #include <math.h> 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 machine￾readable 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 <math.h> #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<=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 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);