# 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
:   `<math>`{=html}

\\begin{align}
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{align} `</math>`{=html} 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}]+fx_{0}, x_{1}+fx_{0}, x_{1}, x_{2}(x-x_{1})+...$$

$$+fx_{0}, x_{1}, x_{2}, ...,x_{n}(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:

!A table for solving the coefficients of a Newton\'s
polynomial.{width="800"}

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{align}y_0&=f[x_0]\\
                                                   &=6\end{align}   & \begin{align}f[x_{0}, x_{1}] &=\frac{f[x_{1}]-f[x_{0}]}{x_{1}-x_{0}}\\ 
                                                                                               &=\frac{11-6}{2-1}\\
                                                                                               &=5 \end{align}& \begin{align} 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{align}& \begin{align}&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{align}\\
\hline
1     & x_{1}=2 & \begin{align}y_1&=f[x_1]\\
                                                   &=11\end{align} & \begin{align}f[x_{1}, x_{2}] &=\frac{f[x_{2}]-f[x_1]}{x_{2}-x_{1}}\\
                                                                                               &=\frac{18-11}{3-2}\\
                                                                                               &=7 \end{align} & \begin{align} 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{align}& \\
\hline
2     & x_{2}=3 & \begin{align}y_2&=f[x_2]\\
                                                   &=18\end{align} & \begin{align}f[x_{2}, x_{3}] &=\frac{f[x_{3}]-f[x_2]}{x_{3}-x_{2}}\\
                                                                                               &=\frac{27-18}{4-3}\\
                                                                                               &=9 \end{align} &    & \\
\hline
3     & x_{3}=4 & \begin{align}y_3&=f[x_3]\\
                                                   &=27\end{align} &    &    & \\
\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{align}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{align}$$

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\le m\le 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{align}
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{align}$$ Obviously the function curve passes both data points. The
first derivative also matches that of the divided difference method:

$$\begin{align}
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{align}$$

# 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.

!Interpolation example
linear{width="400"}

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

!Quadratic spline six
segments{width="400"}

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}]$.

