# Interpolation

### Big Idea.
An **interpolating function** provides a way to estimate values between and beyond a set of known data points. While there are infinitely many ways to create an interpolating function, we use constraints to find a unique, meaningful solution that can be computed efficiently. *Polynomial interpolation* is the simplest approach, while *cubic spline interpolation* offers greater flexibility.

### Interpolating Functions
:::info
**Definition** An *interpolating function* (or *interpolant*) for a set of points $(t_0,y_0),\dots,(t_d,y_d)$ is a function, $f(t)$, that passes through all of those points, meaning $f(t_k) = y_k$ for all $k$ from $0$ to $d$.
:::


## Polynomial Interpolation
:::info
**Definition**
* A **polynomial** of degree $d$ is a function of the form $$p(t)= c_0 + c_1 t + \cdots + c_d t^d,$$ where the coefficients $c_0, c_1, \dots, c_d$ are real numbers.
* The collection of all polynomials of degree (at most) $d$ is denoted by $$\mathbb{P}_d = \{ c_0 + c_1 t + \cdots + c_d t^d : c_0,c_1,\dots,c_d \in \mathbb{R} \}$$
* Given $d+1$ points $(t_0,y_0), \dots,(t_d,y_d)$, **polynomial interpolation** (with respect to the basis $\{\phi_i\}_{i=0}^d$) is a polynomial of the form $$p(t) = \sum_{i=0}^d c_i \phi_i(t)$$ such that $p(t_k) = y_k$ for each $k=0,\dots,d$.
* This process corresponds to a linear system of equations, $A \mathbf{c} = \mathbf{y}$, where
$$
A =
\begin{bmatrix}
\phi_0(t_0) & \phi_1(t_0) & \cdots & \phi_d(t_0) \\
\phi_0(t_1) & \phi_1(t_1) & \cdots & \phi_d(t_1) \\
\vdots & \vdots & \ddots & \vdots \\
\phi_0(t_d) & \phi_1(t_d) & \cdots & \phi_d(t_d)
\end{bmatrix}
\quad
\mathbf{c} = \begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_d \end{bmatrix}
\quad
\mathbf{y} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_d \end{bmatrix}
$$
:::

### Monomial Basis
:::info
**Definition**
* When a polynomial is expressed as $p(t) = \sum_{i=0}^d c_i t^i$, the functions $\{1,t,\dots,t^d\}$ are called the **monomial basis**.
* To find the coefficients, we can set up the linear system $A \mathbf{c} = \mathbf{y}$, where the matrix $$A = \begin{bmatrix} 1 & t_0 & \cdots & t_0^d \\ 1 & t_1 & \cdots & t_1^d \\ \vdots & \vdots & \ddots & \vdots \\ 1 & t_d & \cdots & t_d^d \end{bmatrix} $$ is known as the **Vandermonde matrix**.
:::

:::danger
**Theorem** Let $A$ be the Vandermonde matrix for points $(t_0,y_0), \dots,(t_d,y_d)$. Then
$$
\mathrm{det}(A) = \prod_{0 \leq i < j \leq d} (t_j - t_i)
$$
:::

:::danger
**Theorem**
Consider $d+1$ data points $(t_0,y_0), \dots , (t_d,y_d)$ such that $t_i \not= t_j$ for $i \not= j$. There exists a unique polynomial $p(t)$ of degree (at most) $d$ such that $p(t_k) = y_k$ for each $k=0,\dots,d$.
:::

:::warning
**Example** Find the unique polynomial of degree 3 which interpolates the points
$$
(0,-1),\ (1,-1),\ (2,1),\ (3,-1)
$$

* To find the unique third-degree polynomial, we solve the system $A \mathbf{c} = \mathbf{y}$:
$$
A =
\begin{bmatrix}
1 & 0 & 0 & 0 \\
1 & 1 & 1 & 1 \\
1 & 2 & 4 & 8 \\
1 & 3 & 9 & 27
\end{bmatrix}
\quad
\mathbf{y} = \left[ \begin{array}{r} -1 \\ -1 \\ 1 \\ -1 \end{array} \right]
$$
* The solution for the coefficients is
$$
\mathbf{c} = \left[ \begin{array}{r} -1 \\ -3 \\ 4 \\ -1 \end{array} \right]
$$
* Therefore, the interpolating polynomial is
$$p(t) = -1 - 3t + 4t^2 - t^3$$
:::

**Remark** While simple, the Vandermonde matrix becomes unstable with a large number of points due to a very large condition number, which can lead to computational instability.

### Lagrange Interpolation
:::info
**Definition** The **Lagrange interpolation polynomial** is defined as a sum of basis functions,
$$
l_h(t) = \frac{\prod_{ j \neq h}^d (t - t_j)}{\prod_{j \neq h}^d (t_h - t_j)},
\quad h = 0, \ldots, d
$$

where each basis function is a polynomial of degree $d$. These basis functions are constructed so that

* $l_h(t_i)=1$ when $i=h$ and $0$ when $i\neq h$.
* $p(t) = \sum_{h=0}^d y_j \cdot l_h(t)$
:::

:::warning
**Example**
For example, given the points $(1,1), (2,-4), (4,7)$, the Lagrange basis polynomials are
* $l_0(t) = \frac{(t-2)(t-4)}{(1-2)(1-4)} = \frac{t^2 - 6t + 8}{3}$
* $l_1(t) = \frac{(t-1)(t-4)}{(2-1)(2-4)} = \frac{t^2 - 5t + 4}{-2}$
* $l_2(t) = \frac{(t-1)(t-2)}{(4-1)(4-2)} = \frac{t^2 - 3t + 2}{6}$

The final polynomial is
$$
p(t) = l_0(t) - 4l_1(t) + 7l_2(t)
$$
:::

A significant drawback of this method is that adding a new data point requires recomputing all the basis functions.

### Newton Interpolation
The *Newton interpolation polynomial* addresses the problem of adding new points by using a different basis.

:::info
**Definition**
The **Newton interpolation polynomial** is defined as a sum of basis functions,
$$
\pi_0(t) = 1, \quad \pi_h(t) = (t - t_0)(t - t_1) \cdots (t - t_{h-1}), \; h=1, \cdots, d
$$
:::

The basis functions, $\pi_h(t)$, are constructed such that the matrix of the linear system, $A \mathbf{c} = \mathbf{y}$, is lower-triangular. This makes the system easy to solve using forward substitution.

A key advantage of the Newton method is its iterative nature. When you get a new data point, $(t_{d+1}, y_{d+1})$, you can simply add a new term to the existing polynomial without re-solving the entire system:

$$
p_{d+1}(t) = p_d(t) + c_{d+1} \pi_{d+1}(t)
$$

:::warning
**Example**
For the same points, $(1,1), (2,-4), (4,7)$, the Newton basis functions are:
* $\pi_0(t) = 1$
* $\pi_1(t) = t-1$
* $\pi_2(t) = (t-1)(t-2) = t^2 - 3t + 2$

The corresponding system of equations is:
$$
\begin{bmatrix}
1 & 0 & 0 \\
1 & 1 & 0 \\
1 & 3 & 6
\end{bmatrix}
\begin{bmatrix}
c_0 \\ c_1 \\ c_2
\end{bmatrix}=
\begin{bmatrix}
1 \\ -4 \\ 7
\end{bmatrix}
$$

Solving this system gives the coefficients
$$
c_0=1\ c_1=-5,\ c_2=\frac{7}{2}
$$

The interpolating polynomial is
$$
p(t) = \pi_0(t) - 5 \pi_1(t) + \tfrac{7}{2} \pi_2(t)
$$
:::

:::success
**Exercise**

1. Given $d$ distinct data points, which of the following statements about a unique interpolating polynomial is correct?
* (a) There exists a unique polynomial of degree (at most) $d−1$ which interpolates the data.
* (b) There exists a unique polynomial of degree (at most) $d$ which interpolates the data.
2. Given $d$ distinct data points, which of the following is true?
* (a) There exists a unique polynomial $p(t)$ of degree (at most) $d$ which interpolates the data and also satisfies $p'(t_1) = 0$ and $p''(t_1) = 0$.
* (b) There exists a unique polynomial $p(t)$ of degree (at most) $d$ which interpolates the data and also satisfies $p'(t_1) = 0$.
3. Consider $d+1$ distinct data points. Let $p(t)$, $q(t)$, and $r(t)$ be the interpolating polynomials constructed using the monomial, Lagrange, and Newton bases, respectively. Which of the following is true?
* (a) Then $p(t) = q(t) = r(t)$ for all $t$.
* (b) The polynomials $p(t), q(t), r(t)$ are distinct.
:::

## Cubic Spline Interpolation
A **cubic spline** is a piecewise function made up of cubic polynomials, where the function and its first and second derivatives are continuous. This method provides greater flexibility than polynomial interpolation, especially for larger data sets.

:::info
**Definition**
Consider $N+1$ points $(t_0,y_0),\dots,(t_N,y_N)$. A **cubic spline** is defined by $N$ cubic polynomials, $p_1(t),\dots,p_N(t)$, where:
$$
p_k(t) = a_k(t - t_{k-1})^3 + b_k(t - t_{k-1})^2 + c_k(t - t_{k-1}) + d_k, \quad t \in [t_{k-1},t_k]
$$

such that $p(t)$, $p'(t)$ and $p''(t)$ are continuous functions.
:::

Each polynomial is defined by four coefficients $a_k,b_k,c_k,d_k$, so we need to determine a total of $4N$ unknowns. These unknowns are found by setting up a linear system of equations based on a set of conditions.

* **Interpolation Conditions:** The spline must pass through all the data points.
* $p_k(t_{k-1}) = y_{k-1}$ (for the left endpoint)
* $p_k(t_k) = y_k$ (for the right endpoint)
* **Continuity Conditions:** The first and second derivatives of the polynomials must be continuous at the interior points where they connect.
* $p_k'(t_k) = p_{k+1}'(t_k)$ (continuity of the first derivative)
* $p_k''(t_k) = p_{k+1}''(t_k)$ (continuity of the second derivative)

These conditions provide $4N−2$ equations, so two more conditions are needed to ensure a unique solution.

:::info
**Definition** A **natural cubic spline** satisfies $p''_1(t_0) = p''_N(t_N) = 0$.
:::

:::danger
**Theorem**
Consider $N+1$ points $(t_0,y_0),\dots,(t_N,y_N)$ (with $t_i \not= t_j$ for $i \not= j$). The unique natural cubic spline $p(t)$ which interpolates the points is given by the coefficient matrix
$$
C=
\begin{bmatrix}
a_1 & a_2 & \cdots & a_N \\
b_1 & b_2 & \cdots & b_N \\
c_1 & c_2 & \cdots & c_N \\
d_1 & d_2 & \cdots & d_N
\end{bmatrix}
$$

where $d_k = y_{k-1}$ for $k=1,\dots,N$ and the coefficients $a_1,b_1,c_1,\dots,a_N,b_N,c_N$ are the solution of the linear system
$$
\left[
\begin{array}{c|c|c|c|c}
A(L_1)\vphantom{\bigg|} & B & & & \\ \hline
& A(L_2) & B & & \\ \hline
& & \ddots & \ddots & \phantom{A(L_{N-1})} \\ \hline
\phantom{A(L_{N-1})} & \phantom{A(L_{N-1})} & \phantom{A(L_{N-1})} & A(L_{N-1}) & B \\ \hline
T & & & & V
\end{array}
\right]
\begin{bmatrix} a_1 \\ b_1 \\ c_1 \\ \vdots \\ a_N \\ b_N \\ c_N \end{bmatrix}=
\begin{bmatrix} y_1 - y_0 \\ 0 \\ 0 \\ \vdots \\ y_N - y_{N-1} \\ 0 \\ 0 \end{bmatrix}
$$

where $L_k = t_k - t_{k-1}$ is the length of the subinterval $[t_{k-1},t_k]$ and

\begin{align}
A(L) &= \begin{bmatrix} L^3 & L^2 & L \\ 3L^2 & 2L & 1 \\ 6L & 2 & 0 \end{bmatrix} \quad
B = \left[ \begin{array}{rrr} 0 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & -2 & 0 \end{array} \right] \\
T &= \begin{bmatrix} 0 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad
V = \begin{bmatrix} L_N^3 & L_N^2 & L_N \\ 0 & 0 & 0 \\ 6L_N & 2 & 0 \end{bmatrix}
\end{align}
:::

:::warning
**Example** Consider the natural cubic spline $p(t)$ that interpolates the data points
$$
(0, 1) \ , \ (1, 3) \ , \ (2, 8) \ , \ (3, 10) \ , \ (4, 9) \ , \ (5,-1) \ , \ (6,-17)
$$


Suppose the coefficient matrix of $p(t)$ is
$$
C = \left[ \begin{array}{rrrrrr}
1 & -2 & 1 & \phantom{+}a_4 & 1 & 1 \\
0 & 3 & -3 & b_4 & -6 & -3 \\
1 & 4 & 4 & c_4 & -5 & -14 \\
1 & 3 & 8 & 10 & 9 & -1
\end{array} \right]
$$

Determine the coefficients $a_4, b_4, c_4$, then compute the value $p''(2.5)$.

* The $t$ values are equally spaced, so the length of each subinterval is $L_k=1$.
* The coefficients of the spline are determined by solving a linear system of equations
\begin{align}
p_k(t_k) = p_{k+1}(t_k) \ \ &\Rightarrow \ \ a_k + b_k + c_k + d_k = d_{k+1} \\
p'_k(t_k) = p'_{k+1}(t_k) \ \ &\Rightarrow \ \ 3a_k + 2b_k + c_k = c_{k+1} \\
p''_k(t_k) = p''_{k+1}(t_k) \ \ &\Rightarrow \ \ 6a_k + 2b_k = 2b_{k+1}
\end{align}
* Using the provided coefficient matrix, we can determine the unknown coefficients $a_4, b_4, c_4$.
* From $6a_3 + 2b_3 = 2b_4$, we get $b_4 = 0$.
* From $3a_3 + 2b_3 + c_3 = c_4$, we get $c_4 = 1$.
* From $6a_4 + 2b_4 = 2b_5$, we get $a_4 = -2$.
* To compute the value of $p''(2.5)$, we use the third polynomial, $p_3(t)$, since $2.5$ is in the interval $[t_2,t_3]$.
\begin{align}
p_3''(t) &= 6a_3 (t-t_2) + 2 b_3 \\
p''(2.5) &= p_3''(2.5) = 6(2.5 - 2) + 2(-3) = -3
\end{align}
:::

**Note:** The condition number of the matrix for constructing a natural cubic spline is significantly smaller than for a Vandermonde matrix. This makes cubic splines a much more stable and reliable method for interpolation, especially with a large number of data points. For example, with 11 equally spaced points, the Vandermonde matrix has a condition number of about $10^{12}$, while the cubic spline matrix has a condition number of only about $33$.


:::success
**Exercise**

4. For large values of $N$, do we expect $\mathrm{cond}(A_1) < \mathrm{cond}(A_2)$ or $\mathrm{cond}(A_1) > \mathrm{cond}(A_2)$, where $A_1$ is the matrix for the interpolating polynomial with respect to the monomial basis and $A_2$ is the matrix for the interpolating natural cubic spline?
5. Suppose we have 4 points $(0,y_0),(1,y_1),(2,y_2),(3,y_3)$ and want to interpolate the data using a spline constructed from 3 degree 2 polynomials. The spline requires that $p(t)$ and $p'(t)$ are continuous and $p''(t_0)=0$. Setup a linear system $A \mathbf{x} = \mathbf{b}$ where the solution is $$\mathbf{x} = \begin{bmatrix} a_1 & b_1 & a_2 & b_2 & a_3 & b_3 \end{bmatrix}^T$$
6. Suppose a cubic spline $p(t)$ interpolates the given data $$(0.0,1.0),(1.0,3.0),(2.0,1.0),(3.0,1.0),(4.0,2.0),(5.0,1.0)$$ and has the specified coefficient matrix
$$
\left[\begin{array}{rrrrr}
-1.19 & 1.93 & a_3 & -0.75 & 0.55 \\
0.00 & -3.56 & b_3 & 0.60 & -1.65 \\
3.19 & -0.37 & c_3 & 1.15 & 0.10\\
1.00 & 3.00 & 1.00 & 1.00 & 2.00
\end{array} \right]
$$
Determine the coefficients $a_3,b_3,c_3$.
:::