Numerical Analysis

Math 370  MWF 2:30-3:25
Spring 1997  127 Fowler

(04.09.97) \underline{Composite Numerical Integration; Gaussian Quadrature} Error Term of Composite Trapezoid Starting with the formula for the Trapezoid Rule, we can show that the error formula for the Composite Trapezoid Rule is O(Δx^{2}) and has precision 1.

CT = T(N, f, a, b) = \frac{Δx}{2} [f(a) + 2 ∑_{j=1}^{N−1} f(a + jΔx) + f(b)] − \frac{(b − a)}{12}(Δx)^{2}f′′(μ)

where μ ∈ (a, b) and Δx = (b − a)/N, N being number of subintervals Composite Midpoint Method

CM =M(N, f, a, b) = ∑_{j=1}^{N} f[a + (j − 1/2)Δx]Δx + \frac{(b − a)}{24}(Δx)^{2}f′′(μ)

Composite Simpson’s Method

CS = \frac{2}{3}CM + \frac{1}{3}CT − \frac{(b − a)}{2880}(Δx)^{4}f^{(4)}(μ)

\underline{Example} Now that we know expressions for the error terms of these methods we can use this infor-

mation for error control. Consider ∫_{1}^{3} e^{x^{2}} dx. How many N subintervals would we have to split [1,3] into in order to be confident that you have computed the integral to with 10^{−4} of its exact value? Compare how many subsintervals you need using Midpoint versus using Simpson’s.

Effect of Round-Off Error on Composite Numerical Integration Clearly when we do numerical integration we do quite a bit of function evaluation. As we saw in numerical differentiation, everytime we evaluate a function f(x) the computer is using f̃(x) + e(x). Lets consider the Trapezoid Method and try and compute E(Δx) and see what happens as Δx → 0

Gaussian Quadrature Consider the following pictures of f(x) = 1 − 1/x^{2}. Using two function evaluations how would you estimate the area under the curve from 1 to 3?

Karl Friedrich Gauss had the idea of making the x_{1} and x_{2} variables as well as the coefficients c_{1} and c_{2}. Then since this is 4 unknown parameters it should be able to approximate a polynomial which uses 4 unknown parameters to be completely described. What degree polynomial is that? Thus the N = 2 Gaussian Quadrature rule is designed to have precision .

Find x_{1}, x_{2}, c_{1} and c_{2} such that ∫_{−1}^{1} f(x) dx ≈ c_{1}f(x_{1}) + c_{2}f(x_{2}) Let’s write down the four equations for the four unknowns.

The x_{i} are known as the SAMPLE POINTS and the c_{i} are known as the WEIGHTS. Tables of these constants appear in many books. (See page 209 of Burden & Faires)

N SAMPLE POINTS x_{i} WEIGHTS c_{i}
2 1/\sqrt{3} 1
−1/\sqrt{3}  1
3 \sqrt{3/5} 5/9
0  8/9
−\sqrt{3/5}  5/9

\underline{Example} Use 2-point Gaussian Quadrature to approximate ∫_{−1}^{1} 2x^{2}dx What is the exact value? How accurate do we expect our answer to be?

Use 2-point Gaussian Quadrature to approximate ∫_{1}^{3} 1 − \frac{1}{x^{2}} dx

So, in general, how can we transform the integral on [a,b] to an integral on [-1,1]?

∫_{a}^{b} f(x)dx = ∫_{−1}^{1} N-point Gaussian Quadrature

In general, the N-point Gaussian Quadrature rule will have precision 2N − 1. The formula will look like

∫_{−1}^{1} f(x) dx ≈ ∑_{i=1}^{N} c_{i}f(x_{i}) + \frac{2^{2N+1}(N!)^{4}}{(2N + 1)[(2N)!]^{3}} f^{2N}(ξ)

Weights and Sample Points It is clear that the x_{1} and x_{2} are the roots of the polynomial x^{2} − 1/3. This is better known as the N = 2 Legendre Polynomial. In fact, in general the x_{i} for the N-point Gaussian Quadrature Rule will be the roots of the N-th degree Legendre Polynomials. The w_{i} are given by

w_{i} = ∫_{−1}^{1} ∏_{j≠i^{j=1}} \frac{x − x_{j}}{x_{i} − x_{j}} , where x_{i} are roots of P_{n}(x)

Legendre Polynomials The Legendre Polynomials P_{n} are a set of orthonormal polynomials on [-1,1]. This means

that P_{N}(x) is a polynomial of degree N and that ∫_{−1}^{1} P_{k}P_{N}dx = 0 for all k less than N. The Polynomials are determined recursively, using:

P_{0}(x) = 1
P_{1}(x) = x − α_{1}
P_{m+1}(x) = xP_{m} − α_{m+1}P_{m}(x) − β_{m+1}P_{m−1}(x)

where α_{m+1} = \frac{∫ xP_{m}^{2}dx}{∫ P_{m}^{2}dx} and β_{m+1} = \frac{∫ xP_{m}P_{m−1}dx}{∫ P_{m−1}^{2}dx}

The first few Legendre Polynomials are

P_{0}(x) = 1
P_{1}(x) = x
P_{2}(x) = x^{2} − 1/3
P_{3}(x) = x^{3} − 3x/5