Category: 2. Numerical Integration

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=αx²+βx+γ…………..equation 1

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

The curve passes through the three points(-h,y₀ ),(0,y₁ ),and (h,y₂ ). Then, by an equation, we have:

We can now evaluate the coefficients α,β,γ and express equation2 in terms of h, y₀, y₁and y₂.This is done with the following procedure.

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

αx²-βh=y₀-y₁…..equation4

αx²+βh=y₂-y₁…..equation5

Adding of equation4 with equation 5 yields

2αh²=y₀-2y₁+y₂……equation6

and by substitution into equation2, we obtain

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.

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

Likewise, the area under segment BC is

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

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

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.

Syntax

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

Example1

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.

Solution

The exact value is found from

For the approximation using the trapezoidal rule, we let x₅represent the row vector with n=5, and x₁₀ the vector with n=10, that is, ∆x =1/5 and ∆x=1/10, respectively. The corresponding values are denoted as y₅and y₁₀, 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

Example2

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.

Solution

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  

Example3

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(e^0.2sin3t-1)
• The energy is the integral of the instantaneous is
  

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:

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

To evaluate the definite integral, dx, we divide the interval a≤x≤b into subintervals each of length. Then, the number of points between x₀=a,and x_n=b is x₁=a+∆x,x₂=a+2∆x,…x_n-1=a+(n-1)∆x. Therefore, the integral from a to b is the sum of the integrals from a to x₁, from x₁ to x₂ and so on, and finally from x_n-1 to b.

The total area is:

 0 P₁ x₁that is equal to plus the area of the trapezoid x₁ P₁ P₂ x₂ that is equal to , and so on. Then, the trapezoidal approximation becomes

Example

Using the trapezoidal rule with n=4, estimate the cost of the definite integral

Compare with the exact value and evaluate the percent error.

Solution:

The exact value of this integral is

For the trapezoidal rule approximation, we have

and by substitution into an equation

From equation 3 and equation 4, we find that the percent error is

To evaluate integrals of the form

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.

Example

Let us compute the following 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 =  
       4.0000  

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 provides built-in functions for numerical integration:

quad

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.

quadl

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:

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.

Example

Let us compute the following Integral

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

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.

Function	tol	Answer	% Error	F-evals
quad	default	0.3949073927	1.3907 x 10-6	17
	10-7	0.3949073894	5.5506 x 10-7	25
	10-8	0.3949073873	2.4369 x 10-8	33
quad1	default	0.3949073875	8.0616 x 10-8	18