Lecture Note 3:Polynomial Interpolation Xiaoqun Zhang Shanghai Jiao Tong University Last updated:October 16,2012
Lecture Note 3: Polynomial Interpolation Xiaoqun Zhang Shanghai Jiao Tong University Last updated: October 16, 2012
Lecture note 3 Numerical Analysis 1.1 Introduction We first look at some examples. Lookup table for f(r)=e-dr f(0.520)=0.53790,f(0.521)=0.53876 f(0.522)=0.53962,f(0.523)=0.54048 No f(0.52136)or f(0.52218).How should we do? Census of population: 1940132165 1950151326 1960179323 19702033021980226542 1960249633 There is no record for 1965.How to predict it? General questions:Let (xo,yo),(1,),...,(n,Un)be given.Given a x which is not in fro,x1,...,n},what is the corresponding y? Find a function f(r)such that yo f(o),y1=f(r1),...,Un=f(n). It's hard.There are too "many"functions.One solution is to use polynomials since polynomials have nice properties. Polynomial Interpolation:find a polynomial p(x)such that p(xo)=0,p(r1)=1,,p(rn)=yn 1.1.1 Existence and Uniqueness of polynomial interpolation Questions on polynomial interpolation: Does such a polynomial exist? -Generally,you can not find a linear function to interpolate 3 points. If yes,what is the degree? -2 points,linear;3 points,quadratic;4 points cubic;.. ·Is it unique? -2 points,linear;3 points,quadratic;4 points cubic;... What is a formula for producing p(z)from the given coordinates? -2 points,linear;3 points,quadratic;4 points cubic;... (x0,0):·,(xn,n) 2
Lecture note 3 Numerical Analysis 1.1 Introduction We first look at some examples. • Lookup table for f(x) = √ 2 π R x 0 e −x 2 dx f(0.520) = 0.53790, f(0.521) = 0.53876 f(0.522) = 0.53962, f(0.523) = 0.54048 No f(0.52136) or f(0.52218). How should we do? • Census of population: 1940 132165 1950 151326 1960 179323 1970 203302 1980 226542 1960 249633 There is no record for 1965. How to predict it? General questions: Let (x0, y0), (x1, y1), . . . , (xn, yn) be given. Given a x which is not in {x0, x1, . . . , xn}, what is the corresponding y? • Find a function f(x) such that y0 = f(x0), y1 = f(x1), . . . , yn = f(xn). • It’s hard. There are too “many” functions. One solution is to use polynomials since polynomials have nice properties. • Polynomial Interpolation: find a polynomial p(x) such that p(x0) = y0, p(x1) = y1, . . . , p(xn) = yn 1.1.1 Existence and Uniqueness of polynomial interpolation Questions on polynomial interpolation: • Does such a polynomial exist? – Generally, you can not find a linear function to interpolate 3 points. • If yes, what is the degree? – 2 points, linear; 3 points, quadratic; 4 points cubic;... • Is it unique? – 2 points, linear; 3 points, quadratic; 4 points cubic;... • What is a formula for producing p(x) from the given coordinates? – 2 points, linear; 3 points, quadratic; 4 points cubic;... • (x0, y0), . . . ,(xn, yn). 2
Lecture note 3 Numerical Analysis ● p()=ao ai a2x2+...+amz" .We want to find ao,a1,...,am. m+1 parameters,n+1 conditions (p(zi)=yi) It is reasonable to try first with n=m. So we want to find ao,a1,...,am such that ao a1zo a2o+...+anto y0 do a1z1+a2i+...+anzn=y ao alIn a2zn+...+antn yn Write it into a system of linear equations n To 哈 昭7 ao 「3o1 1 1 “ 理 1 T2 吃 02 .. . . 1 In 品 zn] an Ln」 Xa=可 X-Vandermonde matrix. ·“There exists a unique vector such that Xa=” is equivalent to“det(X)≠0.” Theorem 1 (Theorem 1.1)Given n+1 distinct points xo,11,...,In and the corresponding yo,1,...,yn:there is a polynomial p(x)of degree <n such that p(ro）=o,p(r）=h1,,p(xn)=yn·This polynomial is unique among the set of all polynomials of degree <n. Why degree n and not n? -Why this set? Proof.We will prove that det(X)0.Let To 哈 昭 1 T1 1 T2 Vn(z)= 1 n-1品1 In-1 3
Lecture note 3 Numerical Analysis • p(x) = a0 + a1x + a2x 2 + . . . + amx m • We want to find a0, a1, . . . , am. m + 1 parameters, n + 1 conditions (p(xi) = yi) It is reasonable to try first with n = m. • So we want to find a0, a1, . . . , am such that a0 + a1x0 + a2x 2 0 + . . . + anx n 0 = y0 a0 + a1x1 + a2x 2 1 + . . . + anx n 1 = y1 . . . a0 + a1xn + a2x 2 n + . . . + anx n n = yn • Write it into a system of linear equations 1 x0 x 2 0 . . . xn 0 1 x1 x 2 1 . . . xn 1 1 x2 x 2 2 . . . xn 2 . . . . . . . . . 1 xn x 2 n . . . xn n a0 a1 a2 . . . an = y0 y1 y2 . . . yn X~a = ~y • X — Vandermonde matrix. • “There exists a unique vector such that X~a = ~y” is equivalent to “det(X) 6= 0.” Theorem 1 (Theorem 1.1) Given n + 1 distinct points x0, x1, . . . , xn and the corresponding y0, y1, . . . , yn, there is a polynomial p(x) of degree ≤ n such that p(x0) = y0, p(x1) = y1, . . . , p(xn) = yn. This polynomial is unique among the set of all polynomials of degree ≤ n. – Why degree ≤ n and not n? – Why this set? Proof. We will prove that det(X) 6= 0. Let Vn(x) = 1 x0 x 2 0 . . . xn 0 1 x1 x 2 1 . . . xn 1 1 x2 x 2 2 . . . xn 2 . . . . . . . . . 1 xn−1 x 2 n−1 . . . xn n−1 1 x x2 . . . xn 3
Lecture note 3 Numerical Analysis so that X=Vn(In). We have the claims: 1.Vn(r)is a polynomial of degree n. 2.The roots are zo,1,...,In-1. 3.The coefficients of the term z"is Vn-1(n-1). Proof. 1.Expand the last row by minor. 2.If two row are the same then det(X)=0. 3.Expand the last row by minor. ◇ So Vn(z)=(2n-1)(x-z0)(-z1)...(x-In-1)Vn-1. By induction, ⅓a=dtl到=a-到 V2(x2)=(x2-t0)(x2-x1)(r1-x0) (x3)=(x3-to)(x3-x1)(a3-x2).. Vn(xn)=Πl0≤i<jsn(rj-xi. Since zo,...,n are distinct,det(X)=Vn(Xn)0. 0 1.2 Lagrange Formula Find a polynomial P(z)of order <n passes through (r0,0),(x1,班)，(x2,2),.(rn,h) Method of undetermined coefficients: Pn(z）=a0+a1x+a2x2+..+anx”. Need to solve Xa=灭 Computational cost is high. 4
Lecture note 3 Numerical Analysis so that X = Vn(xn). We have the claims: 1. Vn(x) is a polynomial of degree n. 2. The roots are x0, x1, . . . , xn−1. 3. The coefficients of the term x n is Vn−1(xn−1). Proof. 1. Expand the last row by minor. 2. If two row are the same then det(X) = 0. 3. Expand the last row by minor. So Vn(x) = (xn−1)(x − x0)(x − x1). . .(x − xn−1)Vn−1. By induction, V1(x1) = det 1 x0 1 x1 = (x1 − x0) V2(x2) = (x2 − x0)(x2 − x1)(x1 − x0) V3(x3) = (x3 − x0)(x3 − x1)(x3 − x2). . . Vn(xn) = Π0≤i<j≤n(xj − xi). Since x0, . . . , xn are distinct, det(X) = Vn(Xn) 6= 0. 1.2 Lagrange Formula • Find a polynomial Pn(x) of order ≤ n passes through (x0, y0), (x1, y1), (x2, y2), . . .(xn, yn). • Method of undetermined coefficients: Pn(x) = a0 + a1x + a2x 2 + . . . + anx n . Need to solve X~a = ~y. Computational cost is high. 4
Lecture note 3 Numerical Analysis Computational cost is high: Gaussian Elimination O(n3). Find fast solvers of the linear equation or explicit expression of a. 。Use“bases'”other than{L,x,z2,,z}. .Two points case:(To,y0),(x1,21). Construct two linear functions (polynomials of order 1)L1.o(z)and L11(x) such that: L1.0(z0)=1L1.0(x1)=0 and L1,1(xo)=0L1,1(x1)=1. Then,let 乃(x)=0L1,o(x)+hL1,1(x) We have P1(xo)=0L1,0(xo)+y1L1,1(x0)=0 (c1)=0L1,0(x1)+1L1,1(x1)=h. ·Construct L1,o(r): L1.0(a)=工- T0-E1 Use a graph to illustrate it!. Similarly, L1,1()=-0 x1一T0 Generalize this idea to more points. Lagrange formula. .Given (To,20),(1:11),...(In:Un). Construct polynomials Lo(x),L1(),...,Ln(x)of order n such that Ln.o(o)=1,Ln.o(1)=0,Ln.o(2)=0,...In.o(n)0. Ln,1(xo)=0,Ln,1(x1)=1,Lm,1(x2）=0,.Ln,1(xn)≠0. Ln.i(xi)=1,Ln.i(j)=0,Vj#i. Ln.n(0)=0,Ln,n(1)=0,...Ln.n(In-1)=0,Ln,n(n)#1. 5
Lecture note 3 Numerical Analysis • Computational cost is high: Gaussian Elimination O(n 3 ). • Find fast solvers of the linear equation or explicit expression of a. • Use “bases” other than {1, x, x2 , . . . , xn}. • Two points case: (x0, y0),(x1, y1). Construct two linear functions (polynomials of order 1) L1,0(x) and L1,1(x) such that: L1,0(x0) = 1 L1,0(x1) = 0 and L1,1(x0) = 0 L1,1(x1) = 1. Then, let P1(x) = y0L1,0(x) + y1L1,1(x) We have P1(x0) = y0L1,0(x0)+y1L1,1(x0) = y0 P1(x1) = y0L1,0(x1)+y1L1,1(x1) = y1. • Construct L1,0(x): L1,0(x) = x − x1 x0 − x1 . Use a graph to illustrate it!. Similarly, L1,1(x) = x − x0 x1 − x0 . • Generalize this idea to more points. Lagrange formula. • Given (x0, y0),(x1, y1), . . .(xn, yn). Construct polynomials L0(x), L1(x), . . . , Ln(x) of order n such that Ln,0(x0) = 1, Ln,0(x1) = 0, Ln,0(x2) = 0, . . . Ln,0(xn) 6= 0. Ln,1(x0) = 0, Ln,1(x1) = 1, Ln,1(x2) = 0, . . . Ln,1(xn) 6= 0. . . . Ln,i(xi) = 1, Ln,i(xj ) = 0, ∀j 6= i. . . . Ln,n(x0) = 0, Ln,n(x1) = 0, . . . Ln,n(xn−1) = 0, Ln,n(xn) 6= 1. 5
Lecture note 3 Numerical Analysis Then construct Pn(x)by B.()=∑Ln =0 We have Pn(i)=Vi. Find it！ .1,22,...,In are the n roots of Ln.o(),so Ln,0(x)=C(x-x1)(x-x2).(x-xn) where C is a constant. To find C,we use Ln.o(ro)=1. 1=Ln,0(xo）=C(x-x1)(x-x2).(x-xn) So 1 C=0-10-2)(0-n】） Therefore Ln,0(x)= (c-1)x-2).(红-xn）Π=1(e-x) (0-E1)(z0-x2).(0-xn)Π巧=1(0-) ·Similarly, Ln.i()= Π明=0j≠z-） Π=0，j≠ci-) Use a graph to illustrate it. .Lagrange formula:Given (ro:y0),(z1,),...,(En,Un).Define Π巧=0（c-） =0 Π写=0j≠(-） ·Example: zo 1,1 =2,r2 =4 and the correspondences yo =8,y 1 and y2=5. Find the polynomial p(x)of order 2 such that yo =p(xo),y=p(1),y2= p(r2) B(回=8c-2儿a-0+1g-1g-0+5红-a-2=3x2-16r+21 (1-2)1-4)(2-1)(2-4)(4-1)(4-2) Drawbacks of Lagrange interpolation. 6
Lecture note 3 Numerical Analysis • Then construct Pn(x) by Pn(x) = Xn i=0 yiLn,i(x). We have Pn(xi) = yi . Find it! • x1, x2, . . . , xn are the n roots of Ln,0(x), so Ln,0(x) = C(x − x1)(x − x2). . .(x − xn), where C is a constant. To find C, we use Ln,0(x0) = 1. 1 = Ln,0(x0) = C(x − x1)(x − x2). . .(x − xn). So C = 1 (x0 − x1)(x0 − x2). . .(x0 − xn) Therefore, Ln,0(x) = (x − x1)(x − x2). . .(x − xn) (x0 − x1)(x0 − x2). . .(x0 − xn) = Qn j=1(x − xj ) Qn j=1(x0 − xj ) • Similarly, Ln,i(x) = Qn j=0,j6=i (x − xj ) Qn j=0,j6=i (xi − xj ) • Use a graph to illustrate it. • Lagrange formula: Given (x0, y0), (x1, y1), . . . , (xn, yn). Define Pn(x) = Xn i=0 yi · Qn j=0,j6=i (x − xj ) Qn j=0,j6=i (xi − xj ) . • Example: x0 = 1, x1 = 2, x2 = 4 and the correspondences y0 = 8, y1 = 1 and y2 = 5. Find the polynomial p(x) of order 2 such that y0 = p(x0), y1 = p(x1), y2 = p(x2). P2(x) = 8(x − 2)(x − 4) (1 − 2)(1 − 4) +1 (x − 1)(x − 4) (2 − 1)(2 − 4) +5 (x − 1)(x − 2) (4 − 1)(4 − 2) = 3x 2−16x+21 • Drawbacks of Lagrange interpolation. 6
Lecture note 3 Numerical Analysis 1.Given z,need O(n2)to get Pn(x). 2.No part of the previous calculation can be used. Interpolation for additional values of x requires the same amount of effort as the first value. 3.If we want to add a new point (n+1,Un+1)or delete a point,all the “basis”Ln,i(z)should be changed. 1.3 Neville's method Idea:use previous calculation and generate recursively Lagrange polynomial ap- proximations. Definition 1 Let f be a function defined at ro,1,..,x2,..,In,and suppose that m1,m2,·,mk are k distinct integers,,with0≤mi≤n for each i.The Lagrange polynomial that agrees with f(r)at the k points imi,Im2,.,Imx is denoted Pm1,…,mk(c). Theorem 2 Let f be defined at xo,1,...,rk,and let rj and i be two distinct numbers in this set.Then P(z)= （工-工)P0.1.…j-1+1.…k(）-(r-x)D0.1…i-1i+1…k() 工i-Tj describes the k-th Lagrange polynomial that interpolates f at the k+1 points T0,T1,··,工k Algorithms ·INPUT:numbers z,xo,…,xn;values f(zo),…,f(rn)as thcolumn Q0.0,Q1.0,..,Qn,o of Q. .OUTPUT:the table Q with P()=Qn.n. 。Fori=1,2,…,n forj=1,2,…,i set Qig =Q1-Q11 ·Output(Q).Stop. 7
Lecture note 3 Numerical Analysis 1. Given x, need O(n 2 ) to get Pn(x). 2. No part of the previous calculation can be used. Interpolation for additional values of x requires the same amount of effort as the first value. 3. If we want to add a new point (xn+1, yn+1) or delete a point, all the “basis” Ln,i(x) should be changed. 1.3 Neville’s method Idea: use previous calculation and generate recursively Lagrange polynomial approximations. Definition 1 Let f be a function defined at x0, x1, · · · , x2, · · · , xn, and suppose that m1, m2, · · · , mk are k distinct integers, with 0 ≤ mi ≤ n for each i. The Lagrange polynomial that agrees with f(x) at the k points xm1 , xm2 , · · · , xmk is denoted Pm1,··· ,mk (x). Theorem 2 Let f be defined at x0, x1, · · · , xk, and let xj and xi be two distinct numbers in this set. Then P(x) = (x − xj )P0,1,··· ,j−1,j+1,··· ,k(x) − (x − xi)P0,1,··· ,i−1,i+1,··· ,k(x) xi − xj describes the k-th Lagrange polynomial that interpolates f at the k + 1 points x0, x1, · · · , xk. Algorithms • INPUT: numbers x, x0, · · · , xn; values f(x0), · · · , f(xn) as thcolumn Q0,0, Q1,0, · · · , Qn,0 of Q. • OUTPUT: the table Q with P(x) = Qn,n. • For i = 1, 2, · · · , n for j = 1, 2, · · · , i set Qi,j = (x−xi−j )Qi,j−1−(x−xi)Qi−1,j−1 xi−xj . • Output (Q). Stop. 7
Lecture note 3 Numerical Analysis 1.4 Newton's divided difference .method of undetermined coefficient:bases {1,x2,...,r"}. -Expand Pn(x)=ao +ax+a2x2+...+anz". -Need to solve a system linear equation to get ai. -Each basis doesn't depend on given points. Lagrange:bases {In.i()},Li()are polynomials of degree n and Ln,i(i)=1,Ln,i(j)=0,jti. -n回= -The corresponding coefficient is yi. -Each basis depends on o,1,...,i-1,i+1,...,In. -Less flexibility when a new point(+,)is added. Newton's divided difference -Each basis depends on part of the given points. -More flexibility when a new point (i+1,+1)is added. .Newton's divided difference:bases {x-zo,(x-zo)(x-z1),...,(x-zo)(x- 1)...(x-In-1)}.We assume Pn(x)=c0+c(x-xo)+c2(x-x0)(x-x1)+..+cn(x-x0)(x-x1)..(x-xn-1). We need to find ci. 。For co, 0=Pn(xo)=c0→c0=0 Then, Pn(x)=y0+c1(x-zo)+... ·ForC, 1=Pn(x1)=0+c1(x1-xo） So =欢-0 x1-T0 ·ForC2, 2=Pn(x2)=0+c(x2-x0）+c2(x2-xo)(c2-x1) So c2=2-0-a(2-x0) (x2-x0)(r2-工1) 8
Lecture note 3 Numerical Analysis 1.4 Newton’s divided difference • method of undetermined coefficient: bases {1, x, x2 , . . . , xn}. – Expand Pn(x) = a0 + a1x + a2x 2 + . . . + anx n. – Need to solve a system linear equation to get ai . – Each basis doesn’t depend on given points. • Lagrange: bases {Ln,i(x)}, Li(x) are polynomials of degree n and Ln,i(xi) = 1, Ln,i(xj ) = 0, j 6= i. – Ln,i(x) = Qn j=0,j6=i (x−xj ) Qn j=0,j6=i (xi−xj ) . – The corresponding coefficient is yi . – Each basis depends on x0, x1, . . . , xi−1, xi+1, . . . , xn. – Less flexibility when a new point (xi+1, yi+1) is added. • Newton’s divided difference – Each basis depends on part of the given points. – More flexibility when a new point (xi+1, yi+1) is added. • Newton’s divided difference: bases {x − x0,(x − x0)(x − x1), . . . ,(x − x0)(x − x1). . .(x − xn−1)}. We assume Pn(x) = c0+c1(x−x0)+c2(x−x0)(x−x1)+. . .+cn(x−x0)(x−x1). . .(x−xn−1). We need to find ci . • For c0, y0 = Pn(x0) = c0 =⇒ c0 = y0. Then, Pn(x) = y0 + c1(x − x0) + . . . • For c1, y1 = Pn(x1) = y0 + c1(x1 − x0) So c1 = y1 − y0 x1 − x0 • For c2, y2 = Pn(x2) = y0 + c1(x2 − x0) + c2(x2 − x0)(x2 − x1) So c2 = y2 − y0 − c1(x2 − x0) (x2 − x0)(x2 − x1) 8
Lecture note 3 Numerical Analysis ·And so on. ·Example:(xo,0）=(1,8),(c1,h)=(2,1),(r2,2)=(4,5). C0=0=8. c1=(y1-bo)/(x1-x0）=-7 2=2-h--w=3 (x2-x0)(x2-x1) Pn(x)=8-7(x-1)+3(x-1)(x-2) .There exists a closed form for ci,i=0,...,n. 。Define fz=欢 f儿，x+]= f[z+]-f儿] 工i+1-Ti fz,xi+1,+2= f[i+1,i+2]-f[ri,i+1] 工i+2-工i f儿，+1，x+2,+3l= f[z+1,x+2,x+3-f[x,xi+1,x+2] 工i+3-ri The prove will be delayed. co=f[zo] c1=f[ro,1] ci=f[ro,z1,...,xi] cn f[o,1,...,In] divided difference table: 1f] fr1,2] T2 f(2] fx1,x2,x3] f[z2,x3] f[1,x2,T3,T4] fx2,x3,x4] f[x1,x2,x3,x4,x5 f[3,4] fx2,x3,E4,I6] 4 f4] fz3,x4,T6] f[4,x5] 工5f[zsl 9
Lecture note 3 Numerical Analysis • And so on. • Example: (x0, y0) = (1, 8), (x1, y1) = (2, 1), (x2, y2) = (4, 5). c0 = y0 = 8. c1 = (y1 − b0)/(x1 − x0) = −7 c2 = y2 − y0 − b1(x2 − x0) (x2 − x0)(x2 − x1) = 3. Pn(x) = 8 − 7(x − 1) + 3(x − 1)(x − 2) • There exists a closed form for ci , i = 0, . . . , n. • Define f[xi ] = yi f[xi , xi+1] = f[xi+1] − f[xi ] xi+1 − xi f[xi , xi+1, xi+2] = f[xi+1, xi+2] − f[xi , xi+1] xi+2 − xi f[xi , xi+1, xi+2, xi+3] = f[xi+1, xi+2, xi+3] − f[xi , xi+1, xi+2] xi+3 − xi • The prove will be delayed. • c0 = f[x0] c1 = f[x0, x1] . . . ci = f[x0, x1, . . . , xi ] . . . cn = f[x0, x1, . . . , xn] • divided difference table: x1 f[x1] f[x1, x2] x2 f[x2] f[x1, x2, x3] f[x2, x3] f[x1, x2, x3, x4] x3 f[x3] f[x2, x3, x4] f[x1, x2, x3, x4, x5] f[x3, x4] f[x2, x3, x4, x5] x4 f[x4] f[x3, x4, x5] f[x4, x5] x5 f[x5] 9
Lecture note 3 Numerical Analysis ·Example: 18 -7 21 3 2 45 ·Example:: Give the example of interpolating f(x)=Vr on 2,2.1,2.2,2.3,2.4. 1.4.1 Properties and Derivation of Newton's divided differ- ence We recall Newton's divided difference formula: Pn(z)=f[o]+f[xo;1](x-zo)+f[xo:1,2](x-zo)(x-1)+ ...f[zo,...,In](x zo)(-z1)...(x-In-1) Note that adding a new point (n+1,Un+1)we have to add a new line in the divided difference table to get f[ro,...,n+]. add a new term f[ro,...,n+](-zo)(x-z1)...(x-En)into the expression of Pn(r)to get Pn+1. The derivation of Newton divided difference formula is given int he the following theorem Theorem 3 (Theorem 3.1)The unique polynomial of order n that passes through (xo,y0),(1,11),....(In,yn)is given by Pn()=f[zo]+f[ro,z](-zo)+f[o:21,z2](r-z0)(x-z1) +.+fxo,,xn](x-o)(x-x1).(x-xn-1) Proof.We prove it by induction.The unique polynomial of order 0 that passes through (o,y0)is obviously Po(x)=0=f儿rol. Suppose that the polynomial P()of order <k that passes through(zo,yo),(1,),...,(k,y) is Pk()=f[zo]+f[zo;x1](x-zo)+f[zo,z1,z2](x-zo)(z-z1)+ .+f[z0,,Ek(x-x0)(x-x1).(z-xk-1) Write P+(),the unique polynomial of order <k that passes through (ro,y0),(1,41),...,(k,),(k+ by Pk+1(x)=f[xo]+f[ro,1](x-x0)+f[ro,1,x2](x-xo)(-1) +.+fx0,;xk(x-x0)(x-c1).(x-xk-1） +C(x-x0)(x-x1).(x-Ek-1)(z-Ek) 10
Lecture note 3 Numerical Analysis • Example: 1 8 -7 2 1 3 2 4 5 • Example: Give the example of interpolating f(x) = √ x on 2, 2.1, 2.2, 2.3, 2.4. 1.4.1 Properties and Derivation of Newton’s divided difference We recall Newton’s divided difference formula: Pn(x) = f[x0] + f[x0, x1](x − x0) + f[x0, x1, x2](x − x0)(x − x1)+ . . . + f[x0, . . . , xn](x − x0)(x − x1). . .(x − xn−1) Note that adding a new point (xn+1, yn+1) we have to • add a new line in the divided difference table to get f[x0, . . . , xn+1]. • add a new term f[x0, . . . , xn+1](x−x0)(x−x1). . .(x−xn) into the expression of Pn(x) to get Pn+1. The derivation of Newton divided difference formula is given int he the following theorem Theorem 3 (Theorem 3.1) The unique polynomial of order ≤ n that passes through (x0, y0),(x1, y1), . . . ,(xn, yn) is given by Pn(x) = f[x0] + f[x0, x1](x − x0) + f[x0, x1, x2](x − x0)(x − x1) + . . . + f[x0, . . . , xn](x − x0)(x − x1). . .(x − xn−1) Proof. We prove it by induction. The unique polynomial of order 0 that passes through (x0, y0) is obviously P0(x) = y0 = f[x0]. Suppose that the polynomial Pk(x) of order ≤ k that passes through (x0, y0),(x1, y1),. . . ,(xk, yk) is Pk(x) = f[x0] + f[x0, x1](x − x0) + f[x0, x1, x2](x − x0)(x − x1)+ . . . + f[x0, . . . , xk](x − x0)(x − x1). . .(x − xk−1) Write Pk+1(x), the unique polynomial of order ≤ k that passes through (x0, y0),(x1, y1), . . . ,(xk, yk),(xk+1, yk+1) by Pk+1(x) = f[x0] + f[x0, x1](x − x0) + f[x0, x1, x2](x − x0)(x − x1) + . . . + f[x0, . . . , xk](x − x0)(x − x1). . .(x − xk−1) + C(x − x0)(x − x1). . .(x − xk−1)(x − xk) 10