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:

$$f[x_{0}, x_{1}, \cdots, x_{k}] = \frac{f[x_{1},x_{2},\cdots,x_{k}]-f[x_{0},x_{1},\cdots,x_{k-1}]}{x_{k}-x_{0}}$$

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:

\begin{align} b_{0} = f(x_{0}) && b_{1} = \frac{f(x_{1})-f(x_{0})}{x_{1}-x_{0}} && b_{2} = \frac{\frac{f(x_{2})-f(x_{1})}{x_{2}-x_{1}}-\frac{f(x_{1})-f(x_{0})}{x_{1}-x_{0}}}{x_{2}-x_{0}} \end{align}

We’ll walkthrough the linear and quadratic cases prior to exploring the general case.

## Calculating Interpolation Coefficients Using Divided Differences

### Linear Approximation

The linear form of Newton’s Interpolating polynomial is given by:

$$f_{1}(x)= f_{0}(x) + \frac{f(x_{1})-f(x_{0})}{x_{1}-x_{0}}(x-x_{0})$$

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 the $$\sqrt{x}$$ between the points $$x=0$$ and $$x=3$$ using the points $$(0, 0), (3, 1.732050808)$$.

The intercept is $$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:

$$f_{1}(x) = b_{0}+b_{1}x = .57735027x$$

Estimating $$\sqrt{2.25}$$ with our linear approximation yields $$.57735027*2.25=1.2990381$$, vrs. the actual value $$\sqrt{2.25}=1.5$$, which results in an error of $$\epsilon_{1} = .134$$.

Note that reducing the interval between data points (estimating the linear interpolant between $$x=2.00$$ and $$x=2.50$$ for example) will yield a better approximation with less error.

Next we’ll estimate the quadratic interpolant starting with same two data points taken from $$\sqrt{x}$$.

The second-order polynomial interpolant can be represented as:

$$f_{2}(x) = b_{0}+b_{1}(x - x_{0}) + b_{2}(x-x_{0})(x-x_{1}),$$

where, as above, $$b_{0}$$, $$b_{1}$$ and $$b_{2}$$ are given by:

\begin{align} b_{0} = f(x_{0}) && b_{1} = \frac{f(x_{1})-f(x_{0})}{x_{1}-x_{0}} && b_{2} = \frac{\frac{f(x_{2})-f(x_{1})}{x_{2}-x_{1}}-\frac{f(x_{1})-f(x_{0})}{x_{1}-x_{0}}}{x_{2}-x_{0}} \end{align}

$$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 should 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:

$$f_{2}(x) = b_{0} + b_{1}(x - x_{0}) + b_{2}(x-x_{0})(x-x_{1}) = 0 +.57735027(x-0) -.065068337(x-0)(x-3),$$

Evaluating $$f_{2}(2.25)=1.4088409$$ yields a better approximation of $$\sqrt{2.25}$$ than the linear interpolant, and reduces the error to $$\approx .06077$$.

## Newton’s Interpolating Polynomial: General Form

The $$n^{th}$$-order interpolating polynomial is given by:

$$f_{n}(x) = b_{0} + b_{1}(x-x_{0}) + b_{2}(x-x_{0})(x-x_{1}) \cdots + b_{n}(x-x_{0})(x-x_{1})\cdots(x-x_{n-1})$$

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. Starting with the set of x & y points intended for interpolation, the table looks like:

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 that results 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

$$f_{2}(x) = 0.0 + .5773502(x-0) -.0650683(x-0)(x-3) = \mathbf{.5773502x - .0650683x(x-3)},$$

matching our estimate for the quadratic interpolant above.

Here’s 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 outputting 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 leveraged, to which we’ll compare our estimate shortly:

#===========================================================
# Netwon's Divided Difference Interpolation 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.
"""
assert len(xvals)==len(yvals)
nbr_data_points = len(xvals)

# sort inputs by 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)

# append top-most entries from table to coeffs =>
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.
"""
assert len(xvals)==len(yvals)
nbr_data_points = len(xvals)

# sort inputs by 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)

# append top-most entries in tree to coeffs =>
if i==0: coeffs.append(iterval)

iter_yvals = iterdata
depth+=1

def f(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(f)


This implementation is nothing more than a recreation of Neville’s Algorithm, which is described in greater detail here.

Returning to our attempt at approximating $$\sqrt{x}$$ with various polynomial interpolants, next we demonstrate how to obtain the interpolated estimates via interp_coeffs and interp_val:

import numpy as np

# 1st-degree approximation
x1 = [0, 3]
y1 = np.sqrt(x1).tolist()

# 2nd-degree approximation
x2 = [0, 3, 5]
y2 = np.sqrt(x2).tolist()

# 3rd-degree approximation
x3 = [0, 3, 5, 7]
y3 = np.sqrt(x3).tolist()

>>> 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 exactly 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. Using
# interp_val, we'd proceed as follows:

>>> 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 primary method used for 1-dimensional interpolation in scipy. 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 the same set of results:

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(2.75)
array(1.632447721852294)

>>> est3 = interp1d(x3, y3, kind='cubic')
>>> est3(2.75)
array(1.6442209006974007)


One difference between the two implementations is that interpolate.interp1d returns the estimate in the form of an array. This is due to the fact that interp1d can take an array of inputs in addition to single-valued input, and so returns all estimates in the form of an array. I/O differences aside, note that both functions return identical results for estimates of all degrees tested.

A word of caution: Runge’s phenomenon is a problem of oscillation at the edges of an interval that occurs when using polynomial interpolation with polynomials of high degree 1. If 3rd or 4th-degree interpolating polynomials do not provide a ‘natural’ enough fit to the data, consider instead cubic spline interpolation, which will be covered in a future post. Until then, happy coding!