Writing Fast Matlab Code #9:Numerical Integration
Quadrature formulas are numerical approximation of integrals of the form
where the xk are called the nodes or abscissas and wk are the associated weights. Simpson’s rule is
Simpson’s rule is a quadrature formula with nodes a, a + h, b and node weights h/3,4/3h, h/3 .
Matlab has several functions for quadrature in one dimension:
| quad | adaptive Simpson | Better for low accuracy and nonsmooth integrands |
| quadl | adaptive Gauss-Lobatto | Higher accuracy with smoother integrands |
| quadgk | adaptive Gauss-Kronrod | Oscillatory integrands and high accuracy |
| quadv | adaptive Simpson | Vectorized for multi-valued integrands |
The quad- functions are robust and precise, however, they are not very efficient. They use an adaptive refinement procedure, and while this reduces the number of function calls, they gain little from vectorization and incur significant overhead.
If an application requires approximating an integral whose integrand can be efficiently vectorized, using nonadaptive quadrature may improve speed.
for n = −20:20 % Compute Fourier series coefficients c(n + 21) = quad(@(x)(exp(sin(x).ˆ6).*exp(−1i*x*n)),0,pi,1e−4); end
This code runs in 5.16 seconds. In place of quad, using Simpson’s composite rule with N = 199 nodes yields results with comparable accuracy and allows for vectorized computation. Since the integrals are all over the same interval, the nodes and weights only need to be constructed once.
N = 199; h = pi/(N−1); x = (0:h:pi).'; % Nodes w = ones(1,N); w(2:2:N−1) = 4; w(3:2:N−2) = 2; w = w*h/3; % Weights for n = −20:20 c(n + 21) = w * ( exp(sin(x).ˆ6).*exp(−1i*x*n) ); end
This version of the code runs in 0.02 seconds (200 times faster). The quadrature is performed by the dot product multiplication with w. It can be further optimized by replacing the for loop with one vector-matrix multiply:
[n,x] = meshgrid(−20:20, 0:h:pi); c = w * ( exp(sin(x).ˆ6).*exp(−1i*x.*n) );
For this example, quadv can be used on a multi-valued integrand with similar accuracy and speed,
n = −20:20;
c = quadv(@(x)exp(sin(x).ˆ6).*exp(−1i.*x.*n),0,pi,1e−4);
9.1 One-Dimensional Integration
|
|
is approximated by composite Simpson’s rule with
h = (b − a)/(N−1); x = (a:h:b).'; w = ones(1,N); w(2:2:N−1) = 4; w(3:2:N−2) = 2; w = w*h/3; I = w * f(x); % Approximately evaluate the integral
where N is an odd integer specifying the number of nodes.
A good higher-order choice is 4th-order Gauss-Lobatto [4] (as used by quadl), based on
N = max(3*round((N−1)/3),3) + 1; % Adjust N to the closest valid choice h = (b − a)/(N−1); d = (3/sqrt(5) − 1)*h/2; x = (a:h:b).'; x(2:3:N−2) = x(2:3:N−2) − d; x(3:3:N−1) = x(3:3:N−1) + d; w = ones(1,N); w(4:3:N−3) = 2; w([2:3:N−2,3:3:N−1]) = 5; w = w*h/4; I = w * f(x); % Approximately evaluate the integral
The number of nodes N must be such that (N − 1)/3 is an integer. If not, the first line adjusts N to the closest valid choice. It is usually more accurate than Simpson’s rule when f has six continuous derivatives, f ∈ C 6 (a, b).
A disadvantage of this nonadaptive approach is that the accuracy of the result is only indirectly con- trolled by the parameter N. To guarantee a desired accuracy, either use a generously large value for N or, if possible, determine the error bounds 
[6,7].
where h = (b-a)/(N-1). Note that these bounds are valid only when the integrand is sufficiently differentiable: f N −1 , must have four continuous derivatives for the Simpson’s rule error bound, and six continuous derivatives for 4th-order Gauss-Lobatto.
Composite Simpson’s rule is a fast and effective default choice. But depending on the situation, other methods may be faster and more accurate:
- If the integrand is expensive to evaluate, or if high accuracy is required and the error bounds above are too difficult to compute, use one of Matlab’s adaptive methods. Otherwise, consider composite methods.
- Use higher-order methods (like Gauss-Lobatto/quadl) for very smooth integrands and lower-order methods for less smooth integrands.
- Use the substitution u = 1/(1-x) or Gauss-Laguerre quadrature for infinite-domain integrals like f ∞.
9.2 Multidimensional Integration
An approach for evaluating double integrals of the form 
%%% Construct Simpson nodes and weights over x %%% h = (b − a)/(Nx−1); x = (a:h:b).'; wx = ones(1,Nx); wx(2:2:Nx−1) = 4; wx(3:2:Nx−2) = 2; wx = w*h/3; %%% Construct Simpson nodes and weights over y %%% h = (d − c)/(Ny−1); y = (c:h:d).'; wy = ones(1,Ny); wy(2:2:Ny−1) = 4; wy(3:2:Ny−2) = 2; wy = w*h/3; %%% Combine for two−dimensional integration %%% [x,y] = meshgrid(x,y); x = x(:); y = y(:); w = wy.'*wx; w = w(:).'; I = w * f(x,y); % Approximately evaluate the integral
Similarly for three-dimensional integrals, the weights are combined with
[x,y,z] = meshgrid(x,y,z); x = x(:); y = y(:); z = z(:); w = wy.'*wx; w = w(:)*wz; w = w(:).';
|
|
When the integration region is complicated or of high dimension, Monte Carlo integration techniques are appropriate. The disadvantage is that an N -point Monte Carlo quadrature has error on the order O(1/sqrt(n)), so many points are necessary even for moderate accuracy. Suppose that N points,x1 , x2 , . . . , xN , are uniformly randomly selected in a multidimensional region (or volume) Ω. Then
To integrate a complicated region W that is difficult to sample uniformly, find an easier region Ω that contains W and can be sampled [5]. Then
χW (x) is the indicator function of W : χW (x) = 1 when x is within region W and χW (x) = 0 otherwise. Multiplying the integrand by χW sets contributions from outside of W to zero.
%%% Uniformly randomly sample points (x,y) in Ω %%% x = 4*rand(N,1) − 2; y = 4*rand(N,1) − 2; %%% Restrict the points to region W %%% i = (cos(2*sqrt(x.ˆ2 + y.ˆ2)).*x <= y & x.ˆ2 + y.ˆ2 <= 4); x = x(i); y = y(i); %%% Approximately evaluate the integrals %%% area = 4*4; % The area of rectangle Ω M = (area/N) * length(x); Mx = (area/N) * sum(x); My = (area/N) * sum(y);
Region W sampled with N = 1500. The center of mass is ≈ (0.5, 0.7).
More generally, if W is a two-dimensional region contained in the rectangle defined by a ≤ x ≤ b and
c ≤ y ≤ d, the following code approximates 
x = a + (b−a)*rand(N,1); y = c + (d−c)*rand(N,1); i = logical(indicatorW(x,y)); x = x(i); y = y(i); area = (b−a)*(d−c); I = (area/N) * sum(f(x,y)); % Approximately evaluate the integral
where indicatorW(x,y) is the indicator function χW (x, y) for region W .
For refinements and variations of Monte Carlo integration, see for example [1].
Popularity: 1% [?]
Translator
Follow Us
Video Posts
