Numerical Integration

This section shows how to use MATLAB to calculate values of definite integrals using approximate methods. We show how to use MATLAB to obtain the closed-form solution of some integrals.

Trapezoidal Integration

The simplest way to find the area under a curve is to split the area into rectangles Figure 8.2-1a). If the widths of the rectangles are small enough, the sum of their areas gives the approximate value of the integral. A more sophisticated method is to use trapezoidal elements (Figure 8.2-1b). Each trapezoid is called a panel. It is not necessary to use panels of the same width; to increase the method’s accuracy you can use narrow panels where the function is changing rapidly. When the widths are adjusted according to the function’s behavior, the method is said be adaptive. MATLAB implements trapezoidal integration with the trapz function. Its syntax is trap z (x, y) , where the array y contains the function values at the points contained ‘in the array x. If you want the integral of a single function, then y is a vector. To integrate more than one function, place their values a matrix y; typing trap z (x, y) will compute the integral of each column of y.

You cannot directly specify a function to integrate with the trapz function; you must first compute and store the function’s values ahead·of time in an array. Later we discuss two other integration functions, the quad and functions, that can accept functions directly. However, they cannot handle arrays of values. So the functions complement one another. These functions a,re summarized in Table 8.2-1.

Figure 8.2-1 Illustration of (a) rectangular and (b) trapezoidal numerical integration.

Figure 8.2-1 Illustration of (a) rectangular and (b) trapezoidal numerical integration.

Table 8.2-1 Numerical integration functions

Table 8.2-1 Numerical integration functions

As a simple example of the use of the trapz function, let us compute the integral


The exact answer is


To investigate the effect of panel width, let us first use,10 panels with equal widths .
of π /10. The script file is


The answer is 1.9797, which gives a relative error of 100(2 – 1.9797)/2) = 1%. Now try 100 panels of equal width; replace the array x with x = linspace (0, pi, 1Q 0) . The answer is 1.9998 for a relative error of 100(2 – 1.9998)/2 = 0.01%. If we examine the plot of the integrand sin r, we would see that the function is changing faster near x = 0 and x = n than near x = π /2. Thus we could achieve the same accuracy using fewer panels if narrower panels are used near x = () and x =n .

When numerical integration was done by hand (before World War II), it was important to use as few panels as necessary to achieve the desired accuracy. However, with the speed of modern computers, it is not difficult to find a reasonable number of panels of uniform width to achieve the required accuracy for many problems. The following method usually works: Compute the integral with a reasonable number of panels (say, 100). Then double the number of panels and compare the answers. If they are close, you have the solution. If not, continue increasing the number of panels until the answers converge to a common value. This method does not always give accurate results, but it can be tried before using variable panel widths or methods more sophisticated than “trapezoidal integration.

We can use numerical integration to find the velocity when either the acceleration function a(t) cannot be integrated or the acceleration is given as a table of values. The following example illustrates the latter case.

Quadrature Functions

As we have just seen, when the integrand is a linear function (one having a plot that is a straight line), trapezoidal integration gives the exact answer. However, if the integrand is not a linear function, then the trapezoidal representation will be inexact. We can represent the function’s curve by quadratic functions to obtain more accuracy. This approach is taken with Simpson’s rule, which divides the integration range b – a into an even number of sections and uses a different quadratic for each pair of adjacent panels. A quadratic function has three parameters, and Simpson’s rule computes these parameters by requiring the quadratic to pass through the function’s three points corresponding to the two adjacent panels. To obtain more accuracy, we can use polynomials of a degree higher than two.

MATLAB function quad implements an adaptive version of Simpson’s rule, while the quadl function is based on an adaptive Lobatto integration algorithm. The term quad is an abbreviation of quadrature, which is an old term for the process of measuring areas. The syntax of both functions is identical and is suroma-rized in Table8.2-1. In its basic form, the syntax is quad ( ‘function’ ,a, b), where ‘ function’ is the name of the function representing the integrand, a is the lower integration limit, and b is the upper limit. To illustrate, let us compute an integral we already know.


The session consists of one command:

»A = quad ( ‘sin’ , 0, pi)

The answer given by MATLAB is 2, which is correct. We use quad1 the same way; namely, A = quad1 ( ‘sin’ ,0, pi).

Although the quad and quad1 functions are more accurate than trapz, they are restricted to computing the integrals of functions and cannot be used when the integrand is specified by a set of points. Thus we could not have used quad and quad1 in Example 8.2-1.

Slope Function Singularities

In addition to singularities of the integrand, another condition can cause problems for numerical integration methods. This occurs when the slope of the integrand becomes infinite either on or within the integration limits; that is, the slope function has a singularity. A simple example of this is the square root function’ y = √X, shown in the top graph of Figure 8.2-2. The slope of this function is its derivative, which is dy / dx = O.5/ √X. This slope is plotted in the bottom graph of the figure.
Note that it becomes infinite at x = 0.

Figure 8.2-2 A function having a singularity in its slope function. The top graph shows the function y = √x. The bottom graph shows the derivative of y = √x The slope has a singularity at x = 0.

Figure 8.2-2 A function having a singularity in its slope function. The top graph shows the function y = √x. The bottom graph shows the derivative of y = √x The slope has a singularity at x = 0.

Let us use all three integration functions and compare their performance in computing the following integral:


To use the trapz function, we must evaluate the integrand at a chosen number of points. The number of points determines the panel width. Here we will use a spacing of 0.01, which requires 100 panels. The script file follows.


The answers are A1 = A2 = 0.6667and A3 = 0.6665:The correct answer, which can be obtained analytically, is 2/3 = 0.6667 to four decimal places. The quad and quad1 functions can produce a message warning that a singularity might be present; however, they do not produce such a message here. This indicates that these functions are capable of dealing with this singularity. The trapz is not affected by the slope singularity here, and it gives a reasonably accurate answer.

You can use quad and quad1 to integrate user-defined functions, as shown in the following example. Because the quad and quad1 functions call the integrand function using vector arguments, you must always use array operations when defining the function. The following example shows how this is done.

The quad functions have some optional arguments. For example. one syntax of the quad function is

quad (‘function’,a,b, tol )

where tol indicates the specified error tolerance. The function iterates until the relative error is less than tol. The default value of tol is 15*eps.

Test Your Understanding
T8.2-2Use both the quad and quadl functions to compute the integral
151 -dx . 2 X

and compare the answers with that obtained from the closed-form solution.
T8.2-3Use a tolerance of 0.001 to integrate the square root function from 0 to I by typing quad ( r sqrt I I 0I 1, 0 . 001) .Do the same using the quadl function, and compare with the results obtained with the default tolerance

Share This