3.1 Polynomial Interpolation and Extrapolation 109 The various P's form a“tableau'”with“ancestors'”on the left leading to a single "descendant"at the extreme right.For example,with N=4, E1: h=乃 P12 T2: 2=P2 P123 P23 P1234 (3.1.2) T3: 3=B P234 P34 E4: 4=P Neville's algorithm is a recursive way of filling in the numbers in the tableau 83 a column at a time,from left to right.It is based on the relationship between a “daughter'”P and its two“parents,” 容 P(+1).…(i+m)= (x-xi+m)Pi+1)…(i+m-1)+(zi-t)P+1)t+2.6+m) xi-工i+m RECIPES (3.1.3) This recurrence works because the two parents already agree at points +1... Ti+m-1. 二X分 Press. An improvement on the recurrence (3.1.3)is to keep track of the small differences between parents and daughters,namely to define (for m =1,2,..., 9 N-1)2 Progra Cm.i=Pi...(i+m)-Pi..(i+m-1) Dm,i≡P.(i+m)-P(i+1).(i+m小 (3.1.4) CIENTIFIC Then one can easily derive from(3.1.3)the relations Dm+li =itm+1-)(Cm.it1-Dm.) 工-工i+m+1 9 Cm+=)(Cmit1-Dm.) (3.1.5) Ti-工i+m+1 Numerical 10521 At each level m.the C's and D's are the corrections that make the interpolation one 43106 order higher.The final answer P1...N is equal to the sum of any yi plus a set of C's (outside Recipes and/or D's that form a path through the family tree to the rightmost daughter. Here is a routine for polynomial interpolation or extrapolation from N input points.Note that the input arrays are assumed to be unit-offset.If you have North Software. zero-offset arrays,remember to subtract 1 (see $1.2): #include <math.h> #include "nrutil.h" void polint(float xa[],float ya[],int n,float x,float y,float *dy) Given arrays xa[1..n]and ya[1..n],and given a value x,this routine returns a value y,and an error estimate dy.If P(x)is the polynomial of degree N-1 such that P(xai)=ya;,i= 1,...,n,then the returned value y P(x). int i,m,ns=1; float den,dif,dift,ho,hp,w;3.1 Polynomial Interpolation and Extrapolation 109 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). The various P’s form a “tableau” with “ancestors” on the left leading to a single “descendant” at the extreme right. For example, with N = 4, x1 : y1 = P1 P12 x2 : y2 = P2 P123 P23 P1234 x3 : y3 = P3 P234 P34 x4 : y4 = P4 (3.1.2) Neville’s algorithm is a recursive way of filling in the numbers in the tableau a column at a time, from left to right. It is based on the relationship between a “daughter” P and its two “parents,” Pi(i+1)...(i+m) = (x − xi+m)Pi(i+1)...(i+m−1) + (xi − x)P(i+1)(i+2)...(i+m) xi − xi+m (3.1.3) This recurrence works because the two parents already agree at points xi+1 ... xi+m−1. An improvement on the recurrence (3.1.3) is to keep track of the small differences between parents and daughters, namely to define (for m = 1, 2,..., N − 1), Cm,i ≡ Pi...(i+m) − Pi...(i+m−1) Dm,i ≡ Pi...(i+m) − P(i+1)...(i+m). (3.1.4) Then one can easily derive from (3.1.3) the relations Dm+1,i = (xi+m+1 − x)(Cm,i+1 − Dm,i) xi − xi+m+1 Cm+1,i = (xi − x)(Cm,i+1 − Dm,i) xi − xi+m+1 (3.1.5) At each level m, the C’s and D’s are the corrections that make the interpolation one order higher. The final answer P1...N is equal to the sum of any yi plus a set of C’s and/or D’s that form a path through the family tree to the rightmost daughter. Here is a routine for polynomial interpolation or extrapolation from N input points. Note that the input arrays are assumed to be unit-offset. If you have zero-offset arrays, remember to subtract 1 (see §1.2): #include <math.h> #include "nrutil.h" void polint(float xa[], float ya[], int n, float x, float *y, float *dy) Given arrays xa[1..n] and ya[1..n], and given a value x, this routine returns a value y, and an error estimate dy. If P(x) is the polynomial of degree N − 1 such that P(xai) = yai, i = 1,..., n, then the returned value y = P(x). { int i,m,ns=1; float den,dif,dift,ho,hp,w;