• 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.

Direct Method of Polynomial 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 direct method of interpolation.  We limit this worksheet to using first, second, and third order polynomials.

Introduction

The direct method of interpolation (for detailed explanation, you can read the textbook notes and examples, and see a Power Point Presentation) is based on the following.

Given 'n+1' data points of y vs. x form, fit a polynomial of order 'n' as given below

                            (1)

through the data, where  are real constants.    Since 'n+1' values of y are given at 'n+1' values of x, one can write 'n+1' equations.  Then the 'n+1' constants can be found by solving the 'n+1' simultaneous linear equations.  To find the interpolated value of 'y' at a given value of 'x', simply substitute the value of 'x' in equation (1).
Since the number of data points may be more than you need for a particular order of interpolation, one has to first pick the needed data points, and then one can use the chosen points to interpolate the data.

>	restart; with(LinearAlgebra):

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.  The closeness is based on the least absolute difference between the x data 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 data points 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.

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

We then set up equations to find coefficients of the linear interpolant

>	M := Matrix([[1,xdata[1]],[1,xdata[2]]]);

>	m := MatrixInverse(M);

The Coefficients of the linear interpolant are,

>	a := m . Matrix([[ydata[1]],[ydata[2]]]);

Hence, the equation of the linear interpolant is,

>	flinear:=z->a[1,1]+a[2,1]*z;

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)

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

We then set up equations to find coefficients of the quadratic interpolant

>	M:=Matrix([[1,xdata[1],(xdata[1])^2],[1,xdata[2],(xdata[2])^2],[1,xdata[3],(xdata[3])^2]]);

>	m:=MatrixInverse(M);

The coefficients of quadratic interpolant are,

>	a:= m . Matrix([[ydata[1]],[ydata[2]],[ydata[3]]]);

Hence, the equation of the quadratic interpolant is,

>	fquad:=z->a[1,1]+a[2,1]*z+a[3,1]*z^2;

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

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

>

fnew:=%:

The absolute percentage relative approximate error between the first order and the second order interpolation values 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,eval(fquad(z),z=t),t=ranger(xdata,datapoints)],xy,[[xdesired,eval(fquad(z),z=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 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);

We then set up equations to find coefficients of the cubic interpolant

>	M:=Matrix([[1,xdata[1],(xdata[1])^2,(xdata[1])^3],[1,xdata[2],(xdata[2])^2,(xdata[2])^3],[1,xdata[3],(xdata[3])^2,(xdata[3])^3],[1,xdata[4],(xdata[4])^2,(xdata[4])^3]]);

>	m:=MatrixInverse(M);

The coefficients of cubic interpolant are,

>	a:= m. Matrix([[ydata[1]],[ydata[2]],[ydata[3]],[ydata[4]]]);

Hence, the equation of the cubic interpolant is,

>	fcubic:=z->a[1,1]+a[2,1]*z+a[3,1]*z^2+a[4,1]*z^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 the third 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));

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

>	plot([[t,fcubic(t),t=ranger(xdata,datapoints)],xy,[[xdesired,fcubic(xdesired)]]],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 direct 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_direct.mws

http://numericalmethods.eng.usf.edu/mws/gen/05inp/mws_gen_inp_txt_direct.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.