\documentclass[a4paper]{amsart}
\usepackage{graphicx}
\usepackage{amsmath, amsthm}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{tikz}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup

\newtheorem{thm}{Theorem}[section]
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{defn}{Definition}[section]
\newtheorem{example}{Example}[section]
\newtheorem{exercise}{Exercise}[section]
\newtheorem{algorithm}{Algorithm}[section]


\newcommand\be{\begin{equation}}
\newcommand\ee{\end{equation}}
\newcommand\half{\frac{1}{2}}

\newcommand\CC{{\mathbb C}}
\newcommand\EE{{\mathbb E}}
\newcommand\NN{{\mathbb N}}
\newcommand\OO{{\mathcal O}}
\newcommand\PP{{\mathbb P}}
\newcommand\PPbar{{\overline{\mathbb P}}}
\newcommand\RR{{\mathbb R}}
\newcommand\TT{{\mathbb T}}
\newcommand\ZZ{{\mathbb Z}}

%\renewcommand{\theenumi}{\roman{enumi}}    

\DeclareMathOperator{\conv}{co}
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator{\dist}{dist}
\DeclareMathOperator{\GL}{GL}
\DeclareMathOperator{\Span}{Span}
\DeclareMathOperator{\var}{var}
\DeclareMathOperator{\Vol}{Vol}

\newcommand\fhat{{\widehat{f}}}
\newcommand\ghat{{\widehat{g}}}
\newcommand\Lhat{{\widehat{L}}}
\newcommand\LAhat{{\widehat{L^A}}}
\newcommand\Lhhat{{\widehat{L^h}}}
\newcommand\phihat{{\widehat{\phi}}}
\newcommand\phiAhat{{\widehat{\phi^A}}}
\newcommand\phihhat{{\widehat{\phi^h}}}


\newcommand{\va}{{\bf a}}
\newcommand{\vb}{{\bf b}}
\newcommand{\vc}{{\bf c}}
\newcommand{\vd}{{\bf d}}
\newcommand{\vf}{{\bf f}}
\newcommand{\vh}{{\bf h}}
\newcommand{\ve}{{\bf e}}
\newcommand{\vg}{{\bf g}}
\newcommand{\vp}{{\bf p}}
\newcommand{\vq}{{\bf q}}
\newcommand{\vs}{{\bf s}}
\newcommand{\vu}{{\bf u}}
\newcommand{\vU}{{\bf U}}
\newcommand{\vUhat}{\widehat{\bf U}}
\newcommand{\vv}{{\bf v}}
\newcommand{\vw}{{\bf w}}
\newcommand{\vx}{{\bf x}}
\newcommand{\vy}{{\bf y}}
\newcommand{\vz}{{\bf z}}
\newcommand{\vzhat}{\widehat{\bf z}}
\newcommand{\vO}{{\bf 0}}

\newcommand{\vL}{{\bf L}}
\newcommand{\vT}{{\bf T}}

\newcommand{\Balpha}{{\mbox{\boldmath $\alpha$}}}
\newcommand{\Bbeta}{{\mbox{\boldmath $\beta$}}}
\newcommand{\Bgamma}{{\mbox{\boldmath $\gamma$}}}
\newcommand{\Bell}{{\mbox{\boldmath $\ell$}}}
\newcommand{\Blambda}{{\mbox{\boldmath $\lambda$}}}
\newcommand{\Bmu}{{\mbox{\boldmath $\mu$}}}

%
% Uncomment EXACTLY ONE of the next two commands - 
% Uncomment 1st line to avoid text, 2nd line to get text
%
\newcommand{\ignore}[1]{}
%\newcommand{\ignore}{}
%
%
%





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 

\begin{document}

\title[Problems in Mathematics: Newton and Polynomial Interpolation]
{Problems in Mathematics: Newton and Polynomial Interpolation}
\author[Brad Baxter]{Brad Baxter}


\maketitle



\begin{abstract}
These notes provide a short introduction to divided differences and the Newton
form of the interpolating polynomial for Problems in Mathematics 2021. 
\end{abstract}

\tableofcontents

\vfill
\dedicatory{\bf Version: 202111081142}

\newpage
\noindent
\section{Introduction}
You can
access these notes, and other material,
via my office machine:

\begin{verbatim}
http://econ109.econ.bbk.ac.uk/brad/Problems_in_Maths/
\end{verbatim}

Some students might also be interested in the software packages Matlab and Octave, mentioned
in the lecture.

\begin{verbatim}
http://econ109.econ.bbk.ac.uk/brad/Methods/matlab_intro_notes.pdf
\end{verbatim}

You can find out more about Numerical Analysis in my lecture notes for
a course given at Imperial from 1994--2001:

\begin{verbatim}
http://econ109.econ.bbk.ac.uk/brad/teaching/Methods/nabook.pdf
\end{verbatim}

\newpage
\section{Historical Introduction}
We begin roughly $400$ years ago, at the beginning of the Seventeenth
Century. At that time, applied mathematics was expanding rapidly, from
providing better tables of trigonometric functions for European
oceanic navies, to improved ways to calculate interest on debts. This
might sound mundane, but war and money have always been closely linked
to mathematical applications. Algebra was a Sixteenth century
invention, but almost nothing was known of calculus in
1600. Fortunately, this was about to change!

To gain some idea of the needs of the time we need to understand the
recent history of computation. My phone can compute $10^8$ FLOPS
(floating point operations per second) and cost $\pounds 100$,
but in 1990 a research
workstation would be limited to $10^6$ FLOPS and cost $\pounds
10^4$. Going back to the 1970s, schoolchildren were still taught how
to use trigonometric and logarithmic tables, which would have looked
something like this ($x$ is measured in degrees):

\[
\begin{array}{cc}
  x & \sin x\\
  0 & 0.0000\\
  1 & 0.0175\\
  2 & 0.0349\\
  3 & 0.0523\\
  4 & 0.0698\\
  5 & 0.0872\\
\end{array}
\]

Such four-figure tables
required enormous work to produce in 1700, but what do we
do if we need to calculate $\sin 1.4$? One simple way is linear
interpolation using the table: we let $f(x) = \sin x$ and define
\[
p(x) = f(1) + \left(f(2)-f(1)\right) \left(x-1\right),
\]
that is,
\be
p(x) = 0.0175 + 0.0174 \left(x-1\right).
\label{pp1}\ee
Using \eqref{pp1}, we obtain
\[
p(1.4) = 0.02446.
\]
The true value is $\sin 1.4 = 0.02443 \ldots$, so this isn't too
bad. Can we do better?

One idea is to use three values and fit a quadratic. Let $a = 1$,
$b=2$ and $c=4$ and let's write the quadratic in the form
\be
q(x) = p(x) + Q (x-a)(x-b),
\label{pp2}\ee
where $p(x)$ is given by \eqref{pp1}. The idea here is to avoid
further computation, since we already know that $p(a) = \sin a$ and
$p(b) = \sin b$. Since the quadratic term $(x-a)(x-b)$ vanishes at $a$
and $b$, the quadratic $q(x)$ already satisfies $q(a) = \sin a$ and
$q(b) = \sin b$. To reproduce the value at $c$, we solve
\[
\sin c = p(c) + Q (c-a)(c-b),
\]
or
\[
Q = \frac{\sin c - \sin a - (\sin b - \sin a)(c-a)}{(c-a)(c-b)}.
\]
Working to five significant figures, we find
\[
Q = 1.6667 \times 10^{-5}.
\]





\newpage
\section{Polynomial Interpolation}

Let 
\[
z_0, z_1, \ldots, z_n
\]
be any different complex numbers and let 
\[
f_0,
\ldots, f_n
\]
be any complex numbers (they don't need to be distinct). 
We want to construct a
polynomial $p$ of degree $n$ for which 
\[
p(z_j) = f_j, 
\quad\hbox{ for } 0 \le j \le n.
\]
Such a polynomial is called an {\bf interpolating polynomial}, and
we say that $p$ {\bf interpolates} the data 
\[
(z_0,f_0), (z_1, f_1), \ldots, (z_n, f_n).
\]
We shall let $\PP_n$ denote the vector space of
polynomials of degree $n$.

\begin{example}\label{pex1}
How do we find the quadratic polynomial satisfying $p(0) = \alpha$, $p(1) = \beta$
and $p(4) = \gamma$? We could just substitute $p(x) = p_0 + p_1 x +
p_2 x^2$ and solve the three linear equations to obtain the
coefficients. However, there is a simpler solution: we write
\[
 p(z) = \alpha {(z-1)(z-4) \over (0-1)(0-4)} +
        \beta  {z(z-4) \over (1-0)(1-4)} +
        \gamma {z(z-1) \over (4-0)(4-1)}.
\]
\end{example}

The key point in Example \ref{pex1} is that the polynomial
\[
\ell_0(z) = {(z-1)(z-4) \over (0-1)(0-4)}
= \frac{1}{4}(z-1)(z-4)
\]
satisfies
\[
\ell_0(1) = \ell_0(4) = 0 
\quad\hbox{ and }\quad
\ell_0(0) = 1.
\]
There is a useful notation for this called the {\bf Kronecker delta},
written $\delta_{jk}$. Specifically, if we let $z_0 = 0$, $z_1 = 1$
and $z_2 = 4$, then we have
\[
\ell_0(z_k) = \delta_{0k},
\]
where
\[
\delta_{jk} =
\begin{cases} 1 & \text{if $j = k$,}\\
0 &\text{if $j \ne k$.}
\end{cases}
\label{delta1}
\]
Thus
\[
\ell_1(z) = {z(z-4) \over (1-0)(1-4)}
= -\frac{1}{3}z(z-4) 
\quad\hbox{ satisfies }\quad
\ell_1(z_k) = \delta_{1k}
\]
and
\[
\ell_2(z) = {z(z-1) \over (4-0)(4-1)}
= \frac{1}{12}z(z-1)
\quad\hbox{ satisfies }\quad
\ell_2(z_k) = \delta_{2k}.
\]

\begin{example}\label{pex1a}
Suppose now that $\alpha = \beta = \gamma = 1$ in Example \ref{pex1}.
We have shown that
\begin{align*}
1
&= \ell_0(z) + \ell_1(z) + \ell_2(z)\\
&= \frac{1}{4}(z-1)(z-4) 
+ \frac{-1}{3}z(z-4)
+ \frac{1}{12}z(z-1).
\end{align*}
If we now divide both sides by $z(z-1)(z-4)$, then we obtain
\[
\frac{1}{z(z-1)(z-4)}
= \frac{\left(1/4\right)}{z} +
\frac{\left(-1/3\right)}{z-1} +
\frac{\left(1/12\right)}{z-4}.
\]
\end{example}

With Kronecker delta notation to hand, we can now generalise the trick
used in Example \ref{pex1}.

\begin{lem}\label{pdef1}
Let
\be
\label{p1}
\ell_j(z) = \prod_{k=0, k \ne j}^n {z - z_k \over z_j - z_k},
\qquad 0 \le j \le n.
\ee
Then $\ell_r(z_s) = \delta_{rs}$, $0 \le r, s, \le n$ and $\ell_r \in \PP_n$.
\end{lem}

\begin{proof}
By construction, $\ell_r(z_s) = 0$ when $r \ne s$, because the product
in (\ref{p1}) contains the term $(z-z_s)$. However, $\ell_r(z_r)=1$,
because then every term in (\ref{p1}) occurs in both numerator and denominator.
\end{proof}

These polynomials $\ell_0, \ell_1, \ldots, \ell_n$ are useful because
they allow us to write down a very simple expression for the
polynomial interpolant.

\begin{prop}\label{pprop1}
The interpolating polynomial $p \in \PP_n$ for the data $\{(z_j,f_j):
0 \le j \le n\}$ is given by
\be
 p(z) = \sum_{j=0}^n f_j \ell_j(z), \qquad z \in \CC.
\label{p2}
\ee
\end{prop}

\begin{proof}
Equation \ref{p1} implies $p(z_k) = \sum_{j=0}^n f_j \delta_{jk} =
f_k$, $0 \le k \le n$.
\end{proof}

Equation \ref{p2} is called the {\bf Lagrange form of the
  interpolating polynomial}. Unfortunately, the Lagrange form is
often of limited use in practical work, although it's extremely useful in
theoretical work. 




Uniqueness requires a simple lemma.

\begin{lem}
Let $p(z) = a_0 + a_1 z + a_2 z^2 + \cdots + a_n z^n$, where $z \in
\CC$ and $a_0, a_1, \ldots, a_n \in \CC$. Then $p(z)$ has at most $n$
distinct zeros in $\CC$ unless $a_0 = a_1 = \cdots = a_n = 0$.
\end{lem}

\begin{proof}
The lemma is plainly true when $n=0$ or $n=1$. We then proceed by
induction. Thus let us assume that every polynomial of degree less
than $n$ has at most $n$ different zeros, unless every coefficient is
zero. Given any polynomial $p(z)$ of degree $n+1$, {\bf either} $p$
has a root, say $p(w) = 0$, {\bf or} $p$ has no roots. If the latter
condition is valid, then there's nothing further to demonstrate. If
the former is valid, then $(z-w)$ is a factor of $p(z)$. Thus we can
write $p(z) = q(z) (z-w)$, and the roots of $p$ are $w$ and the roots
of $q$. However, by induction hypothesis, $q$ can have at most $n$
different roots. Thus $p$ can have, in total, at most $n+1$ different roots.
\end{proof}

The last lemma is a very simple version of the
great Fundamental Theorem of Algebra: a polynomial of degree $n$ with
complex coefficients has exactly $n$ complex zeros if we count
multiple zeros multiply. (Thus $(z-2)^2$ has two zeros.) 

\begin{prop}\label{pprop2}
There is exactly one interpolating polynomial $p \in \PP_n$ when the
points $z_0, z_1, \ldots, z_n$ are distinct.
\end{prop}

\begin{proof} Existence was shown in Proposition \ref{pprop1}, so we address
uniqueness. Thus let $p$ and $q$ be interpolating polynomials of
degree $n$. Their difference $p-q$ is a polynomial of degree $n$ that
vanishes at the $n+1$ different points $z_0, \ldots, z_n$. Hence $p-q$
vanishes identically, using the last lemma.
\end{proof}


\begin{example}\label{pex2}
You've already met the Lagrange form when computing partial
fractions. Let $w_0, w_1, \ldots, w_n$ be different complex
numbers. We shall compute the scalars $\alpha_0, \alpha_1, \ldots,
\alpha_n$ in the partial fraction decomposition
\[
 {1 \over (z-w_0)(z-w_1)\cdots(z-w_n)}
 = \sum_{j=0}^n {\alpha_j \over z - w_j}.
\]
Let us set $f(z) \equiv 1$. The Lagrange form of the polynomial
interpolating $f$ at $w_0, w_1, \ldots, w_n$ is
\[
 1 = \sum_{j=0}^n \ell_j(z), 
\]
by \ref{pprop1}. Dividing both sides by $(z-w_0)\cdots(z-w_n)$ yields
the expression
\[
 \alpha_j = \Bigl( \prod_{k=0, k\ne j}^n (w_j - w_k) \Bigr)^{-1},
\qquad j = 0, 1, \ldots, n.
\]
\end{example}

The Lagrange form of the interpolating polynomial is useful when $n$
is small and in theoretical work. However, it is particularly
inconvenient if we have constructed $p_{n-1} \in \PP_{n-1}$
interpolating data $\{ (z_j,f_j) : 0 \le j \le n-1 \}$ and are then
given a new datum $(z_n,f_n)$, because we almost have to start the
calculation from scratch. Fortunately a more compact form is
available. The key idea is to let $p \in \PP_n$ take the form
\be
 p_n(z) = p_{n-1}(z) + C(z-z_0)(z-z_1)\cdots(z-z_{n-1}),
\qquad z \in \CC.
\label{p3}
\ee
We see that $p_n(z_j) = p_{n-1}(z_j) = f_j$, for $0 \le j \le n-1$, so
we do not disturb our previous interpolant at these points. Of course
we choose $C$ to satisfy the equation
\be
 f_n = p_{n-1}(z_n) + C \prod_{k=0}^{n-1}(z_n - z_k).
\label{p4}
\ee
Obviously $C$ depends on $f$ and $z_0, z_1, \ldots, z_n$. A
traditional notation is 
\be
 C = f[z_0, z_1, \ldots, z_n],
\label{p5}
\ee
so that \ref{p3} becomes
\be
 p_n(z) = p_{n-1}(z) 
   + f[z_0,z_1, \ldots, z_n] (z-z_0)(z-z_1)\cdots(z-z_{n-1}).
\label{p6}
\ee
The number $f[z_0, \ldots, z_n]$ is called a {\bf divided difference},
because of the method used to calculate these numbers described
below. Note that the coefficient of highest degree for $p_n$ does not
depend on the order in which we take the points. In other words, if we
replace $z_0, z_1, \ldots, z_n$ by $z_{\pi 0}, z_{\pi 1}, \ldots,
z_{\pi n}$, for any permutation $\pi$ of the numbers $\{0, 1, \ldots,
n\}$, then $f[z_{\pi 0}, \ldots, z_{\pi n}] = f[z_0, \ldots,
z_n]$. Another way to see this is the following explicit expression
for $f[z_0, \ldots, z_n]$, which is sometimes useful in theoretical
work.

\begin{prop}\label{pprop3}
We have
\be
  f[z_0, z_1, \ldots, z_n]
 = \sum_{j=0}^n {f(z_j) \over \prod_{k=0,k\ne j}^n (z_j-z_k)}.
\label{p6a}
\ee
Further, $f[z_0, \ldots, z_n] = 0$ when $f$ is a polynomial of degree
less than $n$.
\end{prop}

\begin{proof}
We just equate the coefficients of $z^n$ in $p_n(z) = \sum_{j=0}^n
f(z_j) \ell_j(z)$, using Proposition \ref{pprop1}. Moreover, if $f(z)
= z^\ell$ and $\ell < n$, then the coefficient of $z^n$ in $p_n$ is
zero. But this highest degree coefficient is $f[z_0, \ldots, z_n]$.
\end{proof}


Recurring equation \ref{p6}, and defining $f[z_0] = f(z_0)$, yields
the explicit expression
\begin{eqnarray*}
 p_n(z) = 
   f[z_0] +
   f[z_0, z_1] (z-z_0) +
   f[z_0, z_1, z_2] (z-z_0)(z-z_1) + \cdots && \\
   + f[z_0, z_1, \ldots, z_n](z-z_0)(z-z_1)\cdots(z-z_{n-1}), &&
\end{eqnarray*}
and this is called the {\bf Newton form} of the interpolating polynomial.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

It's {\bf important} to understand that $f[z_0, \ldots, z_\ell]$ is the
coefficient of highest degree for the polynomial $p_\ell \in \PP_\ell$
interpolating the data $\{ (z_k, f_k) : 0 \le k \le \ell\}$. 

\begin{example}\label{pex3}
The Newton form of the quadratic polynomial satisfying $p(0) =
f(0)$, $p(1) = f(1)$ and $p(4) = f(4)$  is 
\[
 p(z) = f[0] + f[0,1] z + f[0,1,4]z(z-1). 
\]
You'll see how to calculate the coefficients shortly.
\end{example}


The recursion used to calculate divided difference and justifying the
suitability of their name is derived in the following key theorem.

\begin{thm}\label{pthm4}
For any distinct complex numbers $z_0,
z_1, \ldots, z_n, z_{n+1}$ the divided differences satisfy
\be
  f[z_0, \ldots, z_{n+1}] 
 = { f[z_0, \ldots, z_n] - f[z_1, \ldots, z_{n+1}] \over z_0 -
z_{n+1}}. 
\label{p7}\ee
\end{thm}

\begin{proof}
We introduce two polynomials: (i) $p \in \PP_n$ interpolates 
$\{(z_k,f_k):0 \le k \le n\}$, and (ii) $q \in \PP_n$
interpolates $\{(z_k, f_k) : 1 \le k \le n+1\}$. Thus the
coefficients of highest degree for $p$ and $q$ are 
$f[z_0, \ldots, z_n]$ and $f[z_1, \ldots, z_{n+1}]$, respectively. The
key {\bf trick} is now the observation that the polynomial $r \in \PP_{n+1}$
interpolating at all $n+1$ points satisfies
\be
 r(z) = 
 { (z - z_{n+1}) p(z) - (z - z_0) q(z) \over z_0 - z_{n+1}}, 
\label{p8}\ee
because it is unique, by Proposition \ref{pprop2}, and it is easily
checked that the right hand side of (\ref{p8})
interpolates at $z_0, \ldots, z_{n+1}$: an exercise for the reader. 
Now the coefficient of highest degree in $r$ is $f[z_0,
\ldots, z_{n+1}]$, so equating the coefficients of highest degree in
(\ref{p8}) yields (\ref{p7}).
\end{proof}


Proposition \ref{p3} is extremely important. It's the basis of the 
algorithm for divided differences and the Newton form of the interpolating
polynomial. One way to illustrate the divided difference table is as follows.

\begin{center}
\begin{tikzpicture}
  \node (s00) at (0,0){$f(x_2)$};
  \node (s04) at (0,4){$f(x_1)$}; 
  \node (s08) at (0,8){$f(x_0)$}; 
  
  \node (s22) at (2,2){$f[x_1,x_2]$};
  \node (s26) at (2,6){$f[x_0,x_1]$};

  \node (s44) at (4,4){$f[x_1,x_2,x_3]$};


  \draw[red,thick,->] (s00) -- (s22);
  
  \draw[red,thick,->](s04)--(s22);

  \draw[red,thick,->] (s04) -- (s26);
  
  \draw[red,thick,->](s08)--(s26);

  \draw[red,thick,->] (s22) -- (s44);
  
  \draw[red,thick,->](s26)--(s44);

 
\end{tikzpicture}
\end{center}


\newpage
\begin{example}
Let $f(x) = x^3$ and let $x_0=0$, $x_1 = 1$ and $x_2 = 4$. The divided
difference table is as follows.

\begin{center}
\begin{tikzpicture}
  \node (s00) at (0,0){$64$};
  \node (s04) at (0,4){$1$}; 
  \node (s08) at (0,8){$0$}; 
  
  \node (s22) at (2,2){$\frac{64-1}{4-1} = 21$};
  \node (s26) at (2,6){$\frac{1-0}{1-0}=1$};

  \node (s44) at (4,4){$\frac{21-1}{4-0} = 5$};


  \draw[red,thick,->] (s00) -- (s22);
  
  \draw[red,thick,->](s04)--(s22);

  \draw[red,thick,->] (s04) -- (s26);
  
  \draw[red,thick,->](s08)--(s26);

  \draw[red,thick,->] (s22) -- (s44);
  
  \draw[red,thick,->](s26)--(s44);

 
\end{tikzpicture}
\end{center}

Thus $f[x_0, x_1] = 1$, $f[x_0, x_1, x_2] = 5$ and the Newton form of
the quadratic interpolating $f(x) = x^3$ at $0$, $1$ and $4$ is given
by
\[
p(x) = x + 5x(x-1).
\]
\end{example}

In linear algebra terms we form the lower triangular matrix
\[
\left(
\begin{array}{ccccccc}
 x_0 & f(x_0) & & & & &  \\
 x_1 & f(x_1) & f[x_0,x_1] & & & &  \\
 x_2 & f(x_2) & f[x_1,x_2] & f[x_0,x_1,x_2] & & &  \\
 \vdots & \vdots & \vdots & \vdots & \ddots & &  \\
 x_{n-1} & f(x_{n-1}) & f[x_{n-2},x_{n-1}] &
f[x_{n-3},x_{n-2},x_{n-1}] & \ldots &f[x_0,x_1,\ldots,x_{n-1}] &\\
x_n & f(x_n) & f[x_{n-1},x_n] &
f[x_{n-2},x_{n-1},x_n] & \ldots &f[x_1,x_2,\ldots,x_n] & f[x_0, \ldots, x_n]\\
\end{array}
\right).
\]

\noindent
In practice, only the first two columns of this lower triangular
matrix are stored. The diagonal elements are those needed for the
Newton form of the interpolating polynomial, so it's usual for 
the second column to be overwritten by subsequent
columns; that's $\OO(n)$ rather than $\OO(n^2)$ numbers to be stored. At
completion, the second column should contain the diagonal elements of the
matrix, that is 
\[
[f(x_0), f[x_0,x_1], f[x_0,x_1,x_2], f[x_0,x_1,x_2,x_3],
\ldots, f[x_0,x_1,\ldots,x_{n-1}], f[x_0, \ldots, x_n]\  ]^T. 
\]




A polynomial can be evaluated at any point using {\bf Horner's
rule} (nested multiplication), which we state here for completeness.



\medskip\noindent
{\bf Algorithm for polynomial evaluation:} Let $p(z) = a_0 + a_1 z + a_2 z^2 +
\cdots a_n z^n$ and let $\alpha$ be any complex number. The following
algorithm calculates the value $p(\alpha)$ of the polynomial at $z = \alpha$.

\hskip 1in Set $s := a_n$.

\hskip 1in For $k = n-1, n-2, \ldots, 1, 0$, set

\hskip 1.5in $s := s*z + a_k$.
\medskip

\noindent This is faster than naively evaluating $p(\alpha) = a_0 +
a_1*\alpha + a_2*\alpha*\alpha + \cdots + a_n * \alpha * \cdots * \alpha$.

\begin{exercise}
  [{\bf [Unassessed]}]
  Assume you're using a computer for which every arithmetic operation
takes the same amount of time. Calculate (i) the number of operations
required to apply Horner's rule, and (ii) the number of operations for
the naive method.
\end{exercise}

What about the error in polynomial interpolation?

\begin{thm}\label{pthm5}
Let $p \in \PP_n$ interpolate $f$ at $n$
distinct complex numbers $z_0, z_1, \ldots, z_n$. Then the error $e =
f - p$ satisfies the equation
\be
 e(w) = f[z_0, z_1, \ldots, z_n, w] \prod_{k=0}^n (w - z_k), \qquad
w \in \CC. 
\label{p9}\ee
\end{thm}

\begin{proof}
If we add a new interpolation point $z_{n+1}$, then the Newton
interpolating polynomial $q \in \PP_{n+1}$ is given by
\[
 q(z) = p(z) + f[z_0, z_1, \ldots, z_n, z_{n+1}] \prod_{k=0}^n (z -
z_k). 
\]
Hence
\[ f(z_{n+1}) = p(z_{n+1}) + 
f[z_0, z_1, \ldots, z_n, z_{n+1}] \prod_{k=0}^n (z -
z_k). 
\]
Since $z_{n+1}$ can be {\it any} point, we can write $w = z_{n+1}$,
which completes the proof.
\end{proof}

This result is of little use for error bounds unless we can bound
$f[z_0, z_1, \ldots, z_n, w]$ from above in some way.
Now the first mean
value theorem implies the equation
\[
 f[x_0,x_1] = {f(x_1) - f(x_0) \over x_1 - x_0} = f^\prime(\alpha), 
\]
for some point $\alpha \in [x_0,x_1]$. There is an important 
result for divided differences generalising this remark that's essentially a
form of the mean value theorem you'll meet in real analysis. We shall use this
relation to express the error in terms of the maximum modulus of the $(n+1)$st
derivative of $f$.

\begin{thm}\label{pthm6}
Let $f$ have continuous $(n+1)$st
derivative and let $x_0 <  x_1 <  \cdots < x_n$ be 
{\bf real} numbers.
Then there is a point $\alpha \in [x_0,x_n]$ such that
\be
 f[x_0, x_1, \ldots, x_n] = {f^{n}(\alpha) \over n!}. 
\label{p10}\ee
\end{thm}

\begin{proof}
Let $p_n \in \PP_n$ interpolate $f$ at $x_0, \ldots, x_n$. 
Then the error function $e = f - p_n$ has at least $n+1$ zeros in
$[x_0,x_n]$. Hence
its derivative $e^\prime$ has at least $n$ 
zeros in $[x_0,x_n]$, and its second
derivative $e^{\prime\prime}$ has at least $n-1$ zeros. Continuing in
this way, we deduce that $e^{(n)}$ has at last one zero, $\alpha$ say,
in $[x_0,x_n]$. But then
\[
 0 = e^{(n)}(\alpha) = f^{(n)}(\alpha) - f[x_0, \ldots, x_n] n!, 
\]
as required.
\end{proof}

\begin{cor}\label{pcor7}
Let $f$ have continuous $(n+1)$st
derivative and let $x_0, x_1, \ldots, x_n$ be different real numbers.
If $p_n \in \PP_n$ is the interpolating polynomial, then the error 
$e_n = f - p_n$ satisfies
\be
 |e_n(x)| \le {M \prod_{k=0}^n |x - x_k|
\over (n+1)!}, \qquad x \in [a,b], 
\label{p11}\ee
where $M = \max \{ |f^{(n+1)}(t)| : a \le t \le b \}$.
\end{cor}

\begin{proof}
This is immediate from the last two theorems.
\end{proof}

%\begin{example}
%Rework the last example for general $a$ and $b$.
%\end{example}

%\begin{example}
%Let $f(x) = \cos x$. Rework the last example.
%\end{example}

These examples might suggest that increasing the number of
interpolation points always decreases the error. This is {\bf not} so,
as you may see in exercises.

\begin{example}
Let $f(x) = \exp(x)$ and let $a=-1/2$, $b=1/2$. If the
interpolation points are always contained within the interval $[-1/2,1/2]$, then
the error of interpolation satisfies
$$ |e_n(x)| \le {e \over (n+1)!}, \qquad -1/2 \le x \le 1/2. $$
In other words, the error is {\it tiny}, what ever the choice of interpolation
points. In fact, this is true whenever the function being interpolated is
complex differentiable at every point of the complex plane. It is certainly
{\bf not} true for general functions, as we shall shortly see.
\end{example}

Equation (\ref{p11}) suggests the following problem: Find interpolation
points $(x_k)_{k=0}^{n-1}$ minimising
\be
 \max_{-1 \le x \le 1} \left( \prod_{k=0}^{n-1} |x - x_k| \right), 
\label{p12}\ee
which occurs when we want to minimise upper bound (\ref{p11}) on the
interval $[-1,1]$. Equally spaced points are particularly bad.
In fact, the minimum value of (\ref{p12}) occurs when 
\[ 
 \prod_{k=0}^{n-1} \left(x - x_k\right) = 2^{1-n} \cos( n \cos^{-1}x).
\]
This was discovered by the great Russian mathematician Chebyshev and
is called a Chebyshev polynomial.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\begin{thebibliography}{99}

\bibitem{dB}
  de~Boor, C. (2001), ``A Practical Guide to Splines'', Springer.
  
\bibitem{CdB}
Conte, S.~D., and C.~de~Boor (1980), ``Elementary Numerical Analysis'', 
McGraw Hill.

\bibitem{Davis}
  Davis, P.~J. (1963), ``Interpolation and Approximation'', Dover.

\bibitem{HairerWanner}
  Hairer, E. and G.~Wanner (1996), ``Analysis by its History'',
  Springer.

\bibitem{trefethen1}
  Trefethen, L.~N. (2013), ``Approximation Theory and Approximation
  Practice'', SIAM.

\end{thebibliography}

\vfill

\address{Department of Economics, Mathematics and Statistics,
Birkbeck College, University of London,
Malet Street, London WC1E 7HX, 
England}

\email{b.baxter@bbk.ac.uk}




\end{document}