Lab 7 (13/09/2017)
INTERPOLATION
Three methods of interpolation have to be implemented in this lab.
- Piecewise Linear Interpolation
- Vandermonde Polynomial Interpolation
- Lagrange Polynomial Interpolation
Interpolation is the problem of computing values at unknown (or new) points given the data at several known points. Usually, it is done in two steps: first, compute an interpolating function or a set of interpolating functions that pass through the given data points; and second, evaluate the function(s) at the unknown points to compute the required values.
The different interpolation methods refer to the different ways of computing the interpolating function(s).
Preliminaries
We look at how the interpolation functions provided in Python (via the scipy module) work. This is to help in the analysis of the three interpolation functions that you implement in the lab.
First, we create a set of data points using the function $$ y = \sin(2x/5) + \sin(x^2/5) $$ and give the function evaluated at 11 points from -5.0 to +5.0 (including both) as the dataset.
Python code for generating the dataset is given below. The datasets are plotted as green circles to show how these 11 points are distributed in the 2D plane.
# Import the necessary modules
import numpy as np
import scipy.interpolate as sci
from matplotlib import pyplot as plt
from matplotlib import colors as mc
np.set_printoptions(precision=3, suppress=True)
# Create a dataset for testing interpolation functions
X = np.arange(-5, 5.1)
Y = np.sin(X*2/5.) + np.sin(X**2/5.)
print 'Data Points: ', X
print ' ', Y
# Let us create enough points for seeing the actual function
XA = np.arange(-5, 5.1, 0.05)
YA = np.sin(XA*2/5.) + np.sin(XA**2/5.)
plt.grid(True)
plt.plot(X, Y, 'go')
plt.show()
The functions interp1d() and lagrange() are used for interpolating 1-D data, i.e., functions of the form $y = f(x)$. They take two arrays: x- and y-coordinates as input arguments and return an object which stores all the parameters of the interpolating function(s). To compute a value at any $x$, we only need to call these objects with $x$ as the argument. $x$ can be an array of values, in which the case an array of interpolated values is returned (see all the way to the end of the section).
# Perform Linear, Lagrange and Cubic Spline interpolations
L = sci.interp1d(X, Y, kind='linear') # Linear interpolation
La = sci.lagrange(X, Y) # Lagrange interpolation
Cu = sci.interp1d(X, Y, kind='cubic') # Cubic interpolation
After computing the three different interpolating functions above, let us see how they perform. We create a set of points at which we will calculate the interpolated values. These are given in the array $IX$.
# Create a set of points for interpolation
# These points lie between 2 and 4 to see differences clearly
IX = np.array([2.0, 2.25, 2.5, 3.0, 3.333, 3.5, 3.667, 4.0])
Using the data points above, we use Linear Interpolation method to compute the values at those points. They are plotted along with the original data to see how Linear Interpolation performs.
As can be seen from the graph, linear interpolation joins all the points by straight lines (magenta coloured dots and lines) and does not really follow the true function given by the green dashed line.
# Plot for linear integration
plt.grid(True)
plt.plot(X, Y, 'go', XA, YA, 'g--', IX, L(IX), 'm.', IX, L(IX), 'm-')
plt.legend(['Data', 'True f', 'Linear data', 'Linear Fit'], loc='best')
plt.show()
LIV = L(IX) # Save interpolated values
We do the same for Lagrange interpolation. The blue dots and line represent the Lagrange polynomial used as the interpolating function. You can see that the dots are quite close to the true function indicating that the interpolated values are close to the true values.
# Plot for Lagrange interpolation
plt.grid(True)
plt.plot(X, Y, 'go', XA, YA, 'g--', IX, La(IX), 'b.', IX, La(IX), 'c-')
plt.legend(['Data', 'True f', 'Lagrange data', 'Lagrange Fit'], loc='best')
plt.show()
LaIV = La(IX) # Save interpolated values
Finally, we do the same analysis for Cubic Spline interpolation. The orange coloured dots and line represent the cubic spline interpolating function. Again the interpolated values are quite close to the true values (green dashed line).
# Plot for Cubic Spline interpolation
plt.grid(True)
plt.plot(X, Y, 'go', XA, YA, 'g--', IX, Cu(IX), 'r.')
plt.plot(IX, Cu(IX), color='coral')
plt.legend(['data', 'true f', 'CS interp', 'CS data'], loc='best')
plt.show()
CuIV = Cu(IX) # Save interpolated values
Now, we do error analysis to see how these three interpolation methods perform. The True and Interpolated values are all stored in arrays. The mean squared difference between True values and the interpolated values computed from different methods are calculated. The results show that Lagrange Interpolation gives the best result with an error of 5e-5, with Cubic Spline giving the next best performance.
However, remember that Lagrange and Vandermonde interpolations do not scale up well with larger numbers of data values.
TV = np.sin(IX*2/5.) + np.sin(IX**2/5) # Save True values
# Compare the true values with interpolated values
print 'True Values: ', TV
print 'Linear Interpolation: ', LIV
print 'Lagrange Interpolation: ', LaIV
print 'Cubic Spline Interpolation: ', CuIV
# Errors
print 'Mean Squared Error (Linear): ', ((LIV - TV)*(LIV - TV)).mean()
print 'Mean Squared Error (Lagrange): ', ((LaIV - TV)*(LaIV - TV)).mean()
print 'Mean Squared Error (Cubic Spline): ', ((CuIV - TV)*(CuIV - TV)).mean()
Problems
- Implement Piecewise Linear interpolation method and test on the datasets given in preliminaries. Compare your answers with the ones given by scipy methods.
In Linear implementation, you fit $n$ straight lines joining every successive pair of the $n+1$ data points given as input. The straight line parameters $a$ and $b$ ($y = ax + b$) are given by
$$ a_i = \frac{y_{i+1} - y_i}{x_{i+1} - x_i} $$ and $$ b_i = y_i $$ where the subscript $i$ indicates that it is the line joining the data points $(x_i, y_i)$ and $(x_{i+1}, y_{i+1})$. Thus, the equation of the interpolating line for $(x_i, y_i)$ is $$ y = \frac{y_{i+1} - y_i}{x_{i+1} - x_i} (x - x_i) + y_i $$ - Implement Vandermonde Polynomial interpolation. In Vandermonde implementation, we fit the $n$-th degree polynomial $$ P_n(x) = a_0 + a_1x + a_2x^2 + \cdots + a_nx^n $$ through the given dataset containing $n+1$ known values. The problem reduces to solving the linear system of $n+1$ equations in $n+1$ unknown $a$'s $$ \left(\begin{array}{ccccc} 1 & x_0 & x_0^2 & \cdots & x_0^n \\ 1 & x_1 & x_1^2 & \cdots & x_1^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{array}\right) \left(\begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_n \end{array}\right) = \left(\begin{array}{c} y_0 \\ y_1 \\ \vdots \\ y_n \end{array}\right) $$
- Implement Lagrange Polynomial interpolation. Lagrange polynomial is given by $$\begin{eqnarray*} P_n(x) & = & \frac{(x-x_1)(x-x_2)\cdots(x-x_n)}{(x_0-x_1)(x_0-x_2)\cdots(x_0-x_n)} y_0 + \\ & & ~~~~~\frac{(x-x_0)(x-x_2)\cdots(x-x_n)}{(x_1-x_0)(x_1-x_2)\cdots(x_1-x_n)} y_1 + \\ & & ~~~~~ \cdots + \frac{(x-x_0)(x-x_1)\cdots(x-x_{n-1})}{(x_n-x_0)(x_n- x_1)\cdots(x_n-x_{n-1})} y_n \end{eqnarray*}$$