%%  chapter-10.tex
%%  Approximation Theory: Chapter 10
%%  Neal L. Carothers
%%  Bowling Green State University
%%  Bowling Green, Ohio  43403
%%  carother@math.bgsu.edu
%%  http://www.bgsu.edu/~carother/


\input 680-setup.tex

\chaptertitle = {Gaussian Quadrature}

\centerline{\hfil\tf Gaussian Quadrature\hfil}
\vskip-\baselineskip
\line{\sc Math 680 \hfil 7/26/94}

\noindent
Numerical integration, or {\sl quadrature}, is the 
process of approximating the value of a definite integral
$\int_a^bf(x)\,w(x)\,dx$ based only on a finite number of 
values or ``samples'' of $f$ (much like a Riemann sum).
A {\sl linear quadrature formula\/} takes the form
$$\int_a^bf(x)\,w(x)\,dx \ \approx \ \sum_{k=1}^nA_kf(x_k),$$
where the {\sl nodes\/} $(x_k)$ and the {\sl weights\/}
$(A_k)$ are at our disposal.    

\noindent
{\bf Example}.
Consider the quadrature formula
$$I(f) \ = \ \int_{-1}^1f(x)\,dx \ \approx \ 
	{1\over n}\sum_{k=-n}^{n-1}f\left({{2k+1}\over{2n}}\right)
	\ = \ I_n(f).$$
If $f$ is continuous, then we clearly have $I_n(f)\to\int_{-1}^1f$
as $n\to\infty$. (Why?)  But in the particular case $f(x)=x^2$
we have (after some simplification)
$$I_n(f) \ = \ {1\over n}\sum_{k=-n}^{n-1}\left({{2k+1}\over{2n}}\right)^2
	= \ {1\over{2n^3}}\sum_{k=0}^{n-1}(2k+1)^2 \ = \ {2\over3}-{1\over{6n^2}}.$$
That is, $|\,I_n(f)-I(f)\,|=1/6n^2$.  In particular, we
would need
to take $n\ge130$ to get $1/6n^2\le10^{-5}$, for example,
and this would require that we perform over
$250$ evaluations of $f$.  We'd like a method that
converges a bit faster!  In other words, there's no 
shortage of quadrature formulas---we just want 
{\sl faster\/} ones. 

One reasonable requirement for our proposed quadrature 
formula is that it be {\sl exact for polynomials\/}
of low degree.  As it happens, this is easy to come by.

\proclaim Lemma 1.
Given $w(x)$ on $[\,a,b\,]$ and nodes $a\le x_1<\cdots<x_n\le b$, 
there exist unique weights $A_1,\ldots,A_n$ such that 
$$\int_a^bp(x)\,w(x)\,dx \ = \ \sum_{i=1}^nA_i\,p(x_i)$$
for all polynomials $p\in{\cal P}_{n-1}$.

\proof
Let $\ell_1,\ldots,\ell_n$ be the Lagrange interpolating polynomials
of degree $n-1$
associated to the nodes $x_1,\ldots,x_n$, and recall that we have
$p = \sum_{i=1}^np(x_i)\,\ell_i$ for all  $p\in{\cal P}_{n-1}$.
Hence,
$$\int_a^bp(x)\,w(x)\,dx \ = \ 
	\sum_{i=1}^np(x_i)\int_a^b\ell_i(x)\,w(x)\,dx.$$
That is, $A_i=\int_a^b\ell_i(x)\,w(x)\,dx$ works.
To see that this is the only choice, suppose that
$$\int_a^bp(x)\,w(x)\,dx \ = \ \sum_{i=1}^nB_i\,p(x_i)$$
is exact for all $p\in{\cal P}_{n-1}$, and set
$p=\ell_j$:
$$A_j \ = \ \int_a^b\ell_j(x)\,w(x)\,dx \ = \ 
	\sum_{i=1}^nB_i\,\ell_j(x_i) \ = \ B_j.
	\eqno\qed$$

The point here is that $\ell_1,\ldots,\ell_n$ form
a basis for ${\cal P}_{n-1}$ and integration is linear;
thus, integration is completely determined by its
action on the basis.  Said another way, the $n$ point
evaluations $\delta_i(p)=p(x_i)$ satisfy
${\cal P}_{n-1}\cap\left(\bigcap_{i=1}^n\ker\delta_i\right)=\{0\}$,
and it follows that every linear, real-valued function on 
${\cal P}_{n-1}$ must be a linear combination of
the $\delta_i$'s.  Here's why: ${\cal P}_{n-1}$ may be
identified with $\R^n$ by way of the isomorphism
$p\mapsto(p(x_1),\ldots,p(x_n))$.  A linear, real-valued
function on ${\cal P}_{n-1}$, must, then, correspond
to some linear, real-valued function on $\R^n$.  In
other words, it's given by inner product against some
fixed vector $(A_1,\ldots,A_n)$; in particular, we must
have $I(p)=\sum_{i=1}^nA_i\,p(x_i)$.

In any case, we now have our quadrature formula:  For 
$f\in C[\,a,b\,]$ we define 
$$I_n(f) \ = \ \sum_{i=1}^nA_i\,f(x_i) \ = \ 
	\int_a^b(L_{n-1}(f))(x)\,w(x)\,dx \ \approx \ 
	\int_a^bf(x)\,w(x)\,dx \ = \ I(f),$$
where $L_{n-1}$ is the Lagrange interpolating polynomial
of degree $n-1$ that agrees with $f$ at the nodes 
$x_1,\ldots,x_n$.  
Moreover, the formula is exact for $f\in{\cal P}_{n-1}$.

It's easy to give a bound on $|I_n(f)|$ in terms of
$\Vert f\Vert$; indeed,
$$|I_n(f)| \ \le \ \sum_{i=1}^n|A_i|\,|f(x_i)| \ 
	\le \Vert f\Vert\left(\sum_{i=1}^n|A_i|\right).$$
In other words, $\lambda_n=\sum_{i=1}^n|A_i|$ are the
``Lebesgue numbers'' for this process.  As with all
previous settings, we want these numbers to be
uniformly bounded.

If $w(x)\equiv1$ and if $f$ is $n$-times continuously differentiable, 
we even have an error estimate for our quadrature formula: 
$$\left|\int_a^bf-\int_a^bL_{n-1}(f)\right| \ 
	\le \ \int_a^b|f-L_{n-1}(f)|
	\le \ {1\over{n!}}\,\Vert f^{\,(n)}\Vert
		\int_a^b\prod_{i=1}^n|x-x_i|\,dx$$
(recall the Theorem on page 5 of 
``A Brief Introduction to Interpolation'').
As it happens, the integral on the right is minimized when the
$x_i$'s are taken to be the zeros of the Chebyshev
polynomial $U_n$ (see Rivlin, page 72).

The fact that a quadrature formula is exact for
polynomials of low degree does not, by itself,
guarantee that the formula is highly accurate.
The problem is that $\sum_{i=1}^nA_i\,f(x_i)$
may be estimating a very small quantity through
the cancellation of very large quantities.  So,
for example, a positive function may yield a negative
result in this approximate integral.  This wouldn't
happen if the $A_i$'s were all positive---and we've
already seen how useful positivity can be.  Our
goal here is to further 
improve our quadrature formula to have this property.
But we have yet to take advantage of the fact that the
$x_i$'s are at our disposal.  We'll let Gauss show us
the way!  

\proclaim Theorem. {\rm(Gauss)} \ 
Fix a weight $w(x)$ on $[\,a,b\,]$, and let $(Q_n)$ be
the canonical sequence of orthogonal polynomials relative
to $w$.  Given $n$, let $x_1,\ldots,x_n$ be the zeros
of $Q_n$ {\rm(}these all lie in $(a,b)$\/{\rm)}, 
and choose $A_1,\ldots,A_n$ so that the
formula $\sum_{i=1}^nA_if(x_i)\approx\int_a^bf(x)\,w(x)\,dx$
is exact for polynomials of degree less than $n$.
Then, in fact, the formula is exact for all
polynomials of degree less than $2n$.

\proof
Given a polynomial $P$ of degree less than $2n$, we may
divide: $P=Q_nR+S$, where $R$ and $S$ are polynomials
of degree less than $n$.  Thus,
$$\eqalign{
\int_a^bP(x)\,w(x)\,dx \ 
	&= \ \int_a^bQ_n(x)\,R(x)\,w(x)\,dx \ + \ \int_a^bS(x)\,w(x)\,dx\cr
	&= \ \int_a^bS(x)\,w(x)\,dx,\quad\hbox{since deg } R<n\cr
	&= \ \sum_{i=1}^nA_iS(x_i),\quad\hbox{since deg } S<n.\cr
}$$
But $P(x_i)=Q_n(x_i)\,R(x_i)+S(x_i)=S(x_i)$, since $Q_n(x_i)=0$.
Hence, $\int_a^bP(x)\,w(x)\,dx=\sum_{i=1}^nA_iP(x_i)$ for all
polynomials $P$ of degree less than $2n$.~\qed

Amazing!  But, well, not really: ${\cal P}_{2n-1}$ is of
dimension $2n$, and we had $2n$ numbers $x_1,\ldots,x_n$
and $A_1,\ldots,A_n$ that we could choose as we saw fit.
Said another way, the division algorithm tells us that
${\cal P}_{2n-1}\approx Q_n{\cal P}_{n-1}\oplus{\cal P}_{n-1}$.
Since $Q_n{\cal P}_{n-1}\subset\ker(I_n)$, the action of $I_n$
on ${\cal P}_{2n-1}$ is the same as its action on a ``copy''
of ${\cal P}_{n-1}$.

Also not surprising is that this particular choice of
$x_i$'s is unique.

\proclaim Lemma 2.
Suppose that $a\le x_1<\cdots<x_n\le b$ and $A_1,\ldots,A_n$
are given so that the equation $\int_a^bP(x)\,w(x)\,dx=\sum_{i=1}^nA_iP(x_i)$ 
is satisfied for all polynomials $P$ of degree less than $2n$.  Then,
$x_1,\ldots,x_n$ are the zeros of $Q_n$.

\proof
Let $Q(x)=\prod_{i=1}^n(x-x_i)$.  Then, for $k<n$, the polynomial
$Q\cdot Q_k$ has degree $n+k<2n$.  Hence,
$$\int_a^bQ(x)\,Q_k(x)\,w(x)\,dx \ = \ \sum_{i=1}^nA_iQ(x_i)\,Q_k(x_i)=0.$$
Since $Q$ is a monic polynomial of degree $n$ which is orthogonal
to each $Q_k$, $k<n$, we must
have $Q=Q_n$.  Thus, the $x_i$'s are actually the zeros of $Q_n$.~\qed

According to Rivlin, the phrase {\sl Gaussian quadrature\/} 
is usually reserved for the specific 
quadrature formula whereby $\int_{-1}^1f(x)\,dx$ is
approximated by $\int_{-1}^1(L_{n-1}(f))(x)\,dx$, where
$L_{n-1}(f)$ is the Lagrange interpolating polynomial 
to $f$ using the zeros of the $n$-th Legendre polynomial
as nodes.  (What a mouthful!)  What is actually being
described in our version of Gauss's theorem is
{\sl Gaussian-type quadrature}.  

Before computers, Gaussian quadrature was little more
than a curiosity; the roots of $Q_n$ are typically
irrational, and certainly not easy to come by.  By
now, though, it's considered a standard quadrature 
technique.  In any case, we still can't judge the
quality of Gauss's method without a bit more information.

\filbreak

\noindent
{\bf Gaussian-type Quadrature}

\noindent
First, let's summarize our rather cumbersome notation.

$$\vbox{\halign{%
\hfil$#$\hfil & \qquad$#$\hfil & \quad\qquad$#$\hfil & \quad\hfil$#$\hfil\cr
\hbox{orthogonal} &&& \hbox{approximate}\cr
\noalign{\vskip-6true pt}
\hbox{polynomial} & \hbox{zeros} & \hbox{weights} & \hbox{integral}\cr
\noalign{\smallskip}
Q_1 & x_1^{(1)} & A_1^{(1)} & I_1 \cr
Q_2 & x_1^{(2)}, x_2^{(2)} & A_1^{(2)}, A_2^{(2)} & I_2 \cr
Q_3 & x_1^{(3)}, x_2^{(3)}, x_3^{(3)} & A_1^{(3)}, A_2^{(3)}, A_3^{(3)} & I_3 \cr
\vdots & \qquad\vdots & \qquad\vdots & \vdots \cr
}}$$
Hidden here is the Lagrange interpolation formula
$L_{n-1}(f)=\sum_{i=1}^nf(x_i^{(n)})\,\ell_i^{(n-1)}$,
where $\ell_i^{(n-1)}$ denote the Lagrange polynomials 
of degree $n-1$ based on $x_1^{(n)},\ldots, x_n^{(n)}$.
The $n$-th quadrature formula is then
$$I_n(f) \ = \ \int_a^bL_{n-1}(f)(x)\,w(x)\,dx \ 
	= \ \sum_{i=1}^nA_i^{(n)}f(x_i^{(n)}) \ 
	\approx \ \int_a^bf(x)\,w(x)\,dx,$$
which is exact for polynomials of degree
less than $2n$.

By way of one example, Hermite showed that 
$A_k^{(n)}=\pi/n$ for the Chebyshev weight $w(x)=(1-x^2)^{-1/2}$
on $[-1,1\,]$.  Remarkably, $A_k^{(n)}$ doesn't depend
on $k$!  The quadrature formula in this case reads:
$$\int_{-1}^1{{f(x)\,dx}\over\sqrt{1-x^2}} \ \approx \ 
	{\pi\over n}\,\sum_{k=1}^nf\left(\cos\,{{2k-1}\over{2n}}\,\pi\right).$$
You can find full details in Natanson's 
{\it Constructive Function Theory}, Vol.\ III.

The key result, due to Stieltjes, is that $I_n$ is {\sl positive\/}:

\proclaim Lemma 3.
$A_1^{(n)},\ldots, A_n^{(n)}>0$ \ and \  
$\sum_{i=1}^nA_i^{(n)}=\int_a^bw(x)\,dx$.

\proof
The second assertion is obvious (just take $f\equiv1$\/).
For the first, fix $1\le j\le n$ and notice that 
$(\ell_j^{(n-1)})^2$ is of degree $2(n-1)<2n$.  Thus,
$$0 \ < \ \langle\,\ell_j^{(n-1)},\ell_j^{(n-1)}\,\rangle \ 
	= \ \int_a^b\left[\ell_j^{(n-1)}(x)\right]^2w(x)\,dx \ = \ 
	\sum_{i=1}^nA_i^{(n)}\left[\ell_j^{(n-1)}(x_i^{(n)})\right]^2
	\ = \ A_j^{(n)},$$
because $\ell_j^{(n-1)}(x_i^{(n)})=\delta_{i,j}$.~\qed  


Now our last calculation is quite curious; what we've shown is that
$$A_j^{(n)} \ = \ \int_a^b\ell_j^{(n-1)}(x)\,w(x)\,dx
	\ = \ \int_a^b\left[\ell_j^{(n-1)}(x)\right]^2w(x)\,dx.$$

The same calculation as above also proves

\proclaim Corollary.
$\langle\,\ell_i^{(n-1)},\ell_j^{(n-1)}\,\rangle=0$ for
$i\ne j$.

Since $A_1^{(n)},\ldots, A_n^{(n)}>0$, it follows that 
$I_n(f)$ is positive; that
is, $I_n(f)\ge0$ whenever $f\ge0$.  The second assertion
in Lemma 3 tells us that the $I_n$'s are {\sl uniformly bounded\/}:
$$|I_n(f)| \ \le \ \Vert f\Vert\,\sum_{i=1}^nA_i^{(n)}
	\ = \ \Vert f\Vert\int_a^bw(x)\,dx,$$
and this is the same bound that holds for 
$I(f)=\int_a^bf(x)\,w(x)\,dx$ itself.  Given all of this,
proving that $I_n(f)\to I(f)$ is a piece of cake.  The
following result is again due to Stieltjes (\'a la Lebesgue).

\proclaim Theorem.
In the above notation, $|I_n(f)-I(f)|\le 
	2\left(\int_a^bw(x)\,dx\right)E_{2n-1}(f)$.
In particular, $I_n(f)\to I(f)$
for evey $f\in C[\,a,b\,]$.

\proof
Let $p^*$ be the best uniform approximation to $f$ out of
${\cal P}_{2n-1}$.  Then, since $I_n(p^*)=I(p^*)$, we have 
$$\openup1\jot\eqalignno{
|I(f)-I_n(f)| \ &\le \ |I(f-p^*)| \ + \ |I_n(f-p^*)|\cr
	&\le \ \Vert f-p^*\Vert\int_a^bw(x)\,dx \ + \ 
		\Vert f-p^*\Vert\,\sum_{i=1}^nA_i^{(n)}\cr
	&= \ 2\,\Vert f-p^*\Vert\int_a^bw(x)\,dx \ 
		= \ 2E_{2n-1}(f)\int_a^bw(x)\,dx.&\qed\cr
}$$

%\filbreak

\noindent
{\bf Computational Considerations}

\noindent
You've probably been asking yourself: ``How do I find
the $A_i$'s without integrating?''  Well, first let's
recall the definition:  In the case of Gaussian-type
quadrature we have 
$$A_i \ = \ \int_a^b\ell_i(x)\,w(x)\,dx \ 
	= \ \int_a^b{{Q_n(x)}\over{(x-x_i)\,Q_n'(x_i)}}\,w(x)\,dx$$
(because ``$W$'' is the same as $Q_n$ here---the $x_i$'s
are the zeros of $Q_n$\/).  Next, consider the
function
$$\varphi_n(x) \ = \ 
	\int_a^b{{Q_n(t)-Q_n(x)}\over{t-x}}\,w(t)\,dt.$$
Since $t-x$ divides $Q_n(t)-Q_n(x)$, note 
that $\varphi_n$ is actually a {\sl polynomial\/}
and that 
$$\varphi_n(x_i) \ = \ 
	\int_a^b{{Q_n(t)}\over{t-x_i}}\,w(t)\,dt \ 
	= \ A_iQ_n'(x_i).$$
Now $Q_n'(x_i)$ is readily available; we just need to
compute $\varphi_n(x_i)$.

\proclaim Claim.
The $\varphi_n$'s satisfy the same recurrence formula
as the $Q_n$'s
$$\varphi_{n+1}(x)=(x-a_n)\varphi_n(x)-b_n\varphi_{n-1}(x),
	\quad n\ge1,$$
but with different starting values
$$\varphi_0(x)\equiv0,\qquad\hbox{and}\qquad
	\varphi_1(x)\equiv\int_a^bw(x)\,dx.$$

\proof
The formulas for $\varphi_0$ and $\varphi_1$ are
obviously correct, since $Q_0(x)\equiv1$ and 
$Q_1(x)=x-a_0$.  We only need to check the recurrence
formula itself.
$$\openup2\jot\eqalign{
\qquad\varphi_{n+1}(x)  &=  
	\int_a^b{{Q_{n+1}(t)-Q_{n+1}(x)}\over{t-x}}\,w(t)\,dt\cr
	&=  \int_a^b{{(t-a_n)\,Q_n(t)-b_nQ_{n-1}(t)
		-(x-a_n)\,Q_n(x) + b_nQ_{n-1}(x)}\over{t-x}}\,w(t)\,dt\cr
	&=  (x-a_n)\int_a^b{{Q_n(t)-Q_n(x)}\over{t-x}}\,w(t)\,dt 
		\ - \ b_n\int_a^b{{Q_{n-1}(t)-Q_{n-1}(x)}\over{t-x}}\,w(t)\,dt\cr
	&=(x-a_n)\,\varphi_n(x)-b_n\,\varphi_{n-1}(x),\cr
}$$
since $\int_a^bQ_n(t)\,w(t)\,dt=0$.~\qed

Of course, the derivatives $Q_n'$ satisfy a  recurrence
relation of sorts, too:
$$Q_{n+1}'(x) \ = \ Q_n(x) + (x-a_n)\,Q_n'(x) - b_n\,Q_{n-1}'(x).$$
But $Q_n'(x_i)$ can be computed without knowing $Q_n'(x)$.
Indeed, since $Q_n(x)=\prod_{i=1}^n(x-x_i)$, we have 
$Q_n'(x_i)=\prod_{j\ne i}(x_i-x_j)$.

The weights $A_i^{(n)}$, or {\sl Christoffel numbers}, 
together with the zeros of $Q_n$ are tabulated in a
variety of standard cases.  See, for example,
{\it Handbook of Mathematical Functions with Formulas,
Graphs, and Tables}, by Abramowitz and Stegun, eds.
In practice, of course, it's enough to tabulate data
for the case $[\,a,b\,]=[-1,1\,]$.  

%\filbreak

\noindent
{\bf Applications to Interpolation}

\noindent
Although $L_n(f)$ isn't typically a good uniform
approximation to $f$, if we interpolate at the
zeros of an orthogonal polynomial $Q_{n+1}$, then $L_n(f)$
will be a good approximation in the $\Vert\cdot\Vert_1$ or
$\Vert\cdot\Vert_2$ 
norm generated by the corresponding weight $w$.
Specifically, by rewording our earlier results, it's easy to get
estimates for each of the errors $\int_a^b|f-L_n(f)|\,w$ 
and $\int_a^b|f-L_n(f)|^2\,w$.
We use essentially the same notation as before, except
now we take
$$L_n(f) \ = \ \sum_{i=1}^{n+1}f\bigl(x_i^{(n+1)}\bigr)\,\ell_i^{(n)},$$
where $x_1^{(n+1)},\ldots,x_{n+1}^{(n+1)}$ are the roots of 
$Q_{n+1}$ and $\ell_i^{(n)}$ is of degree $n$.  This leads
to a quadrature formula that's exact on polynomials of
degree less than $2(n+1)$.

As we've
already seen, $\ell_1^{(n)},\ldots,\ell_{n+1}^{(n)}$ are 
{\sl orthogonal\/} and so $\Vert L_n(f)\Vert_2$ may be
computed exactly.

\proclaim Lemma.
$\Vert L_n(f)\Vert_2 \ \le \  
	\Vert f\Vert\left(\int_a^bw(x)\,dx\right)^{1/2}$.

\proof
Since $L_n(f)^2$ is a polynomial of degree $\le2n<2(n+1)$,
we have 
$$\openup1\jot\eqalignno{
\Vert L_n(f)\Vert_2^2 \ &= \ \int_a^b[L_n(f)]^2\,w(x)\,dx\cr
	&= \ \sum_{j=1}^{n+1}A_j\left[\sum_{i=1}^{n+1}
		f\bigl(x_i^{(n+1)}\bigr)\,\ell_i^{(n)}\bigl(x_j^{(n+1)}\bigr)\right]^2\cr
	&= \ \sum_{j=1}^{n+1}A_j\left[f\bigl(x_j^{(n+1)}\bigr)\right]^2 \ 
	\le \ \Vert f\Vert^2\sum_{j=1}^{n+1}A_j \ 
		= \ \Vert f\Vert^2\int_a^bw(x)\,dx.&\qed\cr
}$$

As usual, once we have an estimate for the norm of an 
operator, we also have an analogue of Lebesgue's theorem.

\proclaim Theorem.
$\Vert f-L_n(f)\Vert_2 \ \le \  
	2\left(\int_a^bw(x)\,dx\right)^{1/2}E_n(f)$.

\proof
Here we go again!  Let $p^*$ be the best uniform
approximation to $f$ out ${\cal P}_n$ and use
the fact that $L_n(p^*)=p^*$.
$$\eqalignno{
\Vert f-L_n(f)\Vert_2 \ &\le \ \Vert f-p^*\Vert_2+\Vert L_n(f-p^*)\Vert_2\cr
	&\le \ \Vert f-p^*\Vert\left(\int_a^bw(x)\,dx\right)^{1/2}
	+ \ \Vert f-p^*\Vert\left(\int_a^bw(x)\,dx\right)^{1/2}\cr
	&= \ 2E_n(f)\left(\int_a^bw(x)\,dx\right)^{1/2}.&\qed\cr
}$$

Hence, if we interpolate $f\in C[\,a,b\,]$
at the zeros of $(Q_n)$, then 
$L_n(f)\to f$ in $\Vert\cdot\Vert_2$ norm.
The analogous result for the $\Vert\cdot\Vert_1$ 
norm is now easy:

\proclaim Corollary.
$\int_a^b|f(x)-L_n(f)(x)|\,w(x)\,dx \ \le \ 
	2\left(\int_a^bw(x)\,dx\right)E_n(f)$.

\proof
We apply the Cauchy-Schwarz inequality:
$$\eqalignno{
\int_a^b|f(x)-L_n(f)(x)|\,w(x)\,dx \ 
	&= \ \int_a^b|f(x)-L_n(f)(x)|\,\sqrt{w(x)}\,\sqrt{w(x)}\,dx\cr
	&\le \ \left(\int_a^b|f(x)-L_n(f)(x)|^2\,w(x)\,dx\right)^{1/2}
		\left(\int_a^bw(x)\,dx\right)^{1/2}\cr
	&\le \ 2E_n(f)\int_a^bw(x)\,dx.&\qed\cr
}$$

Essentially the same device allows an estimate of
$\int_a^bf(x)\,dx$ in terms of $\int_a^bf(x)\,w(x)\,dx$
(which may be easier to compute).

\proclaim Corollary.
If $\int_a^bw(x)^{-1}\,dx$ is finite, then 
$$\eqalignno{
\int_a^b|f(x)-L_n(f)(x)|\,dx \ 
	&= \ \int_a^b|f(x)-L_n(f)(x)|\,\sqrt{w(x)}\,{1\over{\sqrt{w(x)}}}\,dx\cr
	&\le \ \left(\int_a^b|f(x)-L_n(f)(x)|^2\,w(x)\,dx\right)^{1/2}
		\left(\int_a^b{1\over{w(x)}}\,dx\right)^{1/2}\cr
	&\le \ 2E_n(f)\left(\int_a^bw(x)\,dx\right)^{1/2}
		\left(\int_a^b{1\over{w(x)}}\,dx\right)^{1/2}.&\qed\cr
}$$

In particular, the Chebyshev weight satisfies 
$$\int_{-1}^1{{dx}\over{\sqrt{1-x^2}}}=\pi\qquad\hbox{and}\qquad
	\int_{-1}^1\sqrt{1-x^2}\,dx={\pi\over2}.$$
Thus, interpolation at the zeros of the Chebyshev polynomials
(of the first kind) would provide good, {\sl simultaneous\/} approximation
in each of the norms $\Vert\cdot\Vert_1$, $\Vert\cdot\Vert_2$, 
and $\Vert\cdot\Vert$.

%\filbreak

\noindent
{\bf The Moment Problem}

\noindent
Given a positive, continuous weight function $w(x)$ on $[\,a,b\,]$,
the number
$$\mu_k \ = \ \int_a^bx^k\,w(x)\,dx$$
is called the $k$-th {\sl moment\/} of $w$.  In physical terms,
if we think of $w(x)$ as the density of a thin rod placed
on the interval $[\,a,b\,]$, then $\mu_0$ is the mass of the
rod, $\mu_1/\mu_0$ is its center of mass, $\mu_2$ is its
moment of inertia (about $0$\/), and so on.  In probabilistic
terms, if $\mu_0=1$, then $w$ is the probability density function
for some random variable, $\mu_1$ is the expected value, or mean,
of this random variable, and $\mu_2-\mu_1^2$ is its variance.
The moment problem (or problems, really) concern the inverse
procedure.  What can be measured in real life are the moments---can
the moments be used to find the density function?

\item{}
{\bf Questions}: Do the moments determine $w$?  Do different 
weights have different moment sequences?  If we knew
the sequence $(\mu_k)$, could we find $w$?  How do we tell if 
a given sequence $(\mu_k)$ is the moment sequence for some
positive weight?  Do ``special'' weights give rise to ``special''
sequences?  

Now we've already answered one of these questions: The Weierstrass
theorem tells us that different weights have different moment
sequences.  Said another way, if 
$$\int_a^bx^k\,w(x)\,dx=0\qquad\hbox{for all}\ k=0,1,2,\ldots,$$
then $w\equiv0$.  Indeed, by linearity, this says that 
$\int_a^bp(x)\,w(x)\,dx=0$ for all polynomials $p$ which, in 
turn, tells us that $\int_a^bw(x)^2\,dx=0$.  (Why?)  The 
remaining questions are harder to answer.  We'll settle 
for simply stating a few pertinent results.

Given a sequence of numbers $(\mu_k)$, we define the $n$-th
difference sequence $(\Delta^n\mu_k)$ by
$$\eqalign{
\Delta^0\mu_k&=\mu_k\cr
\Delta^1\mu_k&=\mu_k-\mu_{k+1}\cr
\Delta^n\mu_k&=\Delta^{n-1}\mu_k-\Delta^{n-1}\mu_{k+1},
	\quad n\ge1.\cr
}$$
For example, $\Delta^2\mu_k=\mu_k-2\mu_{k+1}+\mu_{k+2}$.
More generally, induction will show that 
$$\Delta^n\mu_k=\sum_{i=0}^n(-1)^i{n\choose i}\mu_{k+i}.$$
In the case of a weight $w$ on the interval $[\,0,1\,]$, this sum
is easy to recognize as an integral.  Indeed,
$$\int_0^1x^k(1-x)^n\,w(x)\,dx
	=\sum_{i=0}^n(-1)^i{n\choose i}\int_0^1x^{k+i}\,w(x)\,dx
	=\sum_{i=0}^n(-1)^i{n\choose i}\mu_{k+i}.$$
In particular, if $w$ is nonnegative, then we must have 
$\Delta^n\mu_k\ge0$ for every $n$ and $k$.  This 
observation serves as motivation for 

\proclaim Theorem.
The following are equivalent:
\item{\rm(a)} \ 
$(\mu_k)$ is the moment sequence of some nonnegative
weight function $w$ on $[\,0,1\,]$.
\item{\rm(b)} \ 
$\Delta^n\mu_k\ge0$ for every $n$ and $k$.
\item{\rm(c)} \ 
$a_0\mu_0+a_1\mu_1+\cdots+a_n\mu_n\ge0$
whenever $a_0+a_1x+\cdots+a_nx^n\ge0$ for
all\/ $0\le x\le1$.

The equivalence of (a) and (b) is due to Hausdorff.
A real sequence satisfying (b) or (c) is sometimes
said to be {\sl positive definite}.  

Now dozens of mathematicians worked on various aspects
of the moment problem: Chebyshev, Markov, Stieltjes,
Cauchy, Riesz, Fr\'echet, and on and on.  And several
of them, in particular Cauchy and Stieltjes, noticed
the importance of the integral $\int_a^b{{w(t)}\over{x-t}}\,dt$ 
in attacking the problem (compare this to Cauchy's 
integral formula).  It was Stieltjes, however, who gave
the first complete solution to such a problem---developing
his own integral (by considering $\int_a^b{{dW(t)}\over{x-t}}$), 
his own variety of continued fractions,
and planting the seeds for the study of orthogonal polynomials
while he was at it!  We will attempt to at least sketch
a few of these connections.  

To begin, let's fix our notation: To simplifiy things, we
suppose that we're given a nonnegative weight $w(x)$ on
a {\sl symmetric\/} interval $[-a,a\,]$, and that all of the 
moments of $w$ are finite.  We'll stick to our usual notations
for $(Q_n)$, the Gaussian-type quadrature formulas, and so on.  

\proclaim Lemma.
If\/ $x\notin[-a,a\,]$, then 
$\displaystyle\int_{-a}^a{{w(t)}\over{x-t}}\,dt \ 
	= \ \sum_{k=0}^\infty{{\mu_k}\over{x^{k+1}}}$.

\proof
$\displaystyle {1\over{x-t}}={1\over x}\cdot{1\over{1-(t/x)}}
	=\sum_{k=0}^\infty{{t^k}\over{x^{k+1}}}$, and the sum
converges uniformly because $|t/x|\le a/|x|<1$.  Now just
multiply by $w(t)$ and integrate.~\qed

By way of an example, consider the Chebyshev weight
$w(x)=(1-x^2)^{-1/2}$ on $[-1,1\,]$.  For $x>1$ we have 
$$\openup2\jot\eqalign{
\int_{-1}^1{{dt}\over{(x-t)\sqrt{1-t^2}}} \ 
	&= \ {\pi\over\sqrt{x^2-1}}
			\qquad\qquad\bigl(\hbox{set}\ t=2u/(1+u^2)\bigr)\cr
	&={\pi\over x}\left(1-{1\over{x^2}}\right)^{-1/2}\cr
	&={\pi\over x}\left[1+{1\over2}\cdot{1\over x^2}+
		{{1\cdot3}\over{2\cdot2}}\cdot{1\over{2!}}\cdot{1\over{x^4}}+\cdots\ 
		\right],\cr
}$$
using the binomial formula.  Thus, we've found all the moments:
$$\openup1\jot\eqalign{
\mu_0 \ &= \ \int_{-1}^1{{dt}\over\sqrt{1-t^2}} \ = \ \pi\cr
\mu_{2n-1} \ &= \ \int_{-1}^1{{t^{2n-1}dt}\over\sqrt{1-t^2}} \ = \ 0\cr
\mu_{2n} \ &= \ \int_{-1}^1{{t^{2n}dt}\over\sqrt{1-t^2}} \ = \ 
		{{1\cdot3\cdot5\cdots(2n-1)}\over{2^nn!}}\,\,\pi.\cr
}$$

Stieltjes proved much more: The integral $\int_{-a}^a{{w(t)}\over{x-t}}\,dt$
is actually an analytic
function of $x$ in $\C\setminus[-a,a\,]$.  In any case, since 
$x\notin[-a,a\,]$, we know that ${1\over{x-t}}$ is continuous on 
$[-a,a\,]$.  In particular, we can apply our quadrature formulas
(and Stieltjes theorem, p.\ 6) to write 
$$\int_{-a}^a{{w(t)}\over{x-t}}\,dt \ 
	= \ \lim_{n\to\infty}\sum_{i=1}^n{{A_i^{(n)}}\over{x-x_i^{(n)}}},$$
and these sums are recognizable:

\proclaim Lemma. \ 
$\displaystyle \sum_{i=1}^n{{A_i^{(n)}}\over{x-x_i^{(n)}}} \ = \ 
	{{\varphi_n(x)}\over{Q_n(x)}}.$

\proof
Since $\varphi_n$ has degree $<n$ and $\varphi_n(x_i^{(n)})\ne0$
for any $i$, we may appeal to partial-fractions to write
$${{\varphi_n(x)}\over{Q_n(x)}} \ = \ 
	{{\varphi_n(x)}\over{(x-x_1^{(n)})\cdots(x-x_n^{(n)})}} 
	\ = \ \sum_{i=1}^n{{c_i}\over{x-x_i^{(n)}}}$$
where $c_i$ is given by
$$c_i \ = \ 
	\left.{{\varphi_n(x)}\over{Q_n(x)}}\,(x-x_i^{(n)})\right]_{x=x_i^{(n)}}
	\ = \ {{\varphi_n(x_i^{(n)})}\over{Q_n'(x_i^{(n)})}}
	\ = \ A_i^{(n)}.\eqno\qed$$

Now here's where the continued fractions come in: Stieltjes
recognized the fact that 
$${{\varphi_{n+1}(x)}\over{Q_{n+1}(x)}} \ = \ 
	{b_0\over{\textstyle (x-a_0)-{\textstyle b_1\over
		{\textstyle (x-a_1)-\quad\lower2ex\hbox{$\ddots$}
			\hphantom{\strut\textstyle b_n\over\textstyle{(x-a_n)}}}}
		\atop
	{\hphantom{\qquad\qquad\textstyle (x-a_1)-\ddots}\textstyle-
		{{\textstyle b_n\over\textstyle{(x-a_n)}}}}}}$$
%% Zowie!!
(which can be proved by induction).  More generally, induction 
will show that the
$n$-th {\sl convergent\/} of a continued fraction can be
written as 
$${A_n\over B_n} \ = \ 
	{p_1\over{\textstyle q_1-{\textstyle p_2\over
		{\strut\textstyle q_2-\quad\lower2ex\hbox{$\ddots$}
			\hphantom{\qquad\textstyle b_n}}}
		\atop
	{\hphantom{\qquad\quad\textstyle a_1-\ddots}\textstyle-
		{{\textstyle p_n\over\textstyle q_n}}}}}$$
by means of the {\sl recurrence formulas\/}
$$\eqalign{
A_0 &= 0\cr
A_1 &= p_1\cr
A_n &= q_nA_{n-1}+p_nA_{n-2}\cr
}\qquad\qquad
\eqalign{
B_0 &= 1\cr
B_1 &= q_1\cr
B_n &= q_nB_{n-1}+p_nB_{n-2}\cr
}$$
where $n=2,3,4,\ldots$.  Please note that $A_n$ and $B_n$
satisfy the same recurrence formula, but with different
starting values (as is the case with $\varphi_n$ and $Q_n$\/).

Again using the Chebyshev weight as an example, for $x>1$ we have
$${\pi\over\sqrt{x^2-1}} \ = \ 
\int_{-1}^1{{dt}\over{(x-t)\sqrt{1-t^2}}} \ = \ 
	{\pi\over{x-{{\strut\scriptstyle 1/2}\over
		{\textstyle x-{{\strut\scriptstyle 1/4}\over
			{\textstyle x-{{\strut\scriptstyle 1/4}\over\ddots}}}}}}}
$$
since $a_n=0$ for all $n$, $b_1=1/2$, and $b_n=1/4$ for $n\ge2$.
In other words, we've just found a continued fraction expansion
for $(x^2-1)^{-1/2}$. 




\bye


%%  end of chapter-10.tex