\magnification=\magstep1
\centerline{Mathematics 172 -- Numerical Analysis 2}
\centerline{Lab/Problem Set 1 -- Advanced Numerical Integration Methods}
\centerline{January 26, 1999}
\bigskip
\noindent
{\bf Goals}
\bigskip
In this lab we will study the {\it adaptive quadrature\/}
and {\it Gaussian quadrature\/} methods introduced in 
class last week.
\bigskip
\noindent
{\it Gaussian Quadrature}
\bigskip
Recall from class last Thursday 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$(!)

\item{A)}  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{B)} (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 $x_1,\ldots, x_j$.
What can you say about 
$$\int_{-1}^1 (x_1 - x)(x_2 - x) \cdots (x_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
\noindent
{\it Adaptive Quadrature}
\bigskip
\item{C)}  Implement the recursive adaptive quadrature procedure
{\tt Adapt} we discussed in class on Tuesday, January 19.  (Enter
the procedure in a file using the text editor, and read it into
your Maple session.)  Use it to approximate the integrals:
$$\int_2^5 \ln(x^4 - 2x^2 + 1)\ dx, \qquad \int_2^4 \sin(x^4)\ dx$$
accurate to $10^{-5}$.  (Note:  you may need to increase the maximum
recursion level to complete one of these calculations.  Can you
see why?)  Explain what you did, and 
compare with Maple's approximate value.
\bigskip
\noindent
{\bf Due:}  In class, Tuesday February 2.
\bye 
  
    
