Page 180 - Numerical methods for chemical engineering

166 4 Initial value problems

N, and for all p(x) ∈ 2N−1 ,

' b N

w(x)p(x)dx = w j p(x j ) (4.72)
a j=1
That is, if we place the support points at the zeros of the orthogonal polynomial p N (x),we
obtain the exact value of the definite integral for any polynomial integrand of degree 2N –
1 or less. This is much better than we can usually expect, since with only N support points,
we generally can represent exactly only polynomials of degreeN–1or less.

A proof of this theorem is provided in the supplemental material.

Gaussian quadrature with w (x) = 1 and the use of quadl
- b
To compute the integral f (x)dx, with w(x) = 1, we define ξ such that
a
a + b b − a

x(ξ) = + ξ x(−1) = a x(1) = b (4.73)
2 2
and thus
' b b − a ' 1
f (x)dx = f [x(ξ)]dx (4.74)
2
a −1
We use the set of orthogonal polynomials for [a, b] = [1, 1], w(x) = 1, the Legendre poly-
nomials. We select an integration order through our choice of N, and compute the roots
of p N (x), e.g. by the eigenvalue method. Then, we solve for the weights, and evaluate the
integral. The node locations and quadrature weights can be stored for repeated use.
quadl applies Gaussian quadrature adaptively. Actually, quadl uses as a default Lobatto
quadrature, which adds two more support points at a and b to the roots of p N (x) within (a,
b). However, sometimes we want all support points to lie in the interior of the domain. Let
us say that we wish to integrate a function with a singularity at a < x s < b; i.e., f (x s ) blows
up to ±∞.
- b
Let us say that this singularity is integrable so that f (x)dx exists and is finite. Then,
a
we compute the integral by splitting (a, b)as
' b ' x s ' b
f (x)dx = f (x)dx + f (x)dx (4.75)
a a x s
and applying Gaussian quadrature to each integral to avoid evaluating f (x s ).
We demonstrate quadl to integrate f (x) = e −x over [0, 1]. Continuing the calculations
that demonstrated trapz and quad,
int quadl = quadl(‘calc integrand nexp’,0,1),
int quadl = 0.6321
Like quad, quadl expects the integrand routine to return a vector of function values for an
input vector of argument values.