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
The area under this method for the interval -h≤x≤h is
The curve passes through the three points(-h,y0 ),(0,y1 ),and (h,y2 ). Then, by an equation, we have:
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
αx2-βh=y0-y1…..equation4
αx2+βh=y2-y1…..equation5
Adding of equation4 with equation 5 yields
2αh2=y0-2y1+y2……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 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
Z = trapz(Y)
Z = trapz (X, Y)
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 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
x5=linspace (1,2,5);
x10=linspace (1,2,10);
y5=1. /x5; y10=1. /x10;
area5=trapz (x5, y5), area10=trapz (x10, y10)
MATLAB display the following result:
area5 =
0.6970
area10 =
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
t=linspace (0,2,10);
y=exp(-t.^2);
area=trapz (t, y)
MATLAB displays the following result:
area =
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(e0.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
t=sym(‘t’);
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.
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.
t=linspace (0,10,100);
v=sin (3. *t); i=0.1. *(exp (0.2.*v)-1); p=v.*i;
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.
energy=trapz (t, p), plot (t, energy, ‘+’); grid; title (‘Energy vs Time’); …
xlabel(‘seconds’); ylabel(‘joules’)
energy =
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 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:
0 P1 x1that is equal to plus the area of the trapezoid x1 P1 P2 x2 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
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).
function y = erfcousin(x);
% ERFCOUSIN function to evaluate exp(-x^2)
y= exp(-x.^2); % the array operator .^ is used for vector x.
Step2: Let us use quad with its simplest syntax:
>> y=quad (‘erfcousin’, 1/2, 3/2) % here a=1/2, b=3/2
y=
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.