Interpolation of data sets







Numerical Methods

Interpolation of data sets

If a function is only known at a discrete set of data points $(x_{i}, y_{i})$, interpolation can be used to estimate the value of that function at values of $x$ not included in the data. An interpolation function $I(x)$ passes through the points of a discrete data set,

\begin{equation} I(x_{i})=y_{i}. \end{equation}

Usually we desire a function that smoothly connects the data point. One possibility is to use the polynomial of the least degree that passes through all of the points. Lagrange showed that this polynomial function is given by,

\[ \begin{equation} I(x) = \sum\limits_{k=1}^{N}y_k\prod\limits_{i=1\,(i\ne k)}^{N}\frac{x-x_i}{x_k-x_i}. \end{equation} \]

The app below calculates the polynomial fit through the series of points given in the text box. For data set 1, the interpolation function fits the data set well. However for data set 2, the interpolation function deviates far from the points even though it passes through every point. Try changing a data point in the data to see how the interpolation function changes.

Code that calculates this polynomial is shown below.

<script> function I(x,xi,yi) { //interpolation function N = xi.length; sum = 0; for (k=0; k<N; k++) { product = 1; for (i=0; i<N; i++) { if (i!=k) { product = product*(x-xi[i])/(xi[k]-xi[i]); } } sum += yi[k]*product; } return sum; } var xi = [-5,-4,-3,-2,-1,0,1,2,3,4,5]; var yi = [0.038,0.058,0.1,0.2,0.5,1,0.5,0.2,0.1,0.058,0.038]; document.write("I(0.2) = "+I(0.2,xi,yi)); </script>

lagrange_poly.py

<script> function I(x,xi,yi) { //interpolation function N = xi.length; sum = 0; for (k=0; k<N; k++) { product = 1; for (i=0; i<N; i++) { if (i!=k) { product = product*(x-xi[i])/(xi[k]-xi[i]); } } sum += yi[k]*product; } return sum; } var xi = [-5,-4,-3,-2,-1,0,1,2,3,4,5]; var yi = [0.038,0.058,0.1,0.2,0.5,1,0.5,0.2,0.1,0.058,0.038]; document.write("I(0.2) = "+I(0.2,xi,yi)); </script> #!/usr/bin/env python #source: http://pastebin.com/bNVcQt4x import numpy as np from matplotlib import pyplot as plt data_fname = 'knot_points.csv' # x1,y1 # x2,y2 # ... def read_data(fname): X = [] Y = [] with open(fname, 'r') as f: for line in f: (x, y) = line.split(',') X.append(float(x)) Y.append(float(y)) return (X, Y) def langrange_polynomial(X, Y): def L(i): return lambda x: np.prod([(x-X[j])/(X[i]-X[j]) for j in range(len(X)) if i != j]) * Y[i] Sx = [L(i) for i in range(len(X))] # summands return lambda x: np.sum([s(x) for s in Sx]) # F = langrange_polynomial(*read_data(data_fname)) (X, Y) = read_data(data_fname) F = langrange_polynomial(X, Y) x_range = np.linspace(X[0], X[-1], 100) plt.plot(X, Y, 'ro') plt.plot(x_range, map(F, x_range)) plt.xlabel(r'$x$') plt.ylabel(r'$F(x)$') plt.title('Lagrange polynomial interpolation') plt.grid(True) plt.show()

For the polynomial fit, it is the higher order polynomials that sometimes causes the interpolation function to deviate quite far from the data. To avoid this problem, the order of the polynomials is limited. Two commonly used approaches are linear interpolation and cubic splines. These methods of interpolation are discussed in the next two sections.