Introduction to Numerical Methods/Interpolation

Objectives

• understand interpolation
• derive Newton’s divided difference method of interpolation
• derive Lagrangian method of interpolation
• apply the interpolation methods to solve problems
• find derivatives and integrals of discrete functions using interpolation

Resources

• Direct method chapter
• Newton's method chapter
• Lagrange method chapter
• Spline method chapter
• A Chronology of Interpolation
• NIST Digital Library of Mathematical Functions

Interpolation

Interpolation is the process of deriving a simple function from a set of discrete data points so that the function passes through all the given data points (i.e. reproduces the data points exactly) and can be used to estimate data points in-between the given ones. It is necessary because in science and engineering we often need to deal with discrete experimental data. Interpolation is also used to simplify complicated functions by sampling data points and interpolating them using a simpler function. Polynomials are commonly used for interpolation because they are easier to evaluate, differentiate, and integrate - known as polynomial interpolation.

It can be proven that given n+1 data points it is always possible to find a polynomial of order/degree n to pass through/reproduce the n+1 points.

Direct Method

Given n+1 data points the direct method assumes the following polynomial:

$y=f(x)=a_{0}+a_{1}x+a_{2}x^{2}+...+a_{n}x^{n}$

With n+1 values for $x$ and the n+1 corresponding values for $y$ we can solve for $a_{0},a_{1},...,a_{n}$ by solving the n+1 simultaneous linear equations, which is known as the direct method.

For example, given two data points $(x_{0},y_{0})$ and $(x_{1},y_{1})$ we can use a linear function $y=f(x)=a_{0}+a_{1}x$ to pass through the two data points:

$y_{0}=a_{0}+a_{1}x_{0}$

$y_{1}=a_{0}+a_{1}x_{1}$

Once we solve for $a_{0}$and $a_{1}$ (the coefficients of $f(x)$) we can use the function as the basis for interpolation - estimating the missing data points in-between.

Newton's Method

In Newton's method the interpolating function is written in Newton polynomial(a.k.a Newton form). For example, given one data point $(x_{0},y_{0})$ we can only derive a polynomial of order zero: $y=f(x)=a_{0}$. Because $f(x_{0})=y_{0}$ the newton polynomial of order zero is $y=f(x)=y_{0}$.

Given two data points we can write Newton's polynomial in the form of $y=f(x)=a_{0}+a_{1}(x-x_{0})$. Plugging in the two data points gives us

$a_{0}=y_{0}$

$a_{1}={\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}$

Obviously this function passes through both data points (i.e. accurately reproduces the two data points). The first derivative of the function at $x=x_{0}$ is $f'(x_{0})=a_{1}={\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}$, which matches the result from the forward divided difference method.

Given three data points we can write Newton's polynomial in the form of

$y=f(x)=a_{0}+a_{1}(x-x_{0})+a_{2}(x-x_{1})(x-x_{0})$

Plugging in the three data points gives us

As $y_{0}=a_{0}$ we get $a_{0}=y_{0}$

As $y_{1}=a_{0}+a_{1}(x_{1}-x_{0})$ we get $a_{1}={\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}$

As $y_{2}=a_{0}+a_{1}(x_{2}-x_{0})+a_{2}(x_{2}-x_{1})(x_{2}-x_{0})$

we get

${\begin{aligned}a_{2}&={\frac {y_{2}-a_{0}-a_{1}(x_{2}-x_{0})}{(x_{2}-x_{1})(x_{2}-x_{0})}}\\&={\frac {y_{2}-y_{0}-{\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}(x_{2}-x_{0})}{(x_{2}-x_{1})(x_{2}-x_{0})}}\\&={\frac {y_{2}-y_{0}-{\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}((x_{2}-x_{1})+(x_{1}-x_{0}))}{(x_{2}-x_{1})(x_{2}-x_{0})}}\\&={\frac {y_{2}-y_{0}-{\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}(x_{2}-x_{1})-(y_{1}-y_{0})}{(x_{2}-x_{1})(x_{2}-x_{0})}}\\&={\frac {y_{2}-y_{1}-{\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}(x_{2}-x_{1})}{(x_{2}-x_{1})(x_{2}-x_{0})}}\\&={\frac {{\frac {y_{2}-y_{1}}{x_{2}-x_{1}}}-{\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}}{x_{2}-x_{0}}}\\\end{aligned}}$

This third degree polynomial function passes all three data points (the second derivative and the third derivative at $x_{0}$ and $x_{1}$match that from the divided difference method).

From the two examples we can see the coefficients of a Newton polynomial follow a pattern known as divided difference. For example$a_{1}={\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}$ is called divided difference of order 1 (denoted as $f[x_{0},x_{1}]$) because it depends on $x_{0}$ and $x_{1}$. The divided difference notation can be used to write the general order (form) Newton polynomial:

$f(x)=f[x_{0}]+f[x_{0},x_{1}](x-x_{0})+f[x_{0},x_{1},x_{2}](x-x_{0})(x-x_{1})+...$

$+f[x_{0},x_{1},x_{2},...,x_{n}](x-x_{0})(x-x_{1})(x-x_{2})...(x-x_{n-1})$

Where

$a_{0}=f[x_{0}]=y_{0}$ because $f[x_{i}]=y_{i}$ by definition

$a_{1}=f[x_{0},x_{1}]={\frac {f[x_{1}]-f[x_{0}]}{x_{1}-x_{0}}}$

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

$\dots$

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

We can calculate the coefficients of Newton polynomial using the following table

${\begin{matrix}x_{0}&y_{0}=f[x_{0}]&&&\\&&f[x_{0},x_{1}]&&\\x_{1}&y_{1}=f[x_{1}]&&f[x_{0},x_{1},x_{2}]&\\&&f[x_{1},x_{2}]&&f[x_{0},x_{1},x_{2},x_{3}]\\x_{2}&y_{2}=f[x_{2}]&&f[x_{1},x_{2},x_{3}]&\\&&f[x_{2},x_{3}]&&\\x_{3}&y_{3}=f[x_{3}]&&&\\\end{matrix}}$

because

$f[x_{i}]=y_{i}$

$f[x_{i},x_{i+1}]={\frac {f[x_{i+1}]-f[x_{i}]}{x_{i+1}-x_{i}}}$

The following figures shows the dependency between the divided differences:

For example, given data points $(1,6)$, $(2,11)$, $(3,18)$, and $(4,27)$ we can draw the following table:

${\begin{array}{|c||c|c|c|c|c|}i&x_{i}&y_{i}&f[x_{i},x_{i+1}]&f[x_{i},x_{i+1},x_{i+2}]&f[x_{i},x_{i+1},x_{i+2},x_{i+3}]\\\hline 0&x_{0}=1&{\begin{aligned}y_{0}&=f[x_{0}]\\&=6\end{aligned}}&{\begin{aligned}f[x_{0},x_{1}]&={\frac {f[x_{1}]-f[x_{0}]}{x_{1}-x_{0}}}\\&={\frac {11-6}{2-1}}\\&=5\end{aligned}}&{\begin{aligned}f[x_{0},x_{1},x_{2}]&={\frac {f[x_{1},x_{2}]-f[x_{0},x_{1}]}{x_{2}-x_{0}}}\\&={\frac {7-5}{3-1}}\\&=1\end{aligned}}&{\begin{aligned}&f[x_{0},x_{1},x_{2},x_{3}]\\&={\frac {f[x_{1},x_{2},x_{3}]-f[x_{0},x_{1},x_{2}]}{x_{3}-x_{0}}}\\&={\frac {1-1}{4-1}}\\&=0\end{aligned}}\\\hline 1&x_{1}=2&{\begin{aligned}y_{1}&=f[x_{1}]\\&=11\end{aligned}}&{\begin{aligned}f[x_{1},x_{2}]&={\frac {f[x_{2}]-f[x_{1}]}{x_{2}-x_{1}}}\\&={\frac {18-11}{3-2}}\\&=7\end{aligned}}&{\begin{aligned}f[x_{1},x_{2},x_{3}]&={\frac {f[x_{2},x_{3}]-f[x_{1},x_{2}]}{x_{3}-x_{1}}}\\&={\frac {9-7}{4-2}}\\&=1\end{aligned}}&\\\hline 2&x_{2}=3&{\begin{aligned}y_{2}&=f[x_{2}]\\&=18\end{aligned}}&{\begin{aligned}f[x_{2},x_{3}]&={\frac {f[x_{3}]-f[x_{2}]}{x_{3}-x_{2}}}\\&={\frac {27-18}{4-3}}\\&=9\end{aligned}}&&\\\hline 3&x_{3}=4&{\begin{aligned}y_{3}&=f[x_{3}]\\&=27\end{aligned}}&&&\\\end{array}}$

The four data points lie on a polynomial of order 2, which is why the coefficient $a_{3}$($f[x_{0},x_{1},x_{2},x_{3}]$) is zero. Given $[a_{0},a_{1},a_{2}]=[6,5,1]$ the result polynomial is:

${\begin{aligned}y&=a_{0}+a_{1}(x-x_{0})+a_{2}(x-x_{1})(x-x_{0})\\&=6+5\times (x-1)+1\times (x-2)(x-1)\\&=x^{2}+2x+3\end{aligned}}$

Once the coefficients of a Newton's polynomial are solved, we can evaluate the polynomial function for any $x$. A Newton polynomial is often rewritten in a nested form:

$f(x)=a_{0}+(x-x_{0})(a_{1}+(x-x_{1})(a_{2}+(x-x_{2})(a_{3}+...+a_{n}(x-x_{n-1}))...)$,

because this nested form of interpolating polynomial is easier to evaluate because x only appears in the function n times.

For example, the nested form of a third order interpolating polynomial is:

$f(x)=a_{0}+(x-x_{0})(a_{1}+(x-x_{1})(a_{2}+(x-x_{2})a_{3}))$

The algorithm of Newton's method and its implementation can be found in this Jupyter notebook.

Lagrange Form

Lagrange polynomial is another form used for polynomial interpolation. It is called a form because with a given set of distinct points the interpolating polynomial is unique. We can arrive at the same polynomial through different methods.

The Lagrange form specifies the interpolation polynomial as:

$f_{n}(x):=\sum _{i=0}^{n}L_{i}(x)f(x_{i})$

$L_{i}(x):=\prod _{\begin{smallmatrix}0\leq m\leq n\\m\neq i\end{smallmatrix}}{\frac {x-x_{m}}{x_{i}-x_{m}}}={\frac {(x-x_{0})}{(x_{i}-x_{0})}}\cdots {\frac {(x-x_{i-1})}{(x_{i}-x_{i-1})}}{\frac {(x-x_{i+1})}{(x_{i}-x_{i+1})}}\cdots {\frac {(x-x_{n})}{(x_{i}-x_{n})}},$

where $n$ is the order of the polynomial.

For example, given two data points we get $n=1$ and

${\begin{aligned}f(x)&=L_{0}f(x_{0})+L_{1}f(x_{1})\\&={\frac {x-x_{1}}{x_{0}-x_{1}}}f(x_{0})+{\frac {x-x_{0}}{x_{1}-x_{0}}}f(x_{1})\\\end{aligned}}$

Obviously the function curve passes both data points. The first derivative also matches that of the divided difference method:

${\begin{aligned}f'(x)&={\frac {f(x_{0})}{x_{0}-x_{1}}}+{\frac {f(x_{1})}{x_{1}-x_{0}}}\\&={\frac {f(x_{1})-f(x_{0})}{x_{1}-x_{0}}}\\\end{aligned}}$

Spline Interpolation

Spline interpolation uses a number of polynomial functions to interpolate a set of data points with each polynomial for two adjacent data points. The Spline method is necessary because often times when the order of the polynomial become large polynomial interpolation shows oscillatory behavior (instability known as Runge's phenomenon). The following iPython notebook shows an example that suffers this issue: [1]

Spline method is not another method for finding polynomial interpolation of a discrete function, but instead it results in a piecewise polynomial (splines) in order to avoid the oscillatory behavior. The most common spline interpolations are linear, quadratic, and cubic splines. Linear interpolation uses lines to connect each pair of consecutive data points resulting in a piecewise interpolation.

A quadratic spline uses a quadratic polynomial to connect consecutive data points.

A function f(x) is a quadratic spline if the following conditions are true:

1. The domain of $f(x)$ is an interval [a, b].
2. $f(x)$ and $f'(x)$ are continuous on [a, b].
3. The data points $x_{i}$ such that $a=x_{0}<x_{1}<\dots <x_{n}=b$ and $f(x)$ is a polynomial of order at most 2 on each subinterval $[x_{i},x_{i+1}]$.