3.5 Coefficients of the Interpolating Polynomial 121 #include "nrutil.h" void polcoe(float x[],float y[],int n,float cof[]) Given arrays x[0..n]and y[0..n]containing a tabulated function y:=f(xi),this routine returns an array of coefficients cof [0..n],such that y;=,cofj int k,j,i; float phi,ff,b,*s; s=vector(0,n); for(i=0;1<=n;i+)s[i]=cof[i]=0.0; s[n]=-x[o]; for(i=1;i<=n;i++)[ Coefficients si of the master polynomial P(z)are for (j=n-i;j<=n-1;j++) found by recurrence. s[j]-=x[1]*s[j+1]; s[n]-=x[i]; 18881992 for(j-=0;j<=n;j+){ phisn+1; for (k=n;k>=1;k--) The quantity phi=j(j-k)is found as a phi=k*s[k]+x[j]*phi; derivative of P(j). ff=y[j]/phi; RECIPES I b=1.0: Coefficients of polynomials in each term of the La- for (k=n;k>=0;k--){ grange formula are found by synthetic division of cof[k]+b*ff; P(z)by (x-zj).The solution ck is accumu- b=s[k]+x[j]*b; lated. free_vector(s,0,n); Another Method to dir Another technique is to make use of the function value interpolation routine already given (polint $3.1).If we interpolate (or extrapolate)to find the value of SCIENTIFIC COMPUTING (ISBN the interpolating polynomial at z=0,then this value will evidently be co.Now 19841820 we can subtract co from the yi's and divide each by its corresponding zi.Throwing out one point (the one with smallest is a good candidate),we can repeat the procedure to find c1,and so on. Numerical 10-521 It is not instantly obvious that this procedure is stable,but we have generally Recipes 43106 found it to be somewhat more stable than the routine immediately preceding.This method is of order N3,while the preceding one was of order N2.You will find,however,that neither works very well for large N,because of the intrinsic ill-condition of the Vandermonde problem.In single precision,N up to 8 or 10 is North Software. satisfactory:about double this in double precision. #include <math.h> #include "nrutil.h" void polcof(float xa[],float ya[],int n,float cof[]) Given arrays xa[0..n]and ya[0..n]containing a tabulated function ya;=f(xai).this routine returns an array of coefficients cof [..n]such that ya;=cofjxa. f void polint(float xa[],float ya[],int n,float x,float *y,float *dy); int k,j,i; float xmin,dy,*x,*y;3.5 Coefficients of the Interpolating Polynomial 121 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). #include "nrutil.h" void polcoe(float x[], float y[], int n, float cof[]) Given arrays x[0..n] and y[0..n] containing a tabulated function yi = f(xi), this routine returns an array of coefficients cof[0..n], such that yi =  j cofjxj i . { int k,j,i; float phi,ff,b,*s; s=vector(0,n); for (i=0;i<=n;i++) s[i]=cof[i]=0.0; s[n] = -x[0]; for (i=1;i<=n;i++) { Coefficients si of the master polynomial P(x) are for (j=n-i;j<=n-1;j++) found by recurrence. s[j] -= x[i]*s[j+1]; s[n] -= x[i]; } for (j=0;j<=n;j++) { phi=n+1; for (k=n;k>=1;k--) The quantity phi =  j
=k(xj − xk) is found as a phi=k*s[k]+x[j]*phi; derivative of P(xj ). ff=y[j]/phi; b=1.0; Coefficients of polynomials in each term of the La￾grange formula are found by synthetic division of P(x) by (x − xj ). The solution ck is accumu￾lated. for (k=n;k>=0;k--) { cof[k] += b*ff; b=s[k]+x[j]*b; } } free_vector(s,0,n); } Another Method Another technique is to make use of the function value interpolation routine already given (polint §3.1). If we interpolate (or extrapolate) to find the value of the interpolating polynomial at x = 0, then this value will evidently be c 0. Now we can subtract c0 from the yi’s and divide each by its corresponding xi. Throwing out one point (the one with smallest xi is a good candidate), we can repeat the procedure to find c1, and so on. It is not instantly obvious that this procedure is stable, but we have generally found it to be somewhat more stable than the routine immediately preceding. This method is of order N 3, while the preceding one was of order N 2. You will find, however, that neither works very well for large N, because of the intrinsic ill-condition of the Vandermonde problem. In single precision, N up to 8 or 10 is satisfactory; about double this in double precision. #include <math.h> #include "nrutil.h" void polcof(float xa[], float ya[], int n, float cof[]) Given arrays xa[0..n] and ya[0..n] containing a tabulated function yai = f(xai), this routine returns an array of coefficients cof[0..n] such that yai =  j cofjxaj i . { void polint(float xa[], float ya[], int n, float x, float *y, float *dy); int k,j,i; float xmin,dy,*x,*y;