770 Chapter 17.Two Point Boundary Value Problems kmax[j]=km; err /nvars; fac=(err slowc slowc/err 1.0); Reduce correction applied when error is large for (j=1;j<=ne;j++) Apply corrections. jv=indexv[j]; for (k=k1;k<=k2;k++) y[j][k]-fac*c[jv][1][k]; printf("\n%8s %9s %9s\n","Iter.","Error","FAC"); Summary of corrections printf("%6d %12.6f %11.6f\n",it,err,fac); for this step. 常 if (err conv){ Point with largest error for each variable can free_vector(ermax,1,ne); be monitored by writing out kmax and free_ivector(kmax,1,ne); ermax. 83 return; 鱼 19881992 1-800 nrerror("Too many iterations in solvde"); Convergence failed. from NUMERICAL RECIPESI void bksub(int ne,int nb,int jf,int k1,int k2,float **c) (North Backsubstitution,used internally by solvde. 州bMe se int nbf,im,kp,k,j,i; America float xxi ART nbf=ne-nb; 1m=1; Progra for (k=k2;k>=k1;k--){ Use recurrence relations to eliminate remaining de- for their if (k ==k1)im=nbf+1; pendences. kp=k+1; for (j=1;j<=nbf;++){ xx=c[j][if][kp]; to dir for (i=im;i<mne;i++) c[i][jf]k]-=c[1][j]k]*xx; Numerical Recipes books or 18881920 OF SCIENTIFIC COMPUTING(ISBN for(k=k1;k<=k2;k++)( Reorder corrections to be in column 1. kp=k+1; for (i=1;i<=nb;i++)c[i][1][k]=c[i+nbf][if][k]; v@cambri for (i=1;i<=nbf;i++)c[i+nb][1][k]=c[i][jf][kp]; Numerical Recipes 10-621 43108 (outside #include <math.h> #include "nrutil.h" North Amer Software. void pinvs(int ie1,int ie2,int jel,int jsf,int jc1,int k,float ***c, f1oat米*s) Diagonalize the square subsection of the s matrix,and store the recursion coefficients in c; used internally by solvde. f int js1,jpiv,jp,je2,jcoff,j,irow,ipiv,id,icoff,i,*indxr; float pivinv,piv,dum,big,*pscl; indxr=ivector(ie1,ie2); pscl=vector(ie1,ie2); je2=je1+ie2-ie1; js1=je2+1;770 Chapter 17. Two Point Boundary Value Problems 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). kmax[j]=km; } err /= nvars; fac=(err > slowc ? slowc/err : 1.0); Reduce correction applied when error is large. for (j=1;j<=ne;j++) { Apply corrections. jv=indexv[j]; for (k=k1;k<=k2;k++) y[j][k] -= fac*c[jv][1][k]; } printf("\n%8s %9s %9s\n","Iter.","Error","FAC"); Summary of corrections printf("%6d %12.6f %11.6f\n",it,err,fac); for this step. if (err < conv) { Point with largest error for each variable can be monitored by writing out kmax and ermax. free_vector(ermax,1,ne); free_ivector(kmax,1,ne); return; } } nrerror("Too many iterations in solvde"); Convergence failed. } void bksub(int ne, int nb, int jf, int k1, int k2, float ***c) Backsubstitution, used internally by solvde. { int nbf,im,kp,k,j,i; float xx; nbf=ne-nb; im=1; for (k=k2;k>=k1;k--) { Use recurrence relations to eliminate remaining de￾if (k == k1) im=nbf+1; pendences. kp=k+1; for (j=1;j<=nbf;j++) { xx=c[j][jf][kp]; for (i=im;i<=ne;i++) c[i][jf][k] -= c[i][j][k]*xx; } } for (k=k1;k<=k2;k++) { Reorder corrections to be in column 1. kp=k+1; for (i=1;i<=nb;i++) c[i][1][k]=c[i+nbf][jf][k]; for (i=1;i<=nbf;i++) c[i+nb][1][k]=c[i][jf][kp]; } } #include <math.h> #include "nrutil.h" void pinvs(int ie1, int ie2, int je1, int jsf, int jc1, int k, float ***c, float **s) Diagonalize the square subsection of the s matrix, and store the recursion coefficients in c; used internally by solvde. { int js1,jpiv,jp,je2,jcoff,j,irow,ipiv,id,icoff,i,*indxr; float pivinv,piv,dum,big,*pscl; indxr=ivector(ie1,ie2); pscl=vector(ie1,ie2); je2=je1+ie2-ie1; js1=je2+1;