Page 177 - Numerical methods for chemical engineering
P. 177
Gaussian quadrature 163
Use of trapz and quad
We demonstrate the use of these quadrature functions by computing the integral of f (x) =
e −x over [0, 1], for which
' 1 1 0
I F = e −x dx =−e −x = e −x = 1 −e −1 = 0.6321 (4.57)
0 0 1
Using trapz, the integral and the corresponding error are computed by
N = 50;
x = linspace(0,1,N);
f = exp(-x);
int trapz = trapz(x,f),
int trapz = 0.6321
int exact = 1 – exp(-1),
int exact = 0.6321
error trapz = int exact – int trapz,
error trapz = -2.1939e-005
The accuracy of the calculation is improved by increasing the number of support points N.
This integral can also be computed using quad. First, we write a routine that returns the
integrand function,
function f = calc integrand nexp(x);
f = exp(-x);
return;
and then compute the integral using quad,
int quad = quad(‘calc integrand nexp’,0,1),
int quad = 0.6321
error quad = int exact − int quad,
error quad = -1.3768e-009
quad expects the integrand routine to return a vector of function values for an input vector
of argument values. The accuracy can be increased by adjusting the tolerance through an
additional argument,
int quad 2 = quad(‘calc integrand nexp’,0,1,1e-14),
int quad 2 = 0.6321
error quad 2 = int exact – int quad 2,
error quad 2=0
Gaussian quadrature
In Newton–Cotes integration, the support points are equally spaced; however, we might ask
whether there is a more accurate choice of support points. We examine this issue for the