Application Center - Maplesoft

App Preview:

Newton Divided Difference Method of Interpolation

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

newtonDivDiff.mws

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 [Maple OLE 2.0 Object] data points [Maple OLE 2.0 Object]  is given as

[Maple OLE 2.0 Object] where n is the order of the polynomial that approximates the function and,

                [Maple OLE 2.0 Object]

                [Maple OLE 2.0 Object]

                [Maple OLE 2.0 Object]

               ...........

                [Maple OLE 2.0 Object]

                [Maple OLE 2.0 Object]

where the definition of [Maple OLE 2.0 Object] divided difference is [Maple OLE 2.0 Object] [Maple OLE 2.0 Object]

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);

[Maple Plot]

Section III: Linear Interpolation.

Given   [Maple OLE 2.0 Object] [Maple OLE 2.0 Object]    fit a linear interpolant through the data.  Noting   [Maple OLE 2.0 Object]  and   [Maple OLE 2.0 Object] , assume the linear interpolant, [Maple OLE 2.0 Object]   is given by [Maple OLE 2.0 Object] , where [Maple OLE 2.0 Object]  and [Maple OLE 2.0 Object] .

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);

[20], [15]

>    entries(ydata);

[517.35], [362.78]

Calculating coefficients of Newton's Divided difference polynomial

>    b0:=ydata[1];

b0 := 517.35

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

b1 := 30.91400000

Newton's divided difference formula for linear interpolation is

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

flinear := proc (z) options operator, arrow; b0+b1*(z-xdata[1]) end proc

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);

393.6940000

>    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");

[Maple Plot]

Section IV: Quadratic Interpolation (Second order polynomial)

Given   [Maple OLE 2.0 Object] [Maple OLE 2.0 Object] and [Maple OLE 2.0 Object]  fit a quadratic interpolant through the data.  Noting [Maple OLE 2.0 Object] [Maple OLE 2.0 Object] [Maple OLE 2.0 Object]  and [Maple OLE 2.0 Object]  assume the quadratic interpolant [Maple OLE 2.0 Object] given by [Maple OLE 2.0 Object] , where

                    [Maple OLE 2.0 Object]

                    [Maple OLE 2.0 Object]

                    [Maple OLE 2.0 Object] .

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);

[10], [20], [15]

>    entries(ydata);

[227.04], [517.35], [362.78]

Calculating coefficients of Newton's Divided difference polynomial

>    b0:=ydata[1];

b0 := 227.04

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

b1 := 29.03100000

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

b2 := .3766000000

Newton's divided difference formula for quadratic interpolation is

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

fquad := proc (z) options operator, arrow; b0+b1*(z-xdata[1])+b2*(z-xdata[1])*(z-xdata[2]) end proc

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

>    fquad(xdesired);

392.1876000

>    fnew:=%:

>   

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

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

epsilon := .3841018941

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

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

sigdig := 2

>    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");

[Maple Plot]

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);

[10], [20], [15], [22.5]

>    entries(ydata);

[227.04], [517.35], [362.78], [602.97]

Calculating coefficients of Newton's Divided difference polynomial

>    b0:=ydata[1];

b0 := 227.04

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

b1 := 29.03100000

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

b2 := .3766000000

>    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]);

b3 := .5434666560e-2

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);

-63.2700000+29.03100000*z+.3766000000*(z-10)*(z-20)+.5434666560e-2*(z-10)*(z-20)*(z-15)

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

>    fcubic(xdesired);

392.0571680

>    fnew:=%:

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

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

epsilon := .3326861760e-1

The number of significant digits at least correct in the solution

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

sigdig := 3

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");

[Maple Plot]

>   

>   

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.