\documentclass[twoside]{article}

\usepackage{palatino}
\usepackage{camnum}

\newcommand{\ip[1]}{\langle{#1}\rangle}
\def \R {\BB{R}}
\def \X {\BB{X}}
\def \PP {\BB{P}}
\def \l {\ell}


\textheight=25.5cm
\textwidth=15cm
\topmargin=-2cm
\oddsidemargin=0.5cm
\evensidemargin=0.5cm
\parindent=0pt
\parskip=6pt

\setcounter{section}{7}
\setcounter{subsection}{4}
\setcounter{equation}{1}

\begin{document}
\begin{center}
  {\large \textbf{Mathematical Tripos Part IB: Easter 2001}}

  {\Large \textbf{Numerical Analysis -- Lecture
      11\footnote{Corrections and suggestions to these notes should
        be emailed to \texttt{a.shadrin@damtp.cam.ac.uk}. All handouts
        are available on the WWW at the URL 
        \texttt{http://www.damtp.cam.ac.uk/user/na/PartIB/Handouts.html}.}}}
\end{center}

\subsection{Least-squares fitting to discrete function values}

Suppose we are given $m+1$ function values
$f(x_0),\ldots,f(x_m)$ (e.g., the results of observations), 
where the $x_k$s are pairwise distinct, and we expect that
$f$ should coin with a polynomial of degree $n$. If $m \le n$,
then one can use interpolation. If $m > n$, then the least
squares problem can be solved.

{\bf Least squares problem.} Given $(f(x_i))_{i=0}^m$, find  
$
     p^* = {\rm arg} \min_{p\in P_n} 
           \sum_{i=0}^m |f(x_i) - p(x_i)|^2 w(x_i).
$

For the weight function $w$ one usually chooses $w\equiv 1$,
the different values $w(x_i)$ are being used when accuracy
of data depends on $x_i$ (e.g., to supress inaccurate data).
With 
$\|f\|:= \ip[f,f]^{1/2}$ and
$\ip[f,g]_w := \sum_{i=0}^m f(x_i)g(x_i)w(x_i)$,
it is simply the minimization problem 
$p^* = {\rm arg} \min_{p\in P_n} \|f-p\|$.

{\bf Solution.}
One possibility is to express $p$ as $\sum_{j=0}^n c_j x^j$
and find optimal $c_j$s as a solution of a linear least squares
problem (see \S 5), using QR factorization. 
An alternative is to construct orthogonal polynomials 
w.r.t.\ the scalar product $\ip[f,g]_w$. 
The theory is identical to that of \S\S 7.1--4, 
except that we have enough data to evaluate $p_k$ only for $k\le m$.
However, we need just $p_k$ for $k \le n$ and $n < m$. Thus \\
\hspace*{2ex}
1. Employ the recurrence $(7.3)$ to calculate monic $(\tilde p_k)_{k=0}^n$ \\
\hspace*{2ex}
2. Normalize $p_k = \tilde p_k/\ip[\tilde p_k,\tilde p_k]_w^{1/2}$ 
   and  form $p(x)= \sum_{k=0}^n \ip[p_k,f] p_k(x)$.\\
Since the work for each $k$ is bounded by a constant multiple of $m$,
the complete  cost is $\O{mn}$, as compared with $\O{n^2m}$ if QR is
used.

\section{Approximation of linear functionals}

Given a vector space $\X$ (e.g., $\R^n$, or $C^n[a,b]$), 
the {\em linear functional} is a linear mapping 
$$
\mbox{
    $\l: \X \to \R$\qquad (s.t.
    $\l(\alpha f + \beta g) = \alpha \l(f) + \beta \l(g)$, \quad
    $\forall f,g\in \X$).
}
$$
We will treat the space $\X = C^m[a,b]$, and the functionals
$$
  \mbox{
  1) \quad $\l(f) = f^{(k)}(x)$, \quad 
%     $x\in[a,b]$, \quad $0\le k\le n$ \quad (both fixed), \qquad 
      ($x,k$ being fixed), \qquad 
  2) \quad $\l(f)= \int _a^b f(x)w(x) dx$}. 
$$
Suppose we are given $n$ values $(f(x_i))_{i=0}^n$ only. 
Our goal is to find an approximation
$$
\textstyle
   \l(f) \approx \sum_{i=0}^n a_i f(x_i).
\eqno{(8.1)}
$$
For the functionals 1)--2), this is called {\em numerical 
differentiation} and {\em numerical integration}, respectively.
A natural solution is to interpolate $f$ by $p\in\PP_{n}$ and
take $\l(f) \approx \l(p)$. This is called 
the {\em interpolating formula} (of degree $n$).
We have already seen an interpolating formula, 
namely that of Lagrange for the functional $\l(f) = f(x)$:
$$
\textstyle
     f(x) \approx p(x) = \sum_{i=0}^n L_i(x) f(x_i), \qquad
    L_i(x) = \prod_{j=0,j\ne i}^n \frac{(x-x_j)}{(x_i-x_j)}
$$
By linearity, $\l(p) = \sum_{i=0}^n \l(L_i) p(x_i)$, 
thus the interpolating formula has the form
$$
\textstyle
   \l(f) \approx \sum_{i=0}^n \l(L_i) f(x_i).
$$
It is exact on $\PP_{n}$, by definition. The converse is true as well. 
For if (8.1) is exact on $\PP_{n}$, take $f(x)=L_j(x)$ to obtain
$\l(L_j) = \sum_{i=0}^n a_i L_j(x_i) = a_j$. Hence, we have just proved 

{\bf Lemma.} The formula  $\l(f) \approx \sum_{i=0}^n a_i f(x_i)$
is interpolating $\Leftrightarrow$ it is exact for all $p\in\PP_{n}$. 

%Also, we have found the error of approximation 
%$e(f) := \l(f) - \l(p) = \frac{1}{(n+1)!} \omega(x) f^{(n+1)}(\xi)$



\subsection{Numerical integration}

A formula of numerical integration 
$$
\textstyle
    \l(f) := I(f) := \int_a^b f(x) w(x)\D x \approx 
    \sum_{i=0}^n a_i f(x_i),\qquad f\in C[a,b].
$$
is called a {\em quadrature formula} with the {\em nodes} $(x_i)$
and the {\em weights} $(a_i)$. For any $n+1$ fixed $x_i$s,
the interpolating formula $a_i = I(L_i)$ is exact on $\PP_{n}$.
Can one find a {\em quadrature of higher accuracy} 
which is exact for all $p\in\PP_m$ with some $m > n$?

{\bf Claim.} No quadrature formula with $n+1$ nodes
is exact for all $q\in\PP_m$ if $m > 2n+1$.\\
\hspace*{8pt}{\bf Proof.} 
Take $q(x) = \prod_{i=0}^n (x-x_i)^2 \ge 0$. Then $q\in\PP_{2n+2}$,
and $I(q) > 0$, but $\sum_{i=0}^n a_i q(x_i) = 0$ 
for any $a_i$s. Hence the integral and the quadrature do not match. \QED  

Our next goal is to show that $m=2n+1$ can be attained. For this we need
 
{\bf Lemma.} Let $p_{n+1}$ be the $(n+1)$-st orthogonal polynomial.
Then  all the zeros of $p_{n+1}$ are real,
distinct and lie in the interval $(a,b)$.\\
%
\hspace*{8pt}{\bf Proof.} 
Denote by $k$ the number of the sign  changes of $p_{n+1}$ in $(a,b)$ and
assume that $k \le n$. If $k=0$, set $q\equiv 1$. For $1\le k\le n$,
we let $q(x):=\prod_{j=1}^k (x-\xi_j)$ where $\xi_j$s are the points
where a sign change of $p_n$ occurs. 
Then $q\in\BB{P}_k$, $k\le n$, thus $\ip[p_{n+1},q]=0$.
On the other hand, by construction,
$q(x)p_{n+1}(x)$ does not change sign throughout $[a,b]$ 
and vanishes at a finite number of points, hence
$
   |\ip[q,p_{n+1}]| = \left|\int_a^b q(x)p_{n+1}(x)w(x)\D x\right|
=\int_a^b |q(x)p_{n+1}(x)|w(x)\D x>0,
$ 
a contradiction. Thus, $k=n+1$. \QED


{\bf Theorem.}  The interpolating quadrature with $n+1$ nodes
$(x_i)_{i=0}^n$ is exact for all $q\in\PP_{2n+1}$ iff
$(x_i)_{i=0}^n$ are the zeros of $(n+1)$-st orthogonal polynomial $p_{n+1}$.\\
%
\hspace*{8pt}{\bf Proof.} 
Let $(x_i)_{i=0}^n$ be the zeros of $p_{n+1}$. Given any $q\in\PP_{2n+1}$, 
we can represent it uniquely as $q(x)=r(x)p_{n+1}(x)+s(x)$, 
where $r,s\in\PP_{n}$. Since $p_{n+1}$ is orthogonal to $r$,
$$
\textstyle
   I(q) = \int_a^b q(x)w(x)\D x = \int_a^b s(x)w(x)\D x =  I(s).
$$
On the other hand, because $p_{n+1}(x_i) = 0$,
$$
\textstyle
   \sum_{i=0}^n a_i q(x_i) = \sum_{i=1}^n a_i s(x_i).
$$
But $s\in\PP_{n}$ and the interpolating quadrature is exact 
on $\PP_{n}$, hence $I(s) = \sum a_i s(x_i)$, and therefore
$I(q) = \sum a_i q(x_i)$ for all $q \in\PP_{2n+1}$.\\
%
\hspace*{8pt} Conversely, if a quadrature with $n+1$ nodes $(x_i)$ 
is exact for all $p\in\PP_{2n+1}$, then letting 
$p_{n+1}(x) := \prod_{i=0}^n (x-x_i)\in \PP_{n+1}$ 
and taking any $q\in\PP_{n}$ we find 
$$
\textstyle
   \int_a^b q(x) p_n(x) w(x) \D x
 = \sum_{i=0}^n a_i [q(x_i)p_n(x_i)] = 0
$$
i.e., $p_{n+1}$ is orthogonal to all $q\in\PP_{n}$. \QED  

{\small
{\bf Example.} For $[a,b]=[-1,1]$, $w(x)\equiv1$, the underlying
orthogonal polynomials are the {\em Legendre polynomials}.
The corresponding weights and nodes of Gaussian quadratures are 
as follows.
$$
\begin{array}{lll}
 P_1(x) = x  & x_0=0 & a_0 = 2;\\[1ex]
  P_2(x)=\frac{3}{2}(x^2-\frac{1}{3}) 
& x_0 = -\frac{1}{\sqrt{3}}, x_1 = \frac{1}{\sqrt{3}}
& a_0 = 1, a_1 = 1 \\[1ex]
  P_3(x)=\frac{5}{2}(x^3-\frac{3}{5}x) 
& x_0 = -\sqrt\frac{3}{5}, x_1 = 0, x_2 = \sqrt\frac{3}{5}
& a_0 =\frac{5}{9}, a_1 = \frac{8}{9}, a_2 =\frac{5}{9} 
\end{array}
$$
{\bf Example.} For $[a,b]=[-1,1]$, $w(x)=(1-x^2)^{-1/2}$, 
the orthogonal polynomials are the {\em Chebyshev polynomials}
and the qudrature rule is particularly simple
$$
\textstyle
 T_{n+1}(x) = \cos (n+1) \arccos x, \qquad  x_i=\cos\frac{2i+1}{2n+2}\pi, 
 \qquad a_i = \frac{2}{n}, \qquad i=0,\ldots,n.
$$
{\bf Example} ({\em Simplest quadrature formulae} ($w\equiv 1$)).
%
\begin{center}
\begin{tabular}{rcl}
 the rectangle rule:
   & \fbox{$I(f) \approx (b-a)f(a)$}
   & \small \begin{tabular}{l} 1-point, non-Gaussian \\ 
                               exact on constants 
            \end{tabular} \\[2ex]
  the midpoint rule:
   & \fbox{$I(f) \approx (b-a)f(\frac{a+b}{2})$}
   & \small \begin{tabular}{l} 1-point, Gaussian,\\
                         exact on linear f-ns 
            \end{tabular} \\[2ex]
the trapezoidal rule: 
   & \fbox{$I(f) \approx (b-a)[\frac{1}{2}f(a)+\frac{1}{2}f(b)]$}
   & \small \begin{tabular}{l} 2-point, non-Gaussian,\\
                               exact on linear f-ns 
            \end{tabular} \\[2ex]
   the Simpson rule:
   & \fbox{$I(f) \approx (b-a)
  [\frac{1}{6}f(a) + \frac{2}{3}f(\frac{a+b}{2}) + \frac{1}{6}f(b)]$}
   & \small \begin{tabular}{l} 3-point, non-Gaussian, \\
                               but of higher accuracy \\
                               (exact on cubics)
            \end{tabular} 
\end{tabular}
\end{center}
}

\end{document}


