Application Center - Maplesoft
  Submit your work
Application Center Applications Spline Method of Interpolation Preview

App Preview:

Spline Method of Interpolation

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

Learn about Maple
Download Application

 

splineInterp.mws

Spline 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 spline method of interpolation.  We limit this worksheet to linear and quadratic spline interpolation.

Introduction

The Spline me thod of interpolation (for detailed explanation, you can read the textbook notes and examples, and see a Power Point Presentation) is illustrated .   Given 'n' data points of 'y' versus 'x', it is required to find the value of y at a particular value of x using linear and quadratic splines.

>    restart;

>    with(LinearAlgebra):

Section I : 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  a 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.

This function will sort the data matrix into ascending order and puts them into a new matrix.

>    n:=nops(xy):

>    for passnum from 1 to n-1 do
  for i from 1 to n-passnum do
    if xy[i,1]>xy[i+1,1] then
      temp1:=xy[i,1];
      temp2:=xy[i,2];
      xy[i,1]:=xy[i+1,1];
      xy[i,2]:=xy[i+1,2];
      xy[i+1,1]:=temp1;
      xy[i+1,2]:=temp2;
    end if:
  end do:
end do:

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

>    x,y;

Matrix(%id = 18936600), Matrix(%id = 18938880)

>    ranger:=proc(ar,n)
local i,xmin,xmax,xrange;
xmin:=ar[1,1]:
xmax:=ar[1,1]:
for i from 1 to n do
if ar[i,1] > xmax then xmax:=ar[i,1] end if;
if ar[i,1] < xmin then xmin:=ar[i,1] 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 spline interpolation

Given [Maple OLE 2.0 Object]  , fit linear splines to the data.  This simply involves forming the consecutive data through straight lines.  So the data is sorted into an ascending order; the linear splines are given by   [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]

                                                    ..........

                                                    [Maple OLE 2.0 Object]             [Maple OLE 2.0 Object]

>    flinear:=proc(z)
  local d,i;
  if z < x[2,1] then
    d:=y[2,1]+(y[2,1]-y[1,1])/(x[2,1]-x[1,1])*(z-x[2,1]);
  else
    for i from 2 to n do
      if (z > x[i,1]) and (z <= x[i+1,1]) then
        d:=y[i,1]+(y[i+1,1]-y[i,1])/(x[i+1,1]-x[i,1])*(z-x[i,1]);
      end if;
    end do;
  end if;
  d;
end proc:

Value of function at desired value of X is

>    flinear(xdesired);

393.6940000

>    fprev:=%:

Plotting the linear spline interpolant and the value of Y for the desired X

>    plot([xy,[[xdesired,flinear(xdesired)]],flinear],ranger(x,n),style=[POINT,POINT,LINE],color=[BLUE,RED,BLACK],symbol=[CIRCLE,CROSS],symbolsize=[30,40],thickness=2,title="Linear spline interpolation");

[Maple Plot]

Section IV: Quadratic interpolation

In these splines, a quadratic polynomial approximates the data between two consecutive data points.  Given [Maple OLE 2.0 Object]  , fit quadratic splines through the data.  The splines are given by
                      [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] [Maple OLE 2.0 Object]              
There are 3n such coefficients   [Maple OLE 2.0 Object] [Maple OLE 2.0 Object]  1, 2, , n; [Maple OLE 2.0 Object] [Maple OLE 2.0 Object]  1, 2, , n; [Maple OLE 2.0 Object] [Maple OLE 2.0 Object]  1, 2, , n
To find '3n' unknowns, one needs to set up '3n' equations and then simultaneously solve them.

To know more about how to setup these '3n' equations, please click on textbook notes. After setting up these equations, they are solved using the matrix method. The following function assembles the matrix whose inverse is needed to solve for the coefficients of the polynomial splines that fits the data.

>    A:=Matrix(3*(n-1),3*(n-1));

A := Matrix(%id = 19560236)

>    for i from 0 to 3*(n-1)-1 do
  for j from 0 to 3*(n-1)-1 do
    A[i+1,j+1]:=0;
  end do;
end do;

for i from 1 to n-1 do
  for j from 0 to 1 do
    for k from 0 to 2 do
      A[2*i-1+j+1,3*i-3+k+1]:=x[i+j,1]^k;
    end do;
  end do;
end do;

for i from 1 to n-2 do
  for j from 0 to 1 do
    for k from 0 to 1 do
      A[2*(n-1)+i+1,3*i-2+k+j*3+1]:=(-1)^j*(2*x[i+1,1])^k;
    end do;
  end do;
end do;

A[1,3]:=1:

The following gives the matrix which is required to solve for the coefficients of the quadratic splines

A;

Matrix(%id = 19560236)

This assembles the Y matrix also needed to determine the coefficients of the quadratic splines.

>    Y:=Matrix(3*(n-1),1):
for i from 0 to n-2 do
  for j from 0 to 1 do
    Y[2*(i+1)+j,1]:=y[i+j+1,1];
  end do;
end do;
Y;

Matrix(%id = 19612512)

Solving for the coefficients, we get

>    C:=MatrixInverse(A) . Y;

C := Matrix(%id = 19810968)

>    fquad:=proc(z)
local d,i,j:
if z <= x[2,1] then
  d:=C[1,1]+C[2,1]*z+C[3,1]*z^2:
else
  for i from 2 to n-1 do
    if z <= x[i+1,1] and z > x[i,1] then
      d:=0:
      for j from 0 to 2 do
        d:=d+C[3*(i-1)+j+1,1]*z^j:
      end do:
    end if:
  end do:
end if:
d;
end proc:

Value of function at desired value of X is

>    fquad(xdesired);

394.2364000

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

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

[Maple Plot]

Section V: Cubic spline interpolation

The algorithm of cubic spline interpolation is not shown.  However, we are using the Maple function to conduct the cubic spline interpolation,

>    fcubic:=t->spline(convert(x,vector),convert(y,vector),t,cubic):

Computing the value of the function at the desired value of X,

>    fcubic(xdesired);

392.1542016

Plotting the cubic spline interpolant and the value of Y at the desired calue of X,

>    plot([xy,[[xdesired,fcubic(xdesired)]],fcubic],ranger(x,n),style=[POINT,POINT,LINE],color=[BLUE,RED,BLACK],symbol=[CIRCLE,CROSS],symbolsize=[30,40],thickness=2,title="Cubic spline 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 linear and quadratic spline 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_spline.mws

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