Open Educational Resources

Introduction to Numerical Analysis: Piecewise Interpolation

Introduction

In the previous section, it was shown that when the order of the interpolating polynomial increases — which is natural when there is a large number of data points — the interpolating polynomial function can highly oscillate or fluctuate between the data points. With a large number of data points, it is better to keep the order of interpolation low. This can be done using piecewise interpolation using splines. Splines are polynomials on subintervals that are connected in a smooth manner. The order of the polynomials is not associated with the number of data points in this case, so, we can have k data points, and use a spline of degree n where k>>n. Assuming k data points (x_1,y_1), (x_2,y_2),\cdots,(x_k,y_k) such that x_1<x_2<\cdots<x_k, a spline function y of degree n satisfies the following:

  • On each subinterval \left[ x_i,x_{i+1} \right), y is a polynomial of degree n
  • The derivatives of y up to and including the n-1 derivative are continuous on the domain [x_1,x_k].

The data points, when chosen such that the spline function passes through them, are called knots. For example, a spline of degree 0 is a piecewise constant function constructed by simply assuming y to have a constant value on any subinterval \left[x_i,x_{i+1}\right) (Figure 1). A spline of degree 1 is a piecewise linear (piecewise affine) function constructed by connecting a straight line between every two data points (Figure 1). A spline of degree 2 is a piecewise quadratic function on each subinterval. The first derivative exists and is continuous, but it is not smooth. As shown in Figure 2, the first derivative is a broken line at the knots and thus the derivatives at the knots do not exist. A spline of degree 3 is a piecewise cubic function on each subinterval. The first and second derivatives are continuous functions. The second derivative is not smooth. Figure 1 shows the resulting 0 order, linear, quadratic, and cubic spline interpolating functions for the data points generated by the Runge function described in the previous section. Figure 2 on the other hand shows the generated derivatives. The 0 order spline function gives a zero derivative every where. The linear spline function has derivatives that are constant on each subinterval. The quadratic spline gives derivatives that are not smooth at the data points. The cubic spline provides smooth derivatives that are very close to the actual derivative of the Runge function that was used to generate the data!

The chosen order of the spline function depends on the accuracy sought. Linear splines can provide acceptable interpolation between the knots. However, their derivatives are constant. Linear splines are very useful when the data points are very close to each other, and there is no interest in their derivatives. Quadratic and cubic splines provide very smooth interpolation of the data points. However, the derivatives of the cubic splines are “smoother” than the quadratic splines. The choice of the spline order depends on the oscillatory nature of the data points to be interpolated.

The “Interpolation” built-in function generates an interpolation function through the data. To choose splines, the option: Method->”Spline” has to be specified as well as the interpolation order. The following code generates the 0 order, linear, quadratic, and cubic spline functions through the data points of the Runge function shown in the previous section. The code also compares the first derivatives of the spline functions with the actual first derivative of the Runge function.

View Mathematica Code
Data = {{-1, 0.038}, {-0.8, 0.058}, {-0.60, 0.10}, {-0.4, 0.20}, {-0.2, 0.50}, {0, 1}, {0.2, 0.5}, {0.4, 0.2}, {0.6, 0.1}, {0.8, 0.058}, {1, 0.038}};
y0 = Interpolation[Data, Method -> "Spline", InterpolationOrder -> 0];
y1 = Interpolation[Data, Method -> "Spline", InterpolationOrder -> 1];
y2 = Interpolation[Data, Method -> "Spline", InterpolationOrder -> 2];
y3 = Interpolation[Data, Method -> "Spline", InterpolationOrder -> 3];
yactual = 1/(1 + 25 x^2);
a = Plot[{y0[x], yactual}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y0"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y0", "yactual"}, PlotLabel -> "0 order splines"];
b = Plot[{y1[x], yactual}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]},AxesLabel -> {"x", "y1"}, ImageSize -> Medium, PlotRange -> All,PlotLegends -> {"y1", "yactual"}, PlotLabel -> "Linear splines"];
c = Plot[{y2[x], yactual}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]},AxesLabel -> {"x", "y2"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y2", "yactual"}, PlotLabel -> "Quadratic splines"];
d = Plot[{y3[x], yactual}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]},AxesLabel -> {"x", "y3"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y3", "yactual"}, PlotLabel -> "Cubic splines"];
Grid[{{a, b}, {c, d}}]
y0prime = D[y0[x], x];
y1prime = D[y1[x], x];
y2prime = D[y2[x], x];
y3prime = D[y3[x], x];
yactualprime = D[yactual, x];
a = Plot[{y0prime, yactualprime}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y0"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y0'", "y'actual"}, PlotLabel -> "Derivative of 0 order splines"];
b = Plot[{y1prime, yactualprime}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]},AxesLabel -> {"x", "y1"}, ImageSize -> Medium, PlotRange -> All,PlotLegends -> {"y1'", "y'actual"}, PlotLabel -> "Linear splines"];
c = Plot[{y2prime, yactualprime}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]},AxesLabel -> {"x", "y2"}, ImageSize -> Medium, PlotRange -> All,PlotLegends -> {"y2'", "y'actual"}, PlotLabel -> "Quadratic splines"];
d = Plot[{y3prime, yactualprime}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y3"}, ImageSize -> Medium, PlotRange -> All,PlotLegends -> {"y3'", "y'actual"}, PlotLabel -> "Cubic splines"];
Grid[{{a, b}, {c, d}}]


Figure 1. Piecewise interpolation using: 0 order, linear, quadratic, and cubic splines

Figure 1. Piecewise interpolation using: 0 order, linear, quadratic, and cubic splines

Figure 2. The derivatives of piecewise interpolation using 0 order, linear, quadratic, and cubic splines. The derivative of the 0 order splines is the zero function.

Figure 2. The derivatives of piecewise interpolation using 0 order, linear, quadratic, and cubic splines. The derivative of the 0 order splines is the zero function.

0 Degree Spline Interpolation

Given a set of k data points (x_1,y_1),(x_2,y_2),\cdots, (x_k,y_k), a spline of degree 0 can be defined as:

    \[ y=\begin{cases} y_1,&x\in [x_1, x_2) \\ y_2,&x\in [x_2, x_3) \\ \vdots & \vdots \\ y_{k-1},&x\in [x_{k-1}, x_k) \\ y_k, & x=x_k \end{cases} \]

In the previous spline definition, the interpolated value in the interval [x_{i-1},x_i) is equal to the data value y_{i-1}. Alternatively, the interpolated value in the interval (x_{i-1},x_i] can be set equal to the data value y_i. In this case, y can be defined as follows:

    \[ y=\begin{cases} y_1,&x=x_1\\ y_2,&x\in (x_1, x_2] \\ \vdots \vdots \\ y_k,&x\in (x_{k-1}, x_k] \end{cases} \]

Piecewise Linear (Piecewise Affine) Spline Interpolation

Given a set of k data points (x_1,y_1),(x_2,y_2),\cdots, (x_k,y_k), a piecewise linear (piecewise affine) spline can be defined as:

    \[ y=\begin{cases} a_1x+b_1,&x\in [x_1, x_2) \\ a_2x+b_2,&x\in [x_2, x_3) \\ \vdots & \vdots \\ a_{k-1}x+b_{k-1},&x\in [x_{k-1}, x_k] \end{cases} \]

The k data points have k-1 intervals. The linear function for each interval is defined using two coefficients, and therefore, we need to find 2(k-1) coefficients a_1,b_1,a_2,b_2,\cdots,a_{k-1},b_{k-1}. The coefficients a_i and b_i can be found by solving the two equations:

    \[\begin{split} y_i&=a_ix_i + b_i\\ y_{i+1}&=a_ix_{i+1}+b_i \end{split} \]

Therefore:

    \[\begin{split} a_i&=\frac{y_i-y_{i+1}}{x_i-x_{i+1}}\\ b_i&=\frac{y_{i+1}x_i-y_ix_{i+1}}{x_i-x_{i+1}} \end{split} \]

Example

Consider the following 11 data points: (-1.00,0.038),(-0.80,0.058),(-0.60,0.100),(-0.4,0.200),(-0.20,0.500),(0.00,1.000),(0.20,0.500),(0.40,0.200),(0.60,0.100),(0.80,0.0580),(1.00,0.038) which were generated using the Runge function. Find the explicit representation of the linear spline interpolating function for these data points.

Solution

There are 11 data points surrounding 10 intervals. For each interval i, a set of coefficients a_i and b_i are required for the linear interpolation. For the first interval, the coefficients are:

    \[\begin{split} a_1&=\frac{y_1-y_2}{x_1-x_2}=\frac{0.038-0.058}{-1-(-0.8)}=0.1\\ b_1&=\frac{y_2x_1-y_1x_2}{x_1-x_2}=\frac{0.058(-1)-0.038(-0.8)}{-1-(-0.8)}=0.138 \end{split} \]

Similarly:

    \[\begin{split} a_2&=\frac{y_2-y_3}{x_2-x_3}=\frac{0.058-0.1}{-0.8-(-0.6)}=0.21\\ b_2&=\frac{y_3x_2-y_2x_3}{x_2-x_3}=\frac{0.1(-0.8)-0.058(-0.6)}{-0.8-(-0.6)}=0.226 \end{split} \]

Repeating for the other coefficients, the required explicit equation has the form:

    \[ y=\begin{cases} 0.1x+0.138,&x\in [-1, -0.8) \\ 0.21x+0.226,&x\in [-0.8, -0.6) \\ 0.5x+0.4,&x\in [-0.6, -0.4)\\ 1.5x+0.8,&x\in [-0.4, -0.2)\\2.5x+1,&x\in [-0.2, 0.0)\\-2.5x+1,&x\in [0.0, 0.2)\\-1.5x+0.8,&x\in [0.2, 0.4)\\-0.5x+0.4,&x\in [0.4, 0.6)\\-0.21x+0.226,&x\in [0.6, 0.8)\\-0.1x+0.138,&x\in [0.8, 1.0] \end{cases} \]

Notice that for the procedure to work, the x components of the data points have to satisfy: x_1<x_2<\cdots<x_{k}. The following Mathematica code utilizes a user defined procedure “Spline1” that creates the required piecewise linear function.

View Mathematica Code
Data = {{-1, 0.038}, {-0.8, 0.058}, {-0.60, 0.10}, {-0.4,0.20}, {-0.2, 0.50}, {0, 1}, {0.2, 0.5}, {0.4, 0.2}, {0.6,0.1}, {0.8, 0.058}, {1, 0.038}};
Spline1[Data_] := (
k1 = Length[Data] - 1;  
atable = Table[(Data[[i, 2]] - Data[[i + 1, 2]])/(Data[[i, 1]] -Data[[i + 1, 1]]), {i, 1, k1}]; 
btable = Table[(Data[[i + 1, 2]]*Data[[i, 1]] -Data[[i, 2]]*Data[[i + 1, 1]])/(Data[[i, 1]] - Data[[i + 1, 1]]), {i, 1, k1}]; 
pf = Table[{atable[[i]] x + btable[[i]],Data[[i, 1]] <= x <= Data[[i + 1, 1]]}, {i, 1, k1}]; 
Piecewise[pf]
)
y = Spline1[Data]
a = Plot[{y, yactual}, {x, -1, 1}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y1"}, ImageSize -> Medium, PlotRange -> All,  PlotLegends -> {"y1", "yactual"}, PlotLabel -> "Linear splines"]

Quadratic Spline Interpolation

For quadratic spline interpolation, we present two possible quadratic interpolation schemes.

Scheme 1:

In the first scheme, the intervals between the data points are used as intervals on which a quadratic function is defined. On the intermediate nodes, the slope of the parabola on the left is equated to the slope of the parabola on the right (Figure 3). Consider a set of k data points (x_1,y_1),(x_2,y_2),\cdots, (x_k,y_k), with k-1 intervals. We seek to find a quadratic polynomial s_i(x)=a_i+b_ix+c_ix^2 on the i^{th} interval [x_i,x_{i+1}]. In total, there are 3(k-1) coefficients. On each interval, the two equations s_i(x_i)=y_i and s_i(x_{i+1})=y_{i+1} ensure that the spline passes through the data points which also ensure the continuity of the splines. These continuity equations will provide 2(k-1) equations. In addition, the connection of the quadratic polynomials at intermediate nodes will be assumed to be smooth, and so at every intermediate point the first derivative is assumed to be continuous, that is s_i'(x_{i+1})=s_{i+1}'(x_{i+1}). These first-derivative continuity equations provide another k-2 equations. Finally, the second derivative at the starting point can be set to 0 to provide one additional condition so that the number of equations would be equal to 3(k-1) which is equal to the number of unknowns. The following are the resulting equations:

    \[\begin{split} a_1+b_1x_1+c_1x_1^2&=y_1\\ a_1+b_1x_2+c_1x_2^2&=y_2\\ a_2+b_2x_2+c_2x_2^2&=y_2\\ a_2+b_2x_3+c_2x_3^2&=y_3\\ \vdots&\\ a_{k-1}+b_{k-1}x_{k-1}+c_{k-1}x_{k-1}^2&=y_{k-1}\\ a_{k-1}+b_{k-1}x_k+c_{k-1}x_k^2&=y_k\\ b_1+2c_1x_2-b_2-2c_2x_2&=0\\ b_2+2c_2x_3-b_3-2c_3x_3&=0\\ \vdots&\\ b_{k-2}+2c_{k-2}x_{k-1}-b_{k-1}-2c_{k-1}x_{k-1}&=0\\ 2c_1&=0 \end{split} \]

Figure 3. Scheme 1 of piecewise quadratic interpolation with continuity in the interpolating function and its first derivative

Figure 3. Scheme 1 of piecewise quadratic interpolation with continuity in the interpolating function and its first derivative

The following Mathematica code implements this procedure for a set of data (the procedure will only work if x_1<x_2<\cdots<x_k).

View Mathematica Code
Data = {{-1, 0.038}, {-0.8, 0.058}, {-0.60, 0.10}, {-0.4, 0.20}, {-0.2, 0.50}, {0, 1}, {0.2, 0.5}, {0.4, 0.2}, {0.6, 0.1}, {0.8, 0.058}, {1, 0.038}};
Spline2[Data_] := (k1 = Length[Data] - 1;
  y = Table[If[i <= 2 k1,If[EvenQ[i], Data[[i/2 + 1, 2]], Data[[(i + 1)/2, 2]]], 0], {i,1, 3 k1}];
  M = Table[0, {i, 1, 3 k1}, {j, 1, 3 k1}];
  For[i = 1, i <= k1, i = i + 1, 
   M[[2 i - 1, 3 (i - 1) + 1]] = M[[2 i, 3 (i - 1) + 1]] = 1; 
   M[[2 i - 1, 3 (i - 1) + 2]] = Data[[i, 1]]; 
   M[[2 i - 1, 3 (i - 1) + 1 + 2]] = Data[[i, 1]]^2; 
   M[[2 i, 3 (i - 1) + 1 + 1]] = Data[[i + 1, 1]]; 
   M[[2 i, 3 (i - 1) + 1 + 2]] = Data[[i + 1, 1]]^2];
  For[i = 1, i <= k1 - 1, i = i + 1, 
   M[[2 k1 + i, 3 (i - 1) + 1 + 1]] = 1; 
   M[[2 k1 + i, 3 (i - 1) + 1 + 4]] = -1; 
   M[[2 k1 + i, 3 (i - 1) + 1 + 2]] = 2 Data[[i + 1, 1]]; 
   M[[2 k1 + i, 3 (i - 1) + 1 + 5]] = -2 Data[[i + 1, 1]]];
  M[[3 k1, 3]] = 2;
  Coef = Inverse[M].y;
  pf = Table[{Coef[[3 (i - 1) + 1]] + Coef[[3 (i - 1) + 2]]*x + Coef[[3 (i - 1) + 3]]*x^2, Data[[i, 1]] <= x <= Data[[i + 1, 1]]}, {i, 1, k1}];
  Piecewise[pf])
y = Spline2[Data]
M // MatrixForm
yactual=1/(1+25x^2);
a = Plot[{y, yactual}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y1"},ImageSize -> Medium, PlotRange -> All,PlotLegends -> {"y2", "yactual"}, PlotLabel ->"Quadratic spline"]
Example

Use scheme 1 of quadratic interpolation to find the interpolating function given the following 4 data points (-1, 0.038) (-0.8, 0.058), (-0.60, 0.10), (-0.4, 0.20)

Solution

The four data points (k=4) with three intervals (k-1=3), therefore, there are 3(k-1)=9 unknowns. The following are the nine equations according to scheme 1 to be solved to find the unknowns:

    \[\begin{split} a_1+b_1(-1)+c_1(-1)^2&=0.038\\ a_1+b_1(-0.8)+c_1(-0.8)^2&=0.058\\ a_2+b_2(-0.8)+c_2(-0.8)^2&=0.058\\ a_2+b_2(-0.6)+c_2(-0.6)^2&=0.10\\ a_3+b_3(-0.6)+c_3(-0.6)^2&=0.10\\ a_3+b_3(-0.4)+c_3(-0.4)^2&=0.20\\ b_1+2c_1(-0.8)-b_2-2c_2(-0.8)&=0\\ b_2+2c_2(-0.6)-b_3-2c_3(-0.6)&=0\\ 2c_1&=0 \end{split} \]

Solving the above equations produces the following interpolation function:

    \[ y=\begin{cases}s_1(x)=0.138+0.1x,&-1\leq x<-0.8\\ s_2(x)=0.49+0.98x+0.55x^2, &-0.8 \leq x<-0.6\\ s_3(x)=0.616 + 1.4 x + 0.9 x^2, & -0.6 \leq x \leq -0.4 \end{cases} \]

Figure 4 shows the resulting interpolation function along with the four data points. The Mathematica code is provided below the figure.

Figure 4. Scheme 1 quadratic interpolation applied to four data points

Figure 4. Scheme 1 quadratic interpolation applied to four data points

View Mathematica Code
Data = {{-1, 0.038}, {-0.8, 0.058}, {-0.60, 0.10}, {-0.4, 0.20}};
Spline2[Data_] := (k1 = Length[Data] - 1;
  y = Table[
    If[i <= 2 k1, 
     If[EvenQ[i], Data[[i/2 + 1, 2]], Data[[(i + 1)/2, 2]]], 0], {i, 
     1, 3 k1}];
  M = Table[0, {i, 1, 3 k1}, {j, 1, 3 k1}];
  For[i = 1, i <= k1, i = i + 1, 
   M[[2 i - 1, 3 (i - 1) + 1]] = M[[2 i, 3 (i - 1) + 1]] = 1;
   M[[2 i - 1, 3 (i - 1) + 2]] = Data[[i, 1]];
   M[[2 i - 1, 3 (i - 1) + 1 + 2]] = Data[[i, 1]]^2;
   M[[2 i, 3 (i - 1) + 1 + 1]] = Data[[i + 1, 1]];
   M[[2 i, 3 (i - 1) + 1 + 2]] = Data[[i + 1, 1]]^2];
  For[i = 1, i <= k1 - 1, i = i + 1, 
   M[[2 k1 + i, 3 (i - 1) + 1 + 1]] = 1;
   M[[2 k1 + i, 3 (i - 1) + 1 + 4]] = -1;
   M[[2 k1 + i, 3 (i - 1) + 1 + 2]] = 2 Data[[i + 1, 1]];
   M[[2 k1 + i, 3 (i - 1) + 1 + 5]] = -2 Data[[i + 1, 1]]];
  M[[3 k1, 3]] = 2;
  Coef = Inverse[M].y;
  pf = Table[{Coef[[3 (i - 1) + 1]] + Coef[[3 (i - 1) + 2]]*x + Coef[[3 (i - 1) + 3]]*x^2, Data[[i, 1]] <= x <= Data[[i + 1, 1]]}, {i, 1, k1}];
  Piecewise[pf])
y = Spline2[Data]
M // MatrixForm
a = Plot[{y}, {x, -1, -0.4}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y2"},ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y2"}, PlotLabel -> "Quadratic spline"]
Disadvantages

In certain cases, such as fluctuating data, this scheme might lead to some oscillations within the intervals between the data points. As an example, when this scheme is used to fit through the 11 data points of the Runge function given in the previous section, a nice smooth interpolating function is produced for the first few intervals (Figure 5). The condition of having zero second derivative s_1''(x_1)=0 on the initial interval causes the interpolating in the first interval to follow the line connecting the first two data points. However, in the last five intervals, the interpolating function oscillates and deviates away from the expected behaviour.


Figure 5. Behaviour of scheme 1 for quadratic interpolation of the Runge function data points.

Figure 5. Behaviour of scheme 1 for quadratic interpolation of the Runge function data points.

Scheme 2

For the second scheme, the intervals where the quadratic interpolation functions are defined, are chosen different from the intervals between the data points in a manner to ensure symmetry of the quadratic interpolation. The new intervals where the piecewise quadratic functions are defined are bounded by the points x_1<\frac{x_2+x_3}{2}<\frac{x_3+x_4}{2}<\cdots<\frac{x_{k-2}+x_{k-1}}{2}<x_k (Figure 6). Consider a set of k data points (x_1,y_1),(x_2,y_2),\cdots, (x_k,y_k), with k-1 intervals. On each of the k-2 intervals bounded by x_1<\frac{x_2+x_3}{2}<\frac{x_3+x_4}{2}<\cdots<\frac{x_{k-2}+x_{k-1}}{2}<x_k a quadratic polynomial s_i(x)=a_i+b_ix+c_ix^2 is defined. In total, there are 3(k-2) coefficients. First,the quadratic polynomials have to pass through the k data poitns resulting in k equations. These quadratic polynomials have to be continuous and differentiable at the k-3 intermediate points that are the bounds of the intervals resulting in 2(k-3) equations. Therefore, in total, there are 3(k-2) equations which is equal to the number of unknowns. The following are the resulting equations:

    \[\begin{split} a_1+b_1x_1+c_1x_1^2&=y_1\\ a_1+b_1x_2+c_1x_2^2&=y_2\\ a_2+b_2x_3+c_2x_3^2&=y_3\\a_3+b_3x_4+c_3x_4^2&=y_4\\ \vdots&\\ a_{k-2}+b_{k-2}x_{k-1}+c_{k-2}x_{k-1}^2&=y_{k-1}\\ a_{k-2}+b_{k-2}x_k+c_{k-2}x_k^2&=y_k\\ \vdots&\\ a_1+b_1\left(\frac{x_2+x_3}{2}\right)+c_1\left(\frac{x_2+x_3}{2}\right)^2&=a_2+b_2\left(\frac{x_2+x_3}{2}\right)+c_2\left(\frac{x_2+x_3}{2}\right)^2\\ a_2+b_2\left(\frac{x_3+x_4}{2}\right)+c_2\left(\frac{x_3+x_4}{2}\right)^2&=a_3+b_3\left(\frac{x_3+x_4}{2}\right)+c_3\left(\frac{x_3+x_4}{2}\right)^2\\ a_{k-3}+b_{k-3}\left(\frac{x_{k-2}+x_{k-1}}{2}\right)+c_{k-3}\left(\frac{x_{k-2}+x_{k-1}}{2}\right)^2&=\\ a_{k-2}+b_{k-2}\left(\frac{x_{k-2}+x_{k-1}}{2}\right)&+c_{k-2}\left(\frac{x_{k-2}+x_{k-1}}{2}\right)^2\\ b_1+2c_1\left(\frac{x_2+x_3}{2}\right)-b_2-2c_2\left(\frac{x_2+x_3}{2}\right)&=0\\ b_2+2c_2\left(\frac{x_3+x_4}{2}\right)-b_3-2c_3\left(\frac{x_3+x_4}{2}\right)&=0\\ \vdots&\\ b_{k-3}+2c_{k-3}\left(\frac{x_{k-2}+x_{k-1}}{2}\right)-b_{k-2}-2c_{k-2}&\left(\frac{x_{k-2}+x_{k-1}}{2}\right)=0 \end{split} \]

Scheme 2 of piecewise quadratic interpolation with continuity in the interpolating function and its first derivative

Figure 6. Scheme 2 of piecewise quadratic interpolation with continuity in the interpolating function and its first derivative

The following Mathematica code implements this procedure for a set of data (the procedure will only work if x_1<x_2<\cdots<x_k). As shown in Figure 7, the oscillations produced by scheme 1 disappear if this scheme is used for interpolating the data points of the Runge function given in the previous section.

View Mathematica Code
Data = {{-1, 0.038}, {-0.8, 0.058}, {-0.60, 0.10}, {-0.4,0.20}, {-0.2, 0.50}, {0, 1}, {0.2, 0.5}, {0.4, 0.2}, {0.6,0.1}, {0.8, 0.058}, {1, 0.038}};
Spline2[Data_] := (k = Length[Data];
  NewData = Table[0, {i, 1, k - 3}];
  Do[NewData[[i]] = (Data[[i + 1]] + Data[[i + 2]])/2, {i, 1, k - 3}];
  NNewData = Prepend[Append[NewData, Data[[k]]], Data[[1]]];
  y = Table[If[i <= k - 1 + 1, Data[[i, 2]], 0], {i, 1, 3 k - 6}];
  M = Table[0, {i, 1, 3 k - 6}, {j, 1, 3 k - 6}];
  M[[1, 1]] = M[[2, 1]] = 1;
  M[[1, 2]] = Data[[1, 1]];
  M[[1, 3]] = Data[[1, 1]]^2;
  M[[2, 2]] = Data[[2, 1]];
  M[[2, 3]] = Data[[2, 1]]^2;
  M[[k - 1, 3 k - 8]] = M[[k, 3 k - 8]] = 1;
  M[[k - 1, 3 k - 7]] = Data[[k - 1, 1]];
  M[[k - 1, 3 k - 6]] = Data[[k - 1, 1]]^2;
  M[[k, 3 k - 7]] = Data[[k, 1]];
  M[[k, 3 k - 6]] = Data[[k, 1]]^2;
  For[i = 3, i <= k - 2, i = i + 1, M[[i, 3 (i - 2) + 1]] = 1; 
   M[[i, 3 (i - 2) + 2]] = Data[[i, 1]]; 
   M[[i, 3 (i - 2) + 3]] = Data[[i, 1]]^2];
  For[i = 1, i <= k - 3, i = i + 1,
   M[[k + i, 3 (i - 1) + 1]] = 1; M[[k + i, 3 (i - 1) + 4]] = -1; 
   M[[k + i, 3 (i - 1) + 2]] = NewData[[i, 1]]; 
   M[[k + i, 3 (i - 1) + 5]] = -NewData[[i, 1]]; 
   M[[k + i, 3 (i - 1) + 3]] = NewData[[i, 1]]^2; 
   M[[k + i, 3 (i - 1) + 6]] = -1*NewData[[i, 1]]^2;
   M[[k + i + k - 3, 3 (i - 1) + 2]] = 1; 
   M[[k + i + k - 3, 3 (i - 1) + 5]] = -1; 
   M[[k + i + k - 3, 3 (i - 1) + 3]] = 2*NewData[[i, 1]]; 
   M[[k + i + k - 3, 3 (i - 1) + 6]] = -2*NewData[[i, 1]]];
  Coef = Inverse[M].y;
  pf = Table[{Coef[[3 (i - 1) + 1]] + Coef[[3 (i - 1) + 2]]*x + Coef[[3 (i - 1) + 3]]*x^2,NNewData[[i, 1]] <= x <= NNewData[[i + 1, 1]]}, {i, 1, k - 2}];
  Piecewise[pf])
y = Spline2[Data]
yactual = 1/(1 + 25 x^2);
M//MatrixForm
a = Plot[{y, yactual}, {x, -1, 1}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y2"},ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y2", "yactual"}, PlotLabel -> "Quadratic spline"]


Figure 7. Behaviour of scheme 2 of quadratic interpolation when applied to the Runge function data points

Figure 7. Behaviour of scheme 2 of quadratic interpolation when applied to the Runge function data points

Example

Use scheme 2 of quadratic interpolation to find the interpolating function given the following 5 data points (-1, 0.038) (-0.8, 0.058), (-0.60, 0.10), (-0.4, 0.20), (-0.2, 0.5)

Solution

The five data points (k=5) with four intervals (k-1=4). k-2=3 intervals for the quadratic interpolation functions are defined between the points -1<-0.7<-0.5<-0.2. Therefore, there are 3(k-2)=9 unknowns. The following are the nine equations according to scheme 2 to be solved to find the unknowns:

    \[\begin{split} a_1+b_1(-1)+c_1(-1)^2&=0.038\\ a_1+b_1(-0.8)+c_1(-0.8)^2&=0.058\\ a_2+b_2(-0.6)+c_2(-0.6)^2&=0.10\\ a_3+b_3(-0.4)+c_3(-0.4)^2&=0.20\\ a_3+b_3(-0.2)+c_3(-0.2)^2&=0.5\\ a_1+b_1(-0.7)+c_1(-0.7)^2&=a_2+b_2(-0.7)+c_2(-0.7)^2\\ a_2+b_2(-0.5)+c_2(-0.5)^2&=a_3+b_3(-0.5)+c_3(-0.5)^2\\ b_1+2c_1(-0.7)-b_2-2c_2(-0.7)&=0\\ b_2+2c_2(-0.5)-b_3-2c_3(-0.5)&=0\\ \end{split} \]

Solving the above equations produces the following interpolation function:

    \[ y=\begin{cases}s_1(x)=0.336857 + 0.547429 x + 0.248571 x^2,&-1\leq x<-0.7\\ s_2(x)=0.440457 + 0.843429 x + 0.46 x^2, &-0.7 \leq x<-0.5\\ s_3(x)=1.02331 + 3.17486 x + 2.79143 x^2, & -0.5\leq x\leq -0.2 \end{cases} \]

Figure 8 shows the resulting interpolation function along with the five data points. The Mathematica code is provided below the figure.


Figure 8. Scheme 1 quadratic interpolation applied to five data points

Figure 8. Scheme 1 quadratic interpolation applied to five data points

View Mathematica Code
Data = {{-1, 0.038}, {-0.8, 0.058}, {-0.60, 0.10}, {-0.4,0.20}, {-0.2, 0.50}};
Spline2[Data_] := (k = Length[Data];
  NewData = Table[0, {i, 1, k - 3}];
  Do[NewData[[i]] = (Data[[i + 1]] + Data[[i + 2]])/2, {i, 1, k - 3}];
  NNewData = Prepend[Append[NewData, Data[[k]]], Data[[1]]];
  y = Table[If[i <= k - 1 + 1, Data[[i, 2]], 0], {i, 1, 3 k - 6}];
  M = Table[0, {i, 1, 3 k - 6}, {j, 1, 3 k - 6}];
  M[[1, 1]] = M[[2, 1]] = 1;
  M[[1, 2]] = Data[[1, 1]];
  M[[1, 3]] = Data[[1, 1]]^2;
  M[[2, 2]] = Data[[2, 1]];
  M[[2, 3]] = Data[[2, 1]]^2;
  M[[k - 1, 3 k - 8]] = M[[k, 3 k - 8]] = 1;
  M[[k - 1, 3 k - 7]] = Data[[k - 1, 1]];
  M[[k - 1, 3 k - 6]] = Data[[k - 1, 1]]^2;
  M[[k, 3 k - 7]] = Data[[k, 1]];
  M[[k, 3 k - 6]] = Data[[k, 1]]^2;
  For[i = 3, i <= k - 2, i = i + 1, M[[i, 3 (i - 2) + 1]] = 1;
   M[[i, 3 (i - 2) + 2]] = Data[[i, 1]];
   M[[i, 3 (i - 2) + 3]] = Data[[i, 1]]^2];
  For[i = 1, i <= k - 3, i = i + 1, M[[k + i, 3 (i - 1) + 1]] = 1; 
   M[[k + i, 3 (i - 1) + 4]] = -1;
   M[[k + i, 3 (i - 1) + 2]] = NewData[[i, 1]];
   M[[k + i, 3 (i - 1) + 5]] = -NewData[[i, 1]];
   M[[k + i, 3 (i - 1) + 3]] = NewData[[i, 1]]^2;
   M[[k + i, 3 (i - 1) + 6]] = -1*NewData[[i, 1]]^2;
   M[[k + i + k - 3, 3 (i - 1) + 2]] = 1;
   M[[k + i + k - 3, 3 (i - 1) + 5]] = -1;
   M[[k + i + k - 3, 3 (i - 1) + 3]] = 2*NewData[[i, 1]];
   M[[k + i + k - 3, 3 (i - 1) + 6]] = -2*NewData[[i, 1]]];
  Coef = Inverse[M].y;
  pf = Table[{Coef[[3 (i - 1) + 1]] + Coef[[3 (i - 1) + 2]]*x + Coef[[3 (i - 1) + 3]]*x^2, NNewData[[i, 1]] <= x <= NNewData[[i + 1, 1]]}, {i, 1, k - 2}];
  Piecewise[pf])
y = Spline2[Data]
M // MatrixForm
a = Plot[{y}, {x, -1, -0.2}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y2"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y2"},  PlotLabel -> "Quadratic spline"]

Cubic Spline Interpolation

There are different schemes of piecewise cubic spline interpolation functions which vary according to the end conditions. The scheme presented here is sometimes referred to as “Not-a-knot” end condition in which the first cubic spline is defined over the interval x\in [x_1,x_3) and the last cubic spline is defined on the interval x\in[x_{k-2},x_k] (Figure 9). Given k data points with k-1 intervals, k-3 cubic polynomials are defined on the intervals: x_1<x_3<x_4<\cdots<x_{k-3}<x_{k-2}<x_k. \forall i \leq k-3 each cubic polynomial is defined as: s_i(x)=a_i+b_ix+c_ix^2+d_ix^3. The 4(k-3) unknowns can be obtained using the following conditions(Figure 9):

  • Each polynomial has to pass through the data points which are the bounds of the interval, therefore this results in 2(k-3) equations. In addition the first polynomial and the last (k-3) polynomial have to pass through the second and next-to-last data points respectively. In total, this produces 2(k-3)+2 equations.
  • On the k-4 intermediate nodes between the polynomials, the first and second derivatives are equal to each other. This results in 2(k-4) equations

The above system of equations is guaranteed a unique solution as long as the data points are ordered: x_1<x_2<\cdots<x_k.



>

Figure 9. Scheme of piecewise cubic interpolation with continuity in the first and second derivatives.

Figure 9. Scheme of piecewise cubic interpolation with continuity in the first and second derivatives.

The following Mathematica code implements this procedure for a set of data (the procedure will only work if x_1<x_2<\cdots<x_k).
As shown in Figure 10, this cubic interpolation scheme produces a very smooth interpolating function for the data points of the Runge function given in the previous section.

View Mathematica Code
Data = {{-1, 0.038}, {-0.8, 0.058}, {-0.60, 0.10}, {-0.4,0.20}, {-0.2, 0.50}, {0, 1}, {0.2, 0.5}, {0.4, 0.2}, {0.6, 0.1}, {0.8, 0.058}, {1, 0.038}};
Clear[x, a, b, c, d]
Spline3[Data_] := (k = Length[Data];
  atable = Table[Subscript[a, i], {i, 1, k - 3}];
  btable = Table[Subscript[b, i], {i, 1, k - 3}];
  ctable = Table[Subscript[c, i], {i, 1, k - 3}];
  dtable = Table[Subscript[d, i], {i, 1, k - 3}];
  Poly = Table[atable[[i]] + btable[[i]]*x + ctable[[i]]*x^2 + dtable[[i]]*x^3, {i, 1, k - 3}];
  intermediate = Table[Data[[i + 2]], {i, 1, k - 4}];
  Equations = Table[0, {i, 1, 4 (k - 3)}];
  Do[Equations[[i]] = (Poly[[i]] /. x -> Data[[i + 2, 1]]) == Data[[i + 2, 2]], {i, 1, k - 4}];
  Do[Equations[[i + k - 4]] = (Poly[[i + 1]] /. x -> Data[[i + 2, 1]]) == Data[[i + 2, 2]], {i, 1, k - 4}];
  Equations[[2 (k - 4) + 1]] = (Poly[[1]] /. x -> Data[[1, 1]]) == Data[[1, 2]];
  Equations[[2 (k - 4) + 2]] = (Poly[[1]] /. x -> Data[[2, 1]]) == Data[[2, 2]];
  Equations[[2 (k - 4) + 3]] = (Poly[[k - 3]] /. x -> Data[[k - 1, 1]]) == Data[[k - 1, 2]];
  Equations[[2 (k - 4) + 4]] = (Poly[[k - 3]] /. x -> Data[[k, 1]]) == Data[[k, 2]];
  Do[Equations[[2 (k - 3) + 2 + i]] = (D[Poly[[i]], x] /.  x -> intermediate[[i, 1]]) == (D[Poly[[i + 1]], x] /.  x -> intermediate[[i, 1]]), {i, 1, k - 4}];
  Do[Equations[[2 (k - 3) + 2 + k - 4 + i]] = (D[Poly[[i]], {x, 2}] /. x -> intermediate[[i, 1]]) == (D[Poly[[i + 1]], {x, 2}] /.x -> intermediate[[i, 1]]), {i, 1, k - 4}];
  un = Flatten[{atable, btable, ctable, dtable}];
  f = Solve[Equations];
  Poly = Poly /. f[[1]];
  xs = intermediate;
  xs = Prepend[Append[xs, Data[[k]]], Data[[1]]];
  pf = Table[{Poly[[i]], xs[[i, 1]] <= x <= xs[[i + 1, 1]]}, {i, 1, k - 3}];
  Piecewise[pf])
y = Spline3[Data]
yactual = 1/(1 + 25 x^2);
a = Plot[y, {x, -1, 1}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y3"}, ImageSize -> Medium, PlotRange -> All,  PlotLegends -> {"y3"}, PlotLabel -> "cubic spline"]
a = Plot[{y, yactual}, {x, -1, 1},  Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y3"}, ImageSize -> Medium, PlotRange -> All,  PlotLegends -> {"y3", "yactual"}, PlotLabel -> "cubic spline"]


Figure 10. Behaviour of the cubic spline interpolation scheme when applied to the Runge function data points

Figure 10. Behaviour of the cubic spline interpolation scheme when applied to the Runge function data points

Example

Use the cubic spline interpolation scheme presented above to find the interpolating function given the following 5 data points (-1, 0.038) (-0.8, 0.058), (-0.60, 0.10), (-0.4, 0.20), (-0.2, 0.5)

Solution

There are k=5 data points, therefore, there are k-3=5-3=2 cubic polynomials required for interpolation. These are:

    \[ y=\begin{cases}a_1+b_1x+c_1x^2+d_1x^3,&-1\leq x < -0.6\\ a_2+b_2x+c_2x^2+d_2x^3,&-0.6\leq x \leq -0.2 \end{cases} \]

To find the associated 4(2)=8 unknown coefficients, the following system of 8 linear equations is solved:

    \[ \begin{split} a_1-0.6b_1+0.36c_1-0.216d_1&=0.1\\ a_2-0.6b_2+0.36c_2-0.216d_2&=0.1\\ a_1-b_1+c_1-d_1&=0.038\\ a_1-0.8b_1+0.64c_1-0.512d_1&=0.058\\ a_2-0.4b_2+0.16c_2-0.064d_2&=0.2\\ a_2-0.2b_2+0.04c_2-0.008d_2&=0.5\\ b_1-1.2c_1+1.08d_1&=b_2-1.2c_2+1.08d_2\\ 2c_1-3.6d_1&=2c_2-3.6d_2 \end{split} \]

Solving the above system, yields the following interpolation function:

    \[ y=\begin{cases}0.453+0.967083x+0.75x^2+0.197917x^3,&-1\leq x < -0.6\\ 1.1685+4.54458x+6.7125x^2+3.51042x^3,&-0.6\leq x \leq -0.2 \end{cases} \]

Figure 11 shows the resulting interpolation function along with the five data points. The Mathematica code is provided below the figure.


Figure 11. Cubic spline interpolation applied to five data points

Figure 11. Cubic spline interpolation applied to five data points.

View Mathematica Code
Data = {{-1, 0.038}, {-0.8, 0.058}, {-0.60, 0.10}, {-0.4,  0.20}, {-0.2, 0.50}};
Clear[x, a, b, c, d]
Spline3[Data_] := (k = Length[Data];
  atable = Table[Subscript[a, i], {i, 1, k - 3}];
  btable = Table[Subscript[b, i], {i, 1, k - 3}];
  ctable = Table[Subscript[c, i], {i, 1, k - 3}];
  dtable = Table[Subscript[d, i], {i, 1, k - 3}];
  Poly = Table[atable[[i]] + btable[[i]]*x + ctable[[i]]*x^2 + dtable[[i]]*x^3, {i, 1, k - 3}];
  intermediate = Table[Data[[i + 2]], {i, 1, k - 4}];
  Equations = Table[0, {i, 1, 4 (k - 3)}];
  Do[Equations[[i]] = (Poly[[i]] /. x -> Data[[i + 2, 1]]) ==      Data[[i + 2, 2]], {i, 1, k - 4}];
  Do[Equations[[i + k - 4]] = (Poly[[i + 1]] /. x -> Data[[i + 2, 1]]) == Data[[i + 2, 2]], {i, 1, k - 4}];
  Equations[[2 (k - 4) + 1]] = (Poly[[1]] /. x -> Data[[1, 1]]) ==  Data[[1, 2]];
  Equations[[2 (k - 4) + 2]] = (Poly[[1]] /. x -> Data[[2, 1]]) ==  Data[[2, 2]];
  Equations[[2 (k - 4) + 3]] = (Poly[[k - 3]] /.  x -> Data[[k - 1, 1]]) == Data[[k - 1, 2]];
  Equations[[2 (k - 4) + 4]] = (Poly[[k - 3]] /. x -> Data[[k, 1]]) == Data[[k, 2]];
  Do[Equations[[2 (k - 3) + 2 + i]] = (D[Poly[[i]], x] /. x -> intermediate[[i, 1]]) == (D[Poly[[i + 1]], x] /. x -> intermediate[[i, 1]]), {i, 1, k - 4}];
  Do[Equations[[2 (k - 3) + 2 + k - 4 + i]] = (D[Poly[[i]], {x, 2}] /. x -> intermediate[[i, 1]]) == (D[Poly[[i + 1]], {x, 2}] /.  x -> intermediate[[i, 1]]), {i, 1, k - 4}];
  f = Solve[Equations];
  Poly = Poly /. f[[1]];
  xs = intermediate;
  xs = Prepend[Append[xs, Data[[k]]], Data[[1]]];
  pf = Table[{Poly[[i]], xs[[i, 1]] <= x <= xs[[i + 1, 1]]}, {i, 1, k - 3}];
  Piecewise[pf])
Equations // MatrixForm
y = Spline3[Data]
a = Plot[y, {x, -1, 0}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y3"}, ImageSize -> Medium, PlotRange -> All,  PlotLegends -> {"y3"}, PlotLabel -> "cubic spline"]

Leave a Reply

Your email address will not be published.