• Introduction
• Section I : Input Data.
• Section II : Big scary functions.
• Section III: Linear Interpolation.
• Section IV: Quadratic Interpolation (Second order polynomial)
• Section V: Cubic Interpolation (Third order polynomial)
• Section VI: Conclusion.

Newton Divided Difference Method of Interpolation

2003 Nathan Collier, Autar Kaw, Jai Paul , Michael Keteltas, University of South Florida , kaw@eng.usf.edu , http://numericalmethods.eng.usf.edu/mws

NOTE: This worksheet demonstrates the use of Maple to illustrate the Newton's Divided Difference Method of interpolation.  We limit this worksheet to using first, second, and third order polynomials.

Introduction

The Newton's Divided Difference Polynomial method of interpolation (is based on the following.  (For a detailed explanation, you can read the textbook notes and examples, or see a Power Point Presentation) 

The general form of the Newton's divided difference polynomial for data points  is given as

where n is the order of the polynomial that approximates the function and,

               

               

               

               ...........

               

               

where the definition of divided difference is

The examples that follow illustrate the method in a descriptive manner. For further explanations, click on text notes for the method on the website.

>

restart;

Section I : Input Data.

The following is the array of x-y data which is used to interpolate. It is obtained from the physical problem  of velocity of  rocket (y-values) vs. time (x-values) data.   We are asked to find the velocity at an intermediate point of x=16.

>	xy:=[[10,227.04],[0,0],[20,517.35],[15,362.78],[30,901.67],[22.5,602.97]]:

Value of X at which Y is desired

>	xdesired:=16:

Section II : Big scary functions.

>	n:=nops(xy):

>	for i from 1 to n do   x[i]:=xy[i,1];   y[i]:=xy[i,2]; end do:

The following function considers the x data and selects those data points which are close to the desired x value based on the least absolute difference between the x values and the desired x value. This function  selects the two  closest data points that bracket the desired value of x. It first picks the closest data point to the desired x value. It then checks if this value is less than or greater than the desired value.If it is less, then, it selects the data point which is greater than the desired value and also the closest, and viceversa.

Finds the absolute difference between the X values and the desired X value.

>	for i from 1 to n do   co[i]:=abs(x[i]-xdesired); end do:

Identifies the X value with the least absolute difference.
c:=co[1]:
for i from 2 to n do
  if c > co[i] then
    c:=co[i];
    ci:=i;
  end if:
end do:

If the value with the least absolute difference is less than the desired value, then it selects the closest data point greater than the desired value to bracket the desired value.
if x[ci] < xdesired then
  q:=1;
  for i from 1 to n do
    if x[i] > xdesired then
      nex[q]:=x[i];
      q:=q+1;
    end if;
  end do;

  b:=nex[1]:
  for i from 2 to q-1 do
    if b > nex[i] then
      b:=nex[i];
    end if:
  end do:

  for i from 1 to n do
    if x[i]=b then bi:=i end if;
  end do;
end if:

If the value with the least absolute difference is greater than the desired value, then it selects the closest data point less than the desired value to bracket the desired value.
if x[ci] > xdesired then
  q:=1;
  for i from 1 to n do
    if x[i] < xdesired then
      nex[q]:=x[i];
      q:=q+1;
    end if;
  end do;

  b:=nex[1]:
  for i from 2 to q-1 do
    if b < nex[i] then
      b:=nex[i];
    end if:
  end do:

  for i from 1 to n do
    if x[i]=b then bi:=i end if;
  end do;
end if:

firsttwo:=<ci,bi>:

If more than two values are desired, the same procedure as above is followed to choose the 2 datapoints which bracket the desired value. In addition, the following function selects the subsequent values that are closest to the desired value and puts all the values into a matrix, maintaining the original data order.

>

for i from 1 to n do   A[i,2]:=i;   A[i,1]:=co[i]; end do: for passnum from 1 to n-1 do   for i from 1 to n-passnum do     if A[i,1]>A[i+1,1] then       temp1:=A[i,1];       temp2:=A[i,2];       A[i,1]:=A[i+1,1];       A[i,2]:=A[i+1,2];       A[i+1,1]:=temp1;       A[i+1,2]:=temp2;     end if:   end do: end do: for i from 1 to n do   A[i,3]:=i; end do: for passnum from 1 to n-1 do   for i from 1 to n-passnum do     if A[i,2]>A[i+1,2] then       temp1:=A[i,1];       temp2:=A[i,2];       temp3:=A[i,3];       A[i,1]:=A[i+1,1];       A[i,2]:=A[i+1,2];       A[i,3]:=A[i+1,3];       A[i+1,1]:=temp1;       A[i+1,2]:=temp2;       A[i+1,3]:=temp3;     end if:   end do: end do: for i from 1 to n do   d[i]:=A[i,3]; end do: if d[firsttwo[2]]<>2 then   temp:=d[firsttwo[2]];   d[firsttwo[2]]:=1;   for i from 1 to n do     if i <> firsttwo[1] and i <> firsttwo[2] and d[i] <= temp then        d[i]:=d[i]+1;     end if:   end do;   d[firsttwo[1]]:=1; end if:

>	ranger:=proc(ar,n)   local i,xmin,xmax,xrange;   xmin:=ar[1]:   xmax:=ar[1]:   for i from 1 to n do     if ar[i] > xmax then xmax:=ar[i] end if;     if ar[i] < xmin then xmin:=ar[i] end if;   end do;   xrange:=xmin..xmax; end proc:

Plotting the given values of X and Y.

>	plot(xy,ranger(x,n),style=POINT,color=BLUE,symbol=CIRCLE,symbolsize=30);

Section III: Linear Interpolation.

Given      fit a linear interpolant through the data.  Noting    and   , assume the linear interpolant,   is given by , where  and .

Hence, the two closest data points to the desired value are chosen in this method.

>	datapoints:=2:

>	p:=1: for i from 1 to n do   if d[i] <= datapoints then     xdata[p]:=x[i];     ydata[p]:=y[i];     p:=p+1;   end if end do:

>	entries(xdata);

>	entries(ydata);

Calculating coefficients of Newton's Divided difference polynomial

>	b0:=ydata[1];

>	b1:=(ydata[2]-ydata[1])/(xdata[2]-xdata[1]);

Newton's divided difference formula for linear interpolation is

>	flinear:=z->b0+b1*(z-xdata[1]);

Substituting the value of the desired X value for z in the above equation, the corresponding Y value is found.

>	eval(flinear(z),z=xdesired);

>

fprev:=%:

Plotting the Linear interpolant and the value of Y for the desired X

>	plot([[t,eval(flinear(z),z=t),t=ranger(xdata,datapoints)],xy,[[xdesired,eval(flinear(z),z=xdesired)]]],z=ranger(x,n),style=[LINE,POINT,POINT],color=[RED,BLUE,BLUE],symbol=[CROSS,CIRCLE],symbolsize=[40,30],thickness=3,title="Linear interpolation");

Section IV: Quadratic Interpolation (Second order polynomial)

Given   and  fit a quadratic interpolant through the data.  Noting  and  assume the quadratic interpolant given by , where

                   

                   

                    .

Hence, the three closest data points to the desired value are chosen in this method.

>	datapoints:=3:

>	p:=1: for i from 1 to n do   if d[i] <= datapoints then     xdata[p]:=x[i];     ydata[p]:=y[i];     p:=p+1;   end if end do:

>	entries(xdata);

>	entries(ydata);

Calculating coefficients of Newton's Divided difference polynomial

>	b0:=ydata[1];

>	b1:=(ydata[2]-ydata[1])/(xdata[2]-xdata[1]);

>	b2:=((ydata[3]-ydata[2])/(xdata[3]-xdata[2])-(ydata[2]-ydata[1])/(xdata[2]-xdata[1]))/(xdata[3]-xdata[1]);

Newton's divided difference formula for quadratic interpolation is

>	fquad:=z->b0+b1*(z-xdata[1])+b2*(z-xdata[1])*(z-xdata[2]);

Substituting the value of the desired X value for z in the above equation, the corresponding Y value is found.

>	fquad(xdesired);

>

fnew:=%:

>

The absolute percentage relative approximate error between the first order and second order interpolation is

>	epsilon:=abs((fnew-fprev)/fnew)*100;

The number of significant digits at least correct in the solution is

>	sigdig:=floor(2-log10(epsilon/0.5));

>	fprev:=fnew:

Plotting the quadratic interpolant and the value of Y for the desired X

>	plot([[t,fquad(t),t=ranger(xdata,datapoints)],xy,[[xdesired,fquad(xdesired)]]],z=ranger(x,n),style=[LINE,POINT,POINT],color=[RED,BLUE,BLUE],symbol=[CROSS,CIRCLE],symbolsize=[40,30],thickness=3,title="Quadratic interpolation");

Section V: Cubic Interpolation (Third order polynomial)

The same method as above is used for 4 data points and the corresponding coefficients,  are calcualted and evaluated in the formula as shown below to find the corresponding Y value for the desired X value.

Hence, the four closest data points to the desired value are chosen in this method.

>	datapoints:=4:

>	p:=1: for i from 1 to n do   if d[i] <= datapoints then     xdata[p]:=x[i];     ydata[p]:=y[i];   p:=p+1;   end if end do:

>	entries(xdata);

>	entries(ydata);

Calculating coefficients of Newton's Divided difference polynomial

>	b0:=ydata[1];

>	b1:=(ydata[2]-ydata[1])/(xdata[2]-xdata[1]);

>	b2:=((ydata[3]-ydata[2])/(xdata[3]-xdata[2])-(ydata[2]-ydata[1])/(xdata[2]-xdata[1]))/(xdata[3]-xdata[1]);

>	b3:=(((ydata[4]-ydata[3])/(xdata[4]-xdata[3])-(ydata[3]-ydata[2])/(xdata[3]-xdata[2]))/(xdata[4]-xdata[2])-b2)/(xdata[4]-xdata[1]);

Newton's divided difference formula for cubic interpolation is

>	fcubic:=z->b0+b1*(z-xdata[1])+b2*(z-xdata[1])*(z-xdata[2])+b3*(z-xdata[1])*(z-xdata[2])*(z-xdata[3]): fcubic(z);

Substituting the value of the desired X value for z in the above equation, the corresponding Y value is found.

>	fcubic(xdesired);

>

fnew:=%:

The absolute percentage relative approximate error between the second order and third order interpolation is

>	epsilon:=abs((fnew-fprev)/fnew)*100;

The number of significant digits at least correct in the solution

>	sigdig:=floor(2-log10(epsilon/0.5));

Plotting the cubic interpolant and the value of Y for the desired X

>	plot([[t,fcubic(t),t=ranger(xdata,datapoints)],xy,[[xdesired,fcubic(t)]]],z=ranger(x,n),style=[LINE,POINT,POINT],color=[RED,BLUE,BLUE],symbol=[CROSS,CIRCLE],symbolsize=[40,30],thickness=3,title="Cubic interpolation");

>

>

Section VI: Conclusion.

Maple helped us to apply our knowledge of numerical methods of interpolation to find the value of y at a particular value of x using first, second, and third order Newton's Divided Difference Polynomial method of interpolation. Using Maple functions and plotting routines made it easy to illustrate this method.

References

[1] Nathan Collier, Autar Kaw, Jai Paul , Michael Keteltas, Holistic Numerical Methods Institute, See http://numericalmethods.eng.usf.edu/mws/gen/05inp/mws_gen_inp_sim_ndd.mws

http://numericalmethods.eng.usf.edu/mws/gen/05inp/mws_gen_inp_txt_ndd.pdf

Disclaimer :  While every effort has been made to validate the solutions in this worksheet, University of South Florida and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.