Mathematics 172 – Numerical Analysis 2 Lab/Problem Set 1 – Advanced Numerical Integration Methods January 26, 1999

Goals

In this lab we will study the adaptive quadrature and Gaussian quadrature methods introduced in class last week.

Gaussian Quadrature

Recall from class last Thursday that the Gaussian quadrature method on an interval [a, b] depends on knowing an orthogonal basis for the vector space of polynomial functions on [a, b] of degree ≤ n for some n – that is, a collection of polynomials

{p_{0}(x), p_{1}(x), …, p_{n}(x)}

with the properties that

a) (orthogonality)

⟨p_{i}, p_{j}⟩ = ∫_{a}^{b} p_{i}(x)p_{j}(x) dx = \{ \begin{matrix}0 & if i ≠ j, and \\ c > 0 & if i = j \end{matrix}

b) deg p_{i}(x) = i for all i = 0, …, n.

c) The equation p_{i}(x) = 0 has i distinct real roots in the interval [a, b] for all i = 0, …, n.

Knowing such a collection of polynomials, and writing x_{1}, …, x_{n} for the roots of the equation p_{n}(x) = 0 on [a, b], we showed that the numerical integration rule:

(1)  ∫_{a}^{b} f(x) dx = ∑_{j=1}^{n} c_{j}f(x_{j})

(where

(2)  c_{j} = ∫_{a}^{b} \frac{∏_{k≠j}(x − x_{k})}{∏_{k≠j}(x_{j} − x_{k})} dx

for each j = 1, … n) has degree of precision 2n − 1(!) A) In this problem, you will use Maple to construct a collection of such polynomials on

[a, b] = [−1, 1] “from scratch.”

1) Check that p_{0}(x) = 1 and p_{1}(x) = x are orthogonal under the inner product

⟨f, g⟩ = ∫_{−1}^{1} f(x)g(x) dx.

Note that p_{0}(1) = p_{1}(1) = 1 also.

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 ∫_{−1}^{1} p_{0}(x)p_{2}(x) dx and ∫_{−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.

                     solve(eqns,vars);

where eqns is the set of equations to be solved, and vars is the set of variables you want to solve for. For example, to solve

\begin{aligned}3x − 2y & = 1 \\ 4x + 5y & = 7 \end{aligned}

you can use

        solve({3*x - 2*y = 1, 4*x + 5*y = 7},{x,y});

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.

4) The polynomials you are computing here are called the 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

                      with(orthopoly);

The nth Legendre polynomial p_{n}(x) is written as P(n,x) in Maple. List out P(0,x) through P(5,x) and verify that your results from 2 and 3 are correct.

5) Find the roots of each of p_{0}, …, p_{5} using the complex option in the fsolve command (consult the Maple online help listings for the syntax of this one).

6) Compute the weights c_{j}, j = 1, …, n from (1) using the roots of each of p_{i}, i = 2, … 5. (There is one set of such weights for each i.)

7) Use Gaussian quadrature with n = 2, 3, 4, 5 to approximate the integrals

∫_{−1}^{1} cos(sin(x) − 4) dx

and

∫_{−1}^{1} exp(−x^{2}) dx

and compare with Maple’s approximations for these integrals (output from commands evalf(int(f(x),x=-1..1));). What is happening as n increases?

8) How would you use Gaussian quadrature to approximate ∫_{a}^{b} f(x) dx if [a, b] ≠ [−1, 1]? Explain and use your method to approximate

∫_{1}^{2} \frac{cos^{2}(x)}{x} dx

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.

1) First note that ∫_{−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].

2) Aiming for a contradiction, suppose that p_{n} changes sign only some number 1 ≤ j < n times in [−1, 1], say at x_{1}, …, x_{j}. What can you say about

∫_{−1}^{1}(x_{1} − x)(x_{2} − x) ⋯ (x_{j} − x)p_{n}(x) dx?

Show that this leads to a contradiction to the orthogonality of {p_{0}, p_{1}, …, p_{n}}.

Adaptive Quadrature

C) Implement the recursive adaptive quadrature procedure 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:

∫_{2}^{5} ln(x^{4} − 2x^{2} + 1) dx, ∫_{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.

Due: In class, Tuesday February 2.