2. Numerical Integration

MATLAB Simpson’s Rule

The trapezoidal and Simpson’s rules are special cases of the Newton-Cote rules which use higher degree functions for numerical integration.

Let the parabola represent the curve of the figure.

y=αx2+βx+γ…………..equation 1

MATLAB Simpson's Rule

The area under this method for the interval -h≤x≤h is

MATLAB Simpson's Rule

The curve passes through the three points(-h,y0 ),(0,y1 ),and (h,y2 ). Then, by an equation, we have:

MATLAB Simpson's Rule

We can now evaluate the coefficients α,β,γ and express equation2 in terms of h, y0, y1and y2.This is done with the following procedure.

By substitution of (b) of equation3 into (a) and (c) and rearranging we obtain



Adding of equation4 with equation 5 yields


and by substitution into equation2, we obtain

MATLAB Simpson's Rule

Now, we can apply equation8 to successive segments of any curve y=f(x) in the interval a≤x≤b as shown on the curve of the figure.

MATLAB Simpson's Rule

We observe that a parabola can approximate each segment of width 2h of the curve through its ends and its midpoint. Thus, the area under segment AB is

MATLAB Simpson's Rule

Likewise, the area under segment BC is

MATLAB Simpson's Rule

and so on. When the areas under each segment are added, we obtain

MATLAB Simpson's Rule

Since each segment has width 2h, to apply Simpson’s rule of numerical integration, the number n of subdivisions must be even. This restriction does not apply to the trapezoidal rule of numerical integration.

The value of for equation11 is found from

MATLAB Simpson's Rule
2. Numerical Integration

Matlab Trapz

The MATLAB function trapz (x, y, n) where y is the integral for x, approximates the integral of a function y=f(x) using the trapezoidal rule, and n (optional) performs integration along with dimension n.


  1. Z = trapz(Y)  
  2. Z = trapz (X, Y)  
  3. Z = trapz (…, dim)  


Use the MATLAB function trapz (x, y) to approximate the cost of the integral


and by comparison with the exact value, evaluate the percent error when n=5 and n=10.


The exact value is found from


For the approximation using the trapezoidal rule, we let x5represent the row vector with n=5, and x10 the vector with n=10, that is, ∆x =1/5 and ∆x=1/10, respectively. The corresponding values are denoted as y5and y10, and the areas under the curve as area5 and area10, respectively.

Create the following Script

  1. x5=linspace (1,2,5);   
  2. x10=linspace (1,2,10);   
  3. y5=1. /x5; y10=1. /x10;   
  4. area5=trapz (x5, y5), area10=trapz (x10, y10)  

MATLAB display the following result:

  1. area5 =  
  2.           0.6970  
  3. area10 =  
  4.              0.6939  

The percent error when ∆x =1/5 is used is


The percent error when ∆x =1/10 is used is



The integral


where τ is a dummy variable of integration, is called the error function, and it is used extensively in communications theory. Use the MATLAB trapz (x, y) function to find the area under this integral with n=10 when the upper limit of integration is t=2.


Create the following Script

  1. t=linspace (0,2,10);   
  2. y=exp(-t.^2);   
  3. area=trapz (t, y)     

MATLAB displays the following result:

  1. area =  
  2.             0.8818  


The i-v (current-voltage) relation of a non-linear electrical machine is given by


where v(t)=sin3t.

By any means, find

  • The instantaneous power is
    p(t)=v(t)i(t)=0.1 sin3t(e0.2sin3t-1)
  • The energy is the integral of the instantaneous is
    MATLAB Trapz

An analytical solution of the last integral is possible using integration by parts, but it is not easy. We can try the MATLAB int (f, a, b) function where f is a symbolic expression, and a and b are the lower and upper limits of integration, respectively.

When MATLAB cannot found a solution, it returns a warning. For example, MATLAB returns the following message when integration is attempted with the symbolic expression of equation

  1. t=sym(‘t’);  
  2. s=int (0.1*sin(3*t) *(exp (0.2*sin(3*t))-1),0,10)  

When this script is executed, MATLAB show the following message.

  1. Warning: Explicit integral could not be found.   

Next, we will find and sketch the power and energy by the trapezoidal rule using the MATLAB trapz (x, y) function. For this example, we choose n=100, so that ∆x=1/100. The MATLAB script below will compute and plot the power.

  1. t=linspace (0,10,100);  
  2. v=sin (3. *t); i=0.1. *(exp (0.2.*v)-1); p=v.*i;  
  3. plot(t,p); grid; title(‘Power vs Time’); xlabel(‘seconds’); ylabel(‘watts’)  

The power varies in a uniform fashion as shown by the plot of figure


The MATLAB script below computes and plots the energy.

  1. energy=trapz (t, p), plot (t, energy, ‘+’); grid; title (‘Energy vs Time’); …  
  2. xlabel(‘seconds’); ylabel(‘joules’)  
  3. energy =  
  4.                0.1013  

Thus, the value of the energy is 0.1013 joule. The energy is shown in figure:

2. Numerical Integration

MATLAB Trapezoidal Rule

Consider the function y=f(x) for the interval a≤x≤b, shown in figure:

MATLAB Trapezoidal Rule

To evaluate the definite integral, MATLAB Trapezoidal Ruledx, we divide the interval a≤x≤b into subintervals each of lengthMATLAB Trapezoidal Rule. Then, the number of points between x0=a,and xn=b is x1=a+∆x,x2=a+2∆x,…xn-1=a+(n-1)∆x. Therefore, the integral from a to b is the sum of the integrals from a to x1, from x1 to x2 and so on, and finally from xn-1 to b.

The total area is:

MATLAB Trapezoidal Rule

 0 P1 x1that is equal toMATLAB Trapezoidal Rule plus the area of the trapezoid x1 P1 P2 x2 that is equal to MATLAB Trapezoidal Rule, and so on. Then, the trapezoidal approximation becomesMATLAB Trapezoidal Rule


Using the trapezoidal rule with n=4, estimate the cost of the definite integralMATLAB Trapezoidal Rule

Compare with the exact value and evaluate the percent error.


The exact value of this integral isMATLAB Trapezoidal Rule

For the trapezoidal rule approximation, we haveMATLAB Trapezoidal Rule

and by substitution into an equationMATLAB Trapezoidal Rule

From equation 3 and equation 4, we find that the percent error isMATLAB Trapezoidal Rule

2. Numerical Integration

MATLAB Double Integral

To evaluate integrals of the form

MATLAB Double Integral

MATLAB provides a function dblquad. The calling syntax for dblquad is

I=dblquad (‘fxy-fun’,xmin,xmax,ymin,ymax,tol,@method)

where tol and method are optional input arguments. The optional argument tol specified tolerance (default value is 10-6), as previously discussed for 1-D integration, and method determines a choice that the user makes about the purpose of integration to be used, e.g., quad and quad1. The default method is quad. The user-defined integrand function fxy-fun,must be written that it can accept a vector x and a scalar y while evaluating the integrand.


Let us compute the following integral

MATLAB Double Integral

It is merely to verify analytically that I=4. Let us see how dblquad execute on this integral.

>> F = Inline ('1-6*x.^2*y^' );                      // Create the integrand as Inline function.  
>> I = dblquad (F, 0, 2, -1, 1)  
I =  
2. Numerical Integration

Numerical Integration (Quadrature)

A numerical evaluation of the integral ∫ f(x)dx is known as Quadrature.

The general form of numerical integration of a function f (x) over some interval [a, b] is a weighted total of the function values at a finite number (N + 1) of sample points (nodes), indicated to as ‘quadrature’:

MATLAB Numerical Integration

MATLAB provides built-in functions for numerical integration:


It integrates a specified function over specified limits, build on adaptive Simpson’s rule. The adaptive rule follows to improve accuracy by adaptively selecting the estimate of the subintervals within the limits of integrations while evaluating the sums that make up the integral.


It integrates a specified function over specified limits, build on adaptive Lobatto quadrature. This one is more accurate than a quad, but it also uses more function evaluations.

The syntax for both quad and quad1 are as follows:

MATLAB Numerical Integration

To use quad1, you replace quad with quad1 in the syntax. As shown in syntax, both functions require us to supply the integrand as a userwritten function.

The optional input argument tol specifies absolute tolerance (default value is 10-6). A nonzero value of the other optional argument, trace, shows some intermediate calculations at each step.

The optional arguments p1, p2, etc., are passed on to the user-defined function as input arguments in addition to x.

The steps include in numerical integration using these built-in functions are:

Step1: Write a function that returns the value of the integrand f(x) given the value of x. Our function should be able to accept the input value x as a vector and produce the value of the integrand (the output) as a vector.

Step2: It decides which function to use -quad or quad1(quad is faster but less accurate than quad1). The default tolerance value is 10-6 Find the Integral.


Let us compute the following Integral

MATLAB Numerical Integration

This integration is closely related to the error function, erf.

MATLAB Numerical Integration

MATLAB also provides the error function, erf, as a built-in function, we can evaluate our integral in closed form (in terms of the error function) and compare the results from numerical integration. Let us follow the steps outlined previously.

Step1: Here is the function that evaluates the integrand at given x (x is allowed to be a vector).

  1. function y = erfcousin(x);  
  2. % ERFCOUSIN function to evaluate exp(-x^2)  
  3. y= exp(-x.^2);      % the array operator .^ is used for vector x.  

Step2: Let us use quad with its simplest syntax:

  1. >> y=quad (‘erfcousin’, 1/2, 3/2)         % here a=1/2, b=3/2  
  2. y=  
  3.     0.3949  

The exact result for the integral, up to 10 decimal places, is 0.3949073872. In the preceding example, we have used the default tolerance. Here is a table showing the results of a few experiments with the integral. We use both quad and quad1 for integration. For a quad, we tabulate results with different tolerances. In the table, we list the value of the integral, % error, and the number of function evaluations. For quad1 gives quite an accurate solution with just the default tolerance.

FunctiontolAnswer% ErrorF-evals
quaddefault0.39490739271.3907 x 10-617
10-70.39490738945.5506 x 10-725
10-80.39490738732.4369 x 10-833
quad1default0.39490738758.0616 x 10-818