Interpolation is the process of fitting a continuous function to a set of discrete data points for the purpose of estimating intermediate values. Polynomial interpolation involves fitting an \(n^{th}\)-order polynomial that passes through \(n+1\) data points (in order to use an \(n^{th}\)-order interpolating polynomial, \(n+1\) data points are required), and using \(b_i\) as a stand-in for the polynomial coefficients, \(b_{1}, b_{2}, \cdots ,b_{k}\) will be uniquely determined.
There are several methods that can be used to determine the form of the unique interpolating polynomial fit to a collection of discrete data points, but we’re only going to discuss Newton’s divided difference method. With this approach, each new polynomial in the algorithm is obtained from its predecessor by adding a new term. Thus, in its final state, the polynomial exhibits all the previous polynomials as constituents.
Before we begin, some clarification on the terminology used throughout the post: \(f_{n}(x)\) represents an \(n^{th}\)-order interpolating polynomial. For example, \(f_{1}(x)\), represents a first-order interpolating polynomial, \(f_{2}(x)\) a second-order interpolating polynomial, etc. Also, \(b_{k} = f[x_{0}, x_{1}, \cdots, x_{k}]\) represents the finite divided difference of order \(k\). Divided differences obey the formula:
Therefore, we can say \(f[x_{0}, x_{1}, \cdots, x_{k}]\) is the coefficient of \(x^{k}\) in the polynomial \(f_{k}(x)\). For a given set of data points \((x_{0}, f(x_{0})), (x_{1}, f(x_{1})), (x_{2}, f(x_{2}))\), the first three interpolating polynomial coefficients are given by:
Calculating Interpolation Coefficients via Divided Differences
Linear Approximation
The linear form of Newton’s Interpolating polynomial is given by
Note that the term \(\frac{f(x_{1})-f(x_{0})}{x_{1}-x_{0}}\) is a finite divided-difference approximation of the first derivative, which also represents the slope of the line of the first-degree polynomial approximation.
Let’s approximate \(\sqrt{x}\) between \(0 \leq x \leq 3\) using the points \((0,0),(3,1.732050808)\). The intercept \(b_{0}=f(x_{0})=0\), and the corresponding slope is \(b_{1}=\frac{1.732050808 - 0}{3-0} \approx.57735027\). The linear interpolant is therefore
Estimating \(\sqrt{2.25}\) with our linear approximation yields \(.57735027*2.25=1.2990381\), vs. the actual value of \(1.5\), resulting in an error of \(\epsilon_{1} = .134\).
Reducing the interval between data points (estimating the linear interpolant between \(2 \leq x \leq 2.5\)) will yield a better approximation with less error.
Quadratic Approximation
The second-order polynomial interpolant can be represented as
where, as above, \(b_0\), \(b_1\), and \(b_2\) are given by
\(b_1\) still represents the slope of the line connecting points \(x_0\) and \(x_1\). The new term \(b_{2}(x-x_{0})(x-x_{1})\) introduces second-order curvature, which tends to reduce the magnitude of the error in the approximation. Using the points \((0,0), (3,1.732050808), (5,2.236067978)\) to fit our quadratic interpolant results in coefficient estimates \(b_{0}=f(x_{0})=0\), \(b_{1} = .57735027\) and \(b_{2} = -.065068337\) yielding the second-order interpolating polynomial:
Evaluating \(f_{2}(2.25)=1.4088409\) yields a better approximation than the linear interpolant, reducing the error to \(\approx .06077\).
Newton’s Interpolating Polynomial: General Form
The \(n^{th}\)-order interpolating polynomial is given by
One approach used to determine coefficients for the \(n^{th}\) degree interpolating polynomial is to construct a table of finite differences, which for an \(n^{th}\) degree interpolating polynomial will have \(n+1\) distinct levels. Once the table has been populated, the coefficients can be obtained from the top-most diagonal, starting with \(f[x_{0}]\) for \(b_0\) and working toward the endpoint for higher indexed coefficients. Each forward node is determined by calculating \(f[x_{i},x_{j}] = \frac {f(x_{i}) - f(x_{j})}{x_{i}-x_{j}}\).
For example, here’s the grid resulting from estimating coefficients for the quadratic interpolant:
Then \(f[x_{0},x_{1}] = \frac {1.732050808 - 0}{3-0} = .5773502\), and \(f[x_{1},x_{2}] = \frac {2.236067978-1.732050808}{5.0-3.0} = .2520086\) and \(f[x_{0}, x_{1},x_{2}] = \frac {.2520086-.5773502}{5.0-0.0} = -.0650683\).
Extracting the top diagonal of coefficients from the tree above, we obtain
matching our estimate for the quadratic interpolant above.
What follows is a plot of \(\sqrt{x}\) along with the linear and quadratic interpolating polynomials:
Implementation
We can create a Python function which replicates the algorithm described above for an interpolating polynomial of arbitrary degree. In the following section, we present two such functions: The first, interp_coeffs, returns a tuple of interpolating coefficients generated from separate arrays of x and y values. Note that for \(n+1\) data points, interp_coeffs will return the coefficient estimates of an \(n\)-degree polynomial. The second, interp_val, takes the same x and y arrays, but instead of returning a tuple of coefficients returns a function which can then be passed any value within the range of the initial independent variable array, and returns the interpolated estimate evaluated at that point. This is exactly how scipy.interpolate.interp1d is designed, to which we’ll compare our implementation:
def interp_coeffs(xvals, yvals):
"""
Return the divided-difference interpolating polynomial
coefficients based on xvals and yvals. The degree of
the interpolating polynomial will be inferred from data,
which will be 1 less than the number of data points.
"""
nbr_data_points = len(xvals)
data = sorted(zip(xvals, yvals), reverse=False)
xvals, yvals = zip(*data)
depth = 1
coeffs = [yvals[0]]
iter_yvals = yvals
while depth < nbr_data_points:
iterdata = []
for i in range(len(iter_yvals) - 1):
delta_y = iter_yvals[i + 1] - iter_yvals[i]
delta_x = xvals[i + depth] - xvals[i]
iterval = (delta_y / delta_x)
iterdata.append(iterval)
if i==0:
coeffs.append(iterval)
iter_yvals = iterdata
depth+=1
return(coeffs)
def interp_val(xvals, yvals):
"""
Return a function representing the interpolating
polynomial determined by xvals and yvals.
"""
nbr_data_points = len(xvals)
data = sorted(zip(xvals, yvals), reverse=False)
xvals, yvals = zip(*data)
depth = 1
coeffs = [yvals[0]]
iter_yvals = yvals
while depth < nbr_data_points:
iterdata = []
for i in range(len(iter_yvals) - 1):
delta_y = iter_yvals[i+1] - iter_yvals[i]
delta_x = xvals[i + depth] - xvals[i]
iterval = (delta_y / delta_x)
iterdata.append(iterval)
if i==0:
coeffs.append(iterval)
iter_yvals = iterdata
depth+=1
def ff(i):
"""
Evaluate interpolated estimate at x.
"""
terms = []
retval = 0
for j in range(len(coeffs)):
iterval = coeffs[j]
iterxvals = xvals[:j]
for k in iterxvals: iterval*=(i-k)
terms.append(iterval)
retval+=iterval
return(retval)
return(ff)
Next we demonstrate how to obtain the interpolated estimates of \(\sqrt{x}\) using interp_coeffs and interp_val:
import numpy as np
# 1st-degree approximation.
x1 = [0, 3]
y1 = np.sqrt(x1)
# 2nd-degree approximation.
x2 = [0, 3, 5]
y2 = np.sqrt(x2)
# 3rd-degree approximation.
x3 = [0, 3, 5, 7]
y3 = np.sqrt(x3)
interp_coeffs(x1, y1) # [0.0, 0.5773502691896257]
interp_coeffs(x2, y2) # [0.0, 0.5773502691896257, -0.06506833684483389]
interp_coeffs(x3, y3) # [0.0, 0.5773502691896257, -0.06506833684483389, 0.007610943899867132]
The results agree with our earlier derivations.
Note that the output of interp_coeffs returns the intercept first, followed by higher-powered coefficients in order of increasing degree.
In order to obtain the interpolated estimate for some value within the valid range of x-values, we first pass the x and y arrays as we did with interp_coeffs, but we then pass the value of interest to the result of the first evaluation:
import numpy as np
# 1st-degree approximation.
x1 = [0, 3]
y1 = np.sqrt(x1)
# 2nd-degree approximation.
x2 = [0, 3, 5]
y2 = np.sqrt(x2)
# 3rd-degree approximation.
x3 = [0, 3, 5, 7]
y3 = np.sqrt(x3)
# Assume we're interested in obtaining 1st, 2nd and 3rd-
# degree interpolated estimates at x=2.75.
est1 = interp_val(x1, y1)
est1(2.75) # 1.5877132402714706
est2 = interp_val(x2, y2)
est2(2.75) # 1.632447721852294
>>> est3 = interp_val(x3, y3)
est3(2.75) # 1.644220900697401
We compare our results to scipy.interpolate.interp1d. The API for scipy.interpolate.interp1d is identical to interp_val, but takes an additional argument to specify what degree interpolant to estimate. We now verify that both functions return identical estimates:
import numpy as np
from scipy import interpolate
# 1st-degree approximation.
x1 = [0, 3]
y1 = np.sqrt(x1)
# 2nd-degree approximation.
x2 = [0, 3, 5]
y2 = np.sqrt(x2)
# 3rd-degree approximation.
x3 = [0, 3, 5, 7]
y3 = np.sqrt(x3)
# Using scipy's interpolate.interp1d to obtain interpolated estimates:
est1 = interp1d(x1, y1, kind='linear')
est1(2.75) # array(1.5877132402714706)
est2 = interp1d(x2, y2, kind='quadratic')
est2(2.75) # array(1.632447721852294)
est3 = interp1d(x3, y3, kind='cubic')
est3(2.75) # array(1.6442209006974007)