====== Polynomial Interpolation via Newton Divided Differences ======

Newton's Divided Difference algorithm is a slick way to compute an $N$th order polynomial interpolant to a set of $N+1$ data points $(x_0, y_0), (x_1, y_1), \ldots, (x_N, y_N)$ with distinct $x_i$'s. 

It produces a polynomial in the form of [[gibson:teaching:fall-2016:math753:horner|Horner's method]] with base points, e.g. 

\begin{equation*}
y = P(x) = c_0 + (x - x_0) \, [c_1 + (x - x_1) [c_2 + (x - x_2) [c_3 + (x - x_3) \, c_4]]]
\end{equation*}

If we plug in the data points $(x_i, y_i)$ for $i=0,1,\ldots, N$, to the $N$th-order generalization of the above equation, we get a series of $N+1$ equations in the $N+1$ unknowns $c_0, \ldots c_N$.

\begin{eqnarray*}
y_0 &=& c_0 \\
y_1 &=& c_0 + (x_1-x_0) \, c_1 \\
y_2 &=& c_0 + (x_2-x_0) \, c_1 + (x_2-x_0)(x_2-x_1) \, c_2 \\
&\vdots&
\end{eqnarray*}

Lo and behold this is lower-triangular system of equations, which can be written in matrix form like this

\begin{equation*}
\left[\begin{array}{ccccc}
      1 &         & \ldots &        & 0  \\
      1 & x_1-x_0 &        &        &    \\
      1 & x_2-x_0 & (x_2-x_0)(x_2-x_1) &        & \vdots   \\
 \vdots & \vdots  &        & \ddots &    \\
      1 & x_k-x_0 & \ldots & \ldots & \prod_{j=0}^{N-1}(x_N - x_j)
\end{array}\right]
\left[\begin{array}{c}    c_0 \\  c_1 \\  c_2  \\     \vdots \\     \\     c_{N} \end{array}\right] =
\left[\begin{array}{c}    y_0 \\  y_1 \\  y_2 \\  \vdots \\ \\    y_{N} \end{array}\right]
\end{equation*}

Lower-triangular systems can be solved easily via forward substitution. It turns out that for this particular lower-triangular system, the solution can be computed easily by subtracting and dividing
numbers in a table. To see how that works, please refer to [[https://en.wikipedia.org/wiki/Newton_polynomial#Application|Newton Divided Difference Application]] (wikipedia).

Further reading:

  * [[https://en.wikipedia.org/wiki/Newton_polynomial#Application|Newton Polynomial]] (wikipedia).