796 Chapter 18.Integral Equations and Inverse Theory void lubksb(float **a,int n,int *indx,float b[]); void ludcmp(float **a,int n,int *indx,float *d); int i,j,k,l,*indx; float d,sum,**a,*b; indx=ivector(1,m); a=matrix(1,m,1,m); b=vector(1,m); t[1]=t0; for(k=1;k<=m;k+)f[][1]=(*g)(k,t[1]); Initialize. for(1=2;1<=n;i+){ Take a step h. t[i]=t[i-1]+h; for(k=1;k<=m;k++)[ sum=(*g)(k,t[i]); Accumulate right-hand side of linear for(1=1;1<=m;1++){ equations in sum. 83g sum+=0.5*h*(*ak)(k,1,t[i],t[1])*f[1][1]; 18881892 for(j=2;j<1;j++) sum+=h*(*ak)(k,1,t[i],t[j])*f[1][j]; ak][1]=(k=1)-0.5*h*(*ak)(k,1,t[i],t[i]); Left-hand side goes in matrix a. b[k]=sum; from NUMERICAL RECIPESI ludcmp(a,m,indx,&d); Solve linear equations. lubksb(a,m,indx,b); for(k=1;k<=m;k++)f[k][i]=b[k]; THE free_vector(b,1,m); America computer, Press. free_matrix(a,1,m,1,m); free_ivector(indx,1,m); 9 For nonlinear Volterra equations,equation(18.2.4)holds with the product Kfi SCIENTIFIC replaced by Ki(fi),and similarly for the other two products of K's and f's.Thus for each i we solve a nonlinear equation for fi with a known right-hand side. Newton's method(89.4 or 89.6)with an initial guess of fi_1 usually works very well provided the stepsize is not too big Higher-order methods for solving Volterra equations are,in our opinion,not as 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,any method that achieves higher order by operating on several quadrature Numerica 10621 points simultaneously will need a special method to get started,when values at the 43106 first few points are not yet known. Second,stable quadrature rules can give rise to unexpected instabilities in Recipes integral equations.For example,suppose we try to replace the trapezoidal rule in the algorithm above with Simpson's rule.Simpson's rule naturally integrates over Software. an interval 2h,so we easily get the function values at the even mesh points.For the North odd mesh points,we could try appending one panel of trapezoidal rule.But to which end of the integration should we append it?We could do one step of trapezoidal rule followed by all Simpson's rule,or Simpson's rule with one step of trapezoidal rule at the end.Surprisingly,the former scheme 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 h and h/2.Then, assuming the error scales with h2,compute fe 4f(h/2)-f(h) (18.2.5) 3796 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). void lubksb(float **a, int n, int *indx, float b[]); void ludcmp(float **a, int n, int *indx, float *d); int i,j,k,l,*indx; float d,sum,**a,*b; indx=ivector(1,m); a=matrix(1,m,1,m); b=vector(1,m); t[1]=t0; for (k=1;k<=m;k++) f[k][1]=(*g)(k,t[1]); Initialize. for (i=2;i<=n;i++) { Take a step h. t[i]=t[i-1]+h; for (k=1;k<=m;k++) { sum=(*g)(k,t[i]); Accumulate right-hand side of linear for (l=1;l<=m;l++) { equations in sum. sum += 0.5*h*(*ak)(k,l,t[i],t[1])*f[l][1]; for (j=2;j<i;j++) sum += h*(*ak)(k,l,t[i],t[j])*f[l][j]; a[k][l]=(k == l)-0.5*h*(*ak)(k,l,t[i],t[i]); Left-hand side goes } in matrix a. b[k]=sum; } ludcmp(a,m,indx,&d); Solve linear equations. lubksb(a,m,indx,b); for (k=1;k<=m;k++) f[k][i]=b[k]; } free_vector(b,1,m); free_matrix(a,1,m,1,m); free_ivector(indx,1,m); } For nonlinear Volterra equations, equation (18.2.4) holds with the product K iifi replaced by Kii(fi), and similarly for the other two products of K’s and f’s. Thus for each i we solve a nonlinear equation for fi with a known right-hand side. Newton’s method (§9.4 or §9.6) with an initial guess of fi−1 usually works very well provided the stepsize is not too big. Higher-order methods for solving Volterra equations are, in our opinion, not as 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, any method that achieves higher order by operating on several quadrature points 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 the algorithm above with Simpson’s rule. Simpson’s rule naturally integrates over an interval 2h, so we easily get the function values at the even mesh points. For the odd mesh points, we could try appending one panel of trapezoidal rule. But to which end of the integration should we append it? We could do one step of trapezoidal rule followed by all Simpson’s rule, or Simpson’s rule with one step of trapezoidal rule at the end. Surprisingly, the former scheme 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 h and h/2. Then, assuming the error scales with h2, compute fE = 4f(h/2) − f(h) 3 (18.2.5)