Topic : Direct Method
Simulation : Graphical Simulation of the Method
Language : Mathmatica 4.1
Authors : Nathan Collier, Autar Kaw
Date : 15 July 2002
Abstract : This simulation illustrates the direct method of interpolation. Given n data points of y versus x, you are then required to find the value of y at a particular value of x using first, second, and third order interpolation. So on has to first pick the needed data point, and then use those points to interpolate the data

Clear[x, y, xData, yData]

INPUTS: Enter the Following

Array of x data

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

Array of y data

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

xdesired := 16

n := Abs[Dimensions[x]][[1]]

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

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

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

SOLUTION

The following considers the x and y data and selects the two closest data points that bracket the desired value of x.

comp := Abs[x - xdesired]

c := Min[comp]

Do[If[comp[[i]] == c, ci = i], {i, 1, n}]

R = Table[0, {i, 1, n}];

If[x[[ci]] < xdesired, q = 1; 
      Do[If[x[[i]] > xdesired, R[[q]] = x[[i]]; q = q + 1], {i, 
          1, n}]; RR = Table[RR[[i]] = R[[i]], {i, 1, q - 1}]; 
      b = Min[RR]; 
      Do[If[x[[i]] == b, bi = i], {i, 1, n}]]

If[x[[ci]] > xdesired, q = 1; 
      Do[If[x[[i]] < xdesired, R[[q]] = x[[i]]; q = q + 1], {i, 
          1, n}]; RR = Table[RR[[i]] = R[[i]], {i, 1, q - 1}]; 
      b = Max[RR]; 
      Do[If[x[[i]] == b, bi = i], {i, 1, n}]]

firsttwo := {ci, bi}

If more than two values are desired, the following selects the subsequent values and puts all the values into a matrix, maintaining the original data order.

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

Do[A[[i, 2]] = i; 
      A[[i, 1]] = comp[[i]], {i, 1, n}]

A = Sort[A];

Do[A[[i, 3]] = i, {i, 1, n}]

T = A[[All, 1]];

Do[A[[i, 1]] = A[[i, 2]]; 
      A[[i, 2]] = T[[i]], {i, 1, n}]

A = Sort[A];

d = A[[All, 3]];

If[d[[firsttwo[[2]]]] ≠ 2, 
      temp = d[[firsttwo[[2]]]]; d[[firsttwo[[2]]]] = 1; 
      Do[If[i ≠ firsttwo[[2]]  &&  
            i ≠ firsttwo[[1]]  &&  
            d[[i]]  ≤  temp, d[[i]] = d[[i]] + 1]; 
        d[[firsttwo[[1]]]] = 1, {i, 1, n}]]

Linear Interpolation

Pick two data points

datapoints = 2;

xData = Table[0, {i, 1, datapoints}];
yData = Table[0, {i, 1, datapoints}];

p = 1; 
    Do[If[d[[i]] ≤ datapoints, xData[[p]] = x[[i]]; 
        yData[[p]] = y[[i]]; p = p + 1], {i, 1, n}]

xData  // MatrixForm

yData  // MatrixForm

M = {{1, xData[[1]]}, {1, xData[[2]]}} // MatrixForm

Co  =  LinearSolve[%, yData]

f[x_] := Co[[1]] + Co[[2]]*x

f[xdesired]

fprev = %;

lin = Plot[f[x], {x, Min[xData], Max[xData]}];

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

Show[desire, lin];

Show[data, lin, desire, PlotRange -> All];

Quadratic Interpolation

Pick two data points

datapoints = 3;

xData = Table[0, {i, 1, datapoints}];
yData = Table[0, {i, 1, datapoints}];

p = 1; 
    Do[If[d[[i]] ≤ datapoints, xData[[p]] = x[[i]]; 
        yData[[p]] = y[[i]]; p = p + 1], {i, 1, n}]

xData  // MatrixForm

yData  // MatrixForm

M = {{1, xData[[1]], xData[[1]]^2}, {1, xData[[2]], 
            xData[[2]]^2}, {1, xData[[3]], xData[[3]]^2}} // 
        MatrixForm

Co  =  LinearSolve[%, yData]

f[x_] := Co[[1]] + Co[[2]]*x + Co[[3]]*x^2

f[xdesired]

fnew = %;

ϵa = Abs[(fnew - fprev)/fnew*100]

sigdig = Floor[2 - Log[10, (ϵa/0.5)]]

fprev = %%%;

lin = Plot[f[x], {x, Min[xData], Max[xData]}];

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

Show[desire, lin];

Show[data, lin, desire, PlotRange -> All];

Cubic Interpolation

Pick four data points

datapoints = 4;

xData = Table[0, {i, 1, datapoints}];
yData = Table[0, {i, 1, datapoints}];

d

p = 1; 
    Do[If[d[[i]] ≤ datapoints, xData[[p]] = x[[i]]; 
        yData[[p]] = y[[i]]; p = p + 1], {i, 1, n}]

xData  // MatrixForm

yData  // MatrixForm

M = {{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}} // 
        MatrixForm

Co  =  LinearSolve[%, yData]

f[x_] := 
      Co[[1]] + Co[[2]]*x + Co[[3]]*x^2 + 
        Co[[4]]*x^3

f[xdesired]

fnew = %;

ϵa = Abs[(fnew - fprev)/fnew*100]

sigdig = Floor[2 - Log[10, (ϵa/0.5)]]

lin = Plot[f[x], {x, Min[xData], Max[xData]}];

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

Show[desire, lin];

Show[data, lin, desire, PlotRange -> All];