Topic : Spline Interpolation
Simulation : Graphical Simulation of the Method
Language : Mathmatica 4.1
Authors : Nathan Collier, Autar Kaw
Date : 18 July 2002
Abstract : This simulation illustrates the spline interpolation method. Given 'n' data points of y versus x, and then required to find the value of 'y' at a particular 'x', you are asked to used quadrature spline interpolation using Mathematica's internal function are also given.

ClearAll;

INPUTS: Enter the Following

Array of x data

x = {10, 0.000001, 20, 15, 30, 22.5};

Note: If you have a zero in the x matrix, the simulation will not work. Try making the zero a near zero.

Array of y data

y = {227.04, 0, 517.35, 362.78, 901.67, 602.97};

Value of x at which y is desired

xdesired := 16

SOLUTION

The following sorts the x and y arrays in ascending order.

nn := Abs[Dimensions[x]]

n := nn[[1]]

xy = Table[0, {i, n}, {j, 2}];

Do[xy[[i, 1]] = x[[i]]; 
      xy[[i, 2]] = y[[i]], {i, 1, n}]

xy = Sort[xy];

x = xy[[All, 1]];
y = xy[[All, 2]];

The sorted x and y values

x // MatrixForm

y // MatrixForm

data = 
        ListPlot[xy, PlotStyle -> PointSize[0.02], 
          PlotLabel -> "Given y vs x data points", 
          TextStyle -> {FontSize -> 11}, 
          PlotRange -> All];

Linear spline interpolation

m[i_, z_] := 
      y[[i + 1]] + (y[[i + 1]] - y[[i]])/(x[[i + 1]] - 
                x[[i]])*(z - x[[i + 1]])

flinear[z_] := 
      If[z < x[[2]], i = 1; m[i, z], 
        Do[If[z ≤ x[[j + 1]] && z > x[[j]], i = j], {j, 1, 
            n - 1}]; m[i, z]]

Value of function at desired value

flinear[xdesired]

fprev = %;

fline = Plot[flinear[z], {z, Min[x], Max[x]}];

desire = 
        ListPlot[{{xdesired, flinear[xdesired]}}, 
          PlotStyle -> PointSize[0.03]];

Show[data, fline, desire];

Quadratic spline interpolation

The following assembles the matrix whose inverse is needed to solve for the coefficients of the polynomial splines that fit the data.

A = Table[0, {i, 3*n - 3}, {j, 3*n - 3}];

Do[Do[
        Do[A[[2*i - 1 + j + 1, 
                3*i - 3 + 1 + 
                  k]] = (x[[i - 1 + j + 1]])^(k), {k, 0, 
            2}], {j, 0, 1}], {i, 1, n - 1}]

Do[Do[
        Do[A[[2*(n - 1) + i + 1, 
                3*i - 2 + k + j*3 + 1]] = (-1)^
                j*(2*x[[i + 1]])^k, {k, 0, 1}], {j, 0, 1}], {i, 1, 
        n - 2}]

<< LinearAlgebra`MatrixManipulation`

A[[1, 3]] = 1;

A // MatrixForm

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

Y = Table[0, {i, 3*n - 3}];

Do[Do[Y[[2*(i + 1) + j]] = y[[i + j + 1]], {j, 0, 1}], {i, 
        0, n - 2}]

Y // MatrixForm

Co = LinearSolve[A, Y];

Co // MatrixForm

fquadratic[z_] := 
      If[z ≤ x[[2]], 
        Co[[1]] + Co[[2]]*z + Co[[3]]*z^2, 
        Do[If[z ≤ x[[i + 2]] && z > x[[i + 1]], d = 0; 
            Do[d = d + Co[[3*i + j + 1]]*z^j, {j, 0, 2}]], {i, 1, 
            n - 2}]; d]

Value of function at desired value

fquadratic[xdesired]

fprev = %;

fline = Plot[fquadratic[z], {z, Min[x], Max[x]}];

desire = 
        ListPlot[{{xdesired, fquadratic[xdesired]}}, 
          PlotStyle -> PointSize[0.03]];

Show[data, fline, desire];