\magnification=\magstep1
\centerline{Mathematics 371 -- Numerical Analysis}
\centerline{Lab 4 -- Adaptive Quadrature and Gaussian Quadrature}
\centerline{November 12, 2001}
\bigskip
\noindent
{\it Goals}
\bigskip
In this lab we will use the {\it adaptive quadrature\/}
methods introduced in 
class last week, and continue the development of {\it Gaussian quadrature\/}.
\bigskip
\noindent
{\it Adaptive Quadrature}
\bigskip
\item{A)}  The implementation of the recursive adaptive quadrature procedure
{\tt Adapt} we discussed in class on Wednesday, November 7 is available
on the course homepage.  Download it, save to a file in 
your home directory, and read it into
your Maple session.  Use it to approximate the integrals:
$$ \int_0^2 \sin(x^3)\ dx, \qquad \int_1^{10} {1 - \cos(x)\over x^2}\ dx$$
accurate to $10^{-5}$.  Compare with Maple's approximate value
obtained using a command of the form
\medskip
\centerline{\tt evalf(Int(f(x),x = a..b);}
\medskip
\item{B)}  (For after the lab) Draw the ``recursion tree'' of calls to 
the {\tt Adapt} procedure on either one of the integrals from part A -- your
choice -- by having the procedure print out the tracing messages.  
Discuss what you see there relative to the graph of the appropriate function
on the interval of integration.
\bigskip
\noindent
{\it Gaussian Quadrature}
\bigskip
Recall from class last Friday that the Gaussian quadrature
method on an interval $[a,b]$ depends on knowing an {\it orthogonal
basis\/} for the vector space of polynomial functions on $[a,b]$
of degree $\le n$ for some $n$ -- that is, a collection of polynomials
$$\{p_0(x), p_1(x), \ldots, p_n(x)\}$$
with the properties that
\medskip
\item{a)}  (orthogonality)
$$\langle p_i,p_j\rangle = \int_a^b p_i(x) p_j(x)\ dx = \cases{ 0 & if $i \ne j$, and\cr
                                       c > 0 & if $i = j$\cr}$$
\item{b)}  $\deg p_i(x) = i$ for all $i = 0, \ldots, n$.  
\item{c)}  The equation $p_i(x) = 0$ has $i$ distinct real roots
in the interval $[a,b]$ for all $i = 0, \ldots, n$.  
\medskip
\noindent
Knowing such a collection of polynomials, and writing $x_1, \ldots, x_n$
for the roots of the equation $p_n(x) = 0$ on $[a,b]$, 
we showed that the numerical integration rule:
$$\int_a^b f(x)\ dx = \sum_{j=1}^n c_j f(x_j)\leqno(1)$$
(where 
$$c_j = \int_a^b {\prod_{k\ne j}(x - x_k) \over \prod_{k\ne j}(x_j - x_k)}\ dx
\leqno(2)$$
for each $j = 1, \ldots n$) has degree of precision $2n - 1$(!)
\bigskip
\item{C)}  In this problem, you will use Maple to construct a collection
of such polynomials on $[a,b] = [-1,1]$ ``from scratch.''
\itemitem{1)}  Check that $p_0(x) = 1$ and $p_1(x) = x$ are orthogonal
under the inner product 
$$\langle f, g\rangle = \int_{-1}^1 f(x)g(x)\ dx.$$
Note that $p_0(1) = p_1(1) = 1$ also.
\itemitem{2)}  Use the following outline to construct a polynomial $p_2$
of degree 2 that is orthogonal to both $p_0$ and $p_1$.  First, 
let $p_2(x) = ax^2 + bx + c$ with undetermined (variable) coefficients.
Use Maple to compute $\int_{-1}^1 p_0(x)p_2(x)\ dx$ and 
$\int_{-1}^1 p_1(x)p_2(x)\ dx$.  Setting your results equal to zero
yields two linear equations for $a,b,c$.  Solve these together with
the equation obtained from $p_2(1) = 1$ to determine $a,b,c$.  Maple
note:  To solve a system of linear equations in Maple, you can use
a command of the following form.
\medskip
\centerline{\tt solve(eqns,vars);}
\medskip
\itemitem{}
where {\tt eqns} is the set of equations to be solved, and
{\tt vars} is the set of variables you want to solve for.
For example, to solve 
$$\eqalign{3x - 2y &= 1\cr
           4x + 5y &= 7\cr}$$
you can use
\medskip
\centerline{\tt solve($\{$3*x - 2*y = 1, 4*x + 5*y = 7$\}$,$\{$x,y$\}$);}
\medskip
\itemitem{3)}  By similar computations, now compute polynomials
$p_3(x), p_4(x), p_5(x)$ (in that order) so that $\{p_0,p_1,p_2,p_3,p_4,p_5\}$
is orthogonal, $p_i(1) = 1$, and $\deg p_i = i$ for all $i$.
\itemitem{4)}  The polynomials you are computing here are called
the {\it Legendre polynomials\/}.  Maple has a package that computes
these by a different method, so you can use this to check your work
from parts 2 and 3:  First execute the command
\medskip
\centerline{\tt with(orthopoly);}
\medskip
\itemitem{}  The $n$th Legendre polynomial $p_n(x)$ is written
as {\tt P(n,x)} in Maple.  List out {\tt P(0,x)} through {\tt P(5,x)}
and verify that your results from 2 and 3 are correct.
\itemitem{5)}  Find the roots of each of $p_0,\ldots, p_5$
using the {\tt complex} option in the {\tt fsolve} command (consult
the Maple online help listings for the syntax of this one).
\itemitem{6)}  Compute the weights $c_j$, $j = 1,\ldots,n$ from
(1) using the roots of each of $p_i$, $i = 2, \ldots 5$.  (There
is one set of such weights for each $i$.)
\itemitem{7)}  Use Gaussian quadrature with $n = 2,3,4,5$
to approximate the integrals
$$\int_{-1}^1 \cos(\sin(x) - 4)\ dx$$
and 
$$\int_{-1}^1 {\rm exp}(-x^2)\ dx$$
and compare with Maple's approximations for these integrals
(output from commands {\tt evalf(Int(f(x),x=-1..1));}).  What
is happening as $n$ increases?
\itemitem{8)}  How would you use Gaussian quadrature to 
approximate $\int_a^b f(x)\ dx$ if $[a,b] \ne [-1,1]$?  Explain
and use your method to approximate
$$\int_1^2 {\cos^2(x) \over x}\ dx$$
\bigskip
\item{D)} (For after the lab)  In this problem, you will prove
that the Legendre polynomial $p_n(x)$ actually has $n$ distinct
real roots in $[-1,1]$ for all $n$.  
\itemitem{1)}  First note that $\int_{-1}^1 p_n(x)\ dx = 0$ (why?).
Use this to deduce that $p_n(x) = 0$ has at least one solution in 
$[-1,1]$.  
\itemitem{2)} Aiming for a contradiction, suppose that $p_n$ changes sign 
only some
number $1 \le j < n$ times in $[-1,1]$, say at $r_1,\ldots, r_j$.
What can you say about 
$$\int_{-1}^1 (r_1 - x)(r_2 - x) \cdots (r_j - x)p_n(x)\ dx?$$
Show that this leads to a contradiction to the orthogonality
of $\{p_0,p_1,\ldots,p_n\}$.
\bigskip
\item{E)} (Also for after the lab) Problem 7 from Section 4.7 in Burden 
and Faires.
\bigskip
\noindent
{\it Assignment}
\bigskip
{\bf Due:}  In class, Monday, November 19.
\bye 
  
    
