CHAPTER 4: MATHEMATICAL MODELING WITH MATLAB
Lecture 4.2: Summation rules
for numerical integration
Trapezoidal, Simpson and midpoint rules for
integrals:
Problem: Given a set of data points:
(x_{0},y_{0});
(x_{1},y_{1}); (x_{2},y_{2}); …(x_{n},y_{n})
Suppose
the data points represent a function y = f(x). Suppose also that
the data points are equally spaced with constant step size h = x_{1 }
x_{0}. Find a numerical approximation for the integral, which
is the signed area under the curve y = f(x) between the end
points (x_{0},y_{0}) and (x_{n},y_{n}).
Example:
Linear
electrical circuits can be easily miniaturized if they do not include large and
bulky inductors. When a current I = I(t) is applied to the input
port of a simple resistorcapacitor oneport network, the voltage V = V(t) develops across the port
terminals. At the time instance t = T, the voltage output is
determined as a sum of the voltage drop across the resistor (which is R
I(T)) and of the voltage drop across the capacitor (which is V_{0}
+ _{}/ C), where V_{0} is an initial
voltage. If the input current I(t) can be measured at different
times t = t_{k} for k = 0,1,2,…,n, such
that t_{0} = 0 and t_{n} = T, them
the integral is to be evaluated numerically from the given data set.
Solution:
The
function y = f(x) is either analytically defined or given in a
tabular form. The numerical integration is based on the use of numerical
interpolation y = P_{n}(x) fitted to the given data
points and analytical integration of the polynomial P_{n}(x). This
is a socalled NewtonCotes integration algorithm. The most important
NewtonCotes integration formulas are trapezoidal, Simpson and midpoint rules.

Trapezoidal rule: A
linear interpolation between the points (x_{0},y_{0})
and (x_{1},y_{1}) approximates the area under
the curve y = f(x) by the area of the trapezoid: _{} _{} I_{trapezoidal}(f;x_{0},x_{1})
= _{}( y_{1 }+ y_{0 }) Trapezoidal
rule is popular in numerical integration as it is a simple method. Although
its accuracy is low, the accuracy can be controled by doubling the number of
elementary subintervals (trapezoids). 

Simpson rule: A
quadratic interpolation between the points (x_{0},y_{0})
(x_{1},y_{1}), and (x_{2},y_{2})
approximates the area under the curve y = f(x) by the
area under the interpolant: _{} _{} I_{Simpson}(f;x_{0},x_{2})
= _{}( y_{0 }+ 4y_{1 }+ y_{2 }) Simpson rule is popular because of high accuracy of
numerical integration compared to the trapezoidal rule. 

Midpoint rule: A
constant interpolation between the point (x_{1},y_{1}),
centered in the interval between (x_{0},y_{0})
and (x_{2},y_{2}), approximates the area under
the curve y = f(x) by the area of a rectangle centered at the
midpoint: _{} _{} I_{midpoint}(f;x_{0},x_{2})
= 2h y_{1} Midpoint
rule is popular in numerical integration of functions with singularities at
the end of the interval. It has the same accuracy as the trapezoidal rule and
is often used in combination with the trapezoidal rule for computations of
integrals near singularities. 
h = 0.1; x0 = 0; x1 = x0+h; x2 = x0+2*h;
% the three data points are
taken on [0,1] with equal step size
y0 = sqrt(1x0^2); y1 = sqrt(1x1^2); y2 = sqrt(1x2^2);
yIexact = quad('sqrt(1x.^2)',x0,x2); % 'exact answer' is computed
by MATLAB
yItrap = h*(y0+y1+y1+y2)/2; % trapezoidal rule for two
subintervals with h
yIsimp = h*(y0+4*y1+y2)/3;
% Simpson rule for two subintervals with h
yImid = 2*h*y1; %
midpoint rule for two subintervals with h
fprintf('Exact = %6.6f\nTrapezoidal = %6.6f\nSimpson = %6.6f\nMidpoint = %6.6f',yIexact,yItrap,yIsimp,yImid);
Exact = 0.198659
Trapezoidal = 0.198489
Simpson = 0.198658
Midpoint = 0.198997
Composite summation rules:
The
summation rules are extended to multiple intervals, when the function y =
f(x) is represented by (n+1) data points with constant
step size h. The composite rule is obtained by summating areas of
all n individual areas.
·
Composite trapezoidal rule:
_{} _{} I_{trapezoidal}(f;x_{0},x_{1},…,x_{n})
= _{}( y_{0 }+ 2 y_{1
}+ 2 y_{2} + … + 2 y_{n1 }+ y_{n} )
·
Composite Simpson rule:
_{} _{} I_{simpson}(f;x_{0},x_{1},…,x_{n})
= _{}( y_{0 }+ 4 y_{1
}+ 2 y_{2} + 4 y_{3 }+ 2 y_{4} + … + 2 y_{n2 }+ 4 y_{n1} +
y_{n} )
·
Composite midpoint rule:
_{} _{} I_{midpoint}(f;x_{0},x_{1},…,x_{n})
= 2h ( y_{1 }+ y_{3 }+
… + y_{n3 }+ y_{n1} )
For
the composite Simpson and midpoint rules, the total interval between x _{}[x_{0},x_{n}] has to be divided into
even number of subintervals.
Errors of numerical integration:
Numerical
integration is much more reliable process compared to numerical
differentiation. Rounding errors in computing sums in numerical integration are
always constant, which are independent neither of the integration rule nor from
the number of subintervals for summation. Truncation errors can be reduced with
the use of more accurate summation rules for numerical integration.
Consider
the equally spaced data points with constant step size: h = x_{2}
– x_{1 }= x_{1} –_{ } x_{0}. The theory based on the Taylor expansion
method shows the following local truncation errors:
·
Trapezoidal rule:
_{}f(x)dx – I_{trapezoidal}(f;x_{0},x_{1 })
= _{}f''(x), x _{} [x_{0},x_{1}]
The
truncation error of the trapezoidal rule is proportional to h^{3},
i.e. it has the order of O(h^{3}). The error is also
proportional to the second derivative of the function f(x) at an
interior point x of the integration interval. Thus, the
trapezoidal rule is exact for a linear function f(x).
·
Simpson rule:
_{}f(x)dx – I_{Simpson}(f;x_{0},x_{2 })
=  _{} f''''(x), x _{} [x_{0},x_{2}]
The
truncation error of the Simpson rule is proportional to h^{5} rather
than h^{3}, i.e. it has the order of O(h^{5}).
The error is also proportional to the fourth derivative of the function f(x)
at an interior point x of the integration interval. The Simpson
rule is exact for polynomial functions f(x) of order m =
0,1,2,3.
·
Midpoint rule:
_{}f(x)dx – I_{midpoint}(f;x_{0},x_{2 })
= _{} f''(x), x _{} [x_{0},x_{2}]
The
truncation error of the midpoint rule is as bad as that of the trapezoidal
rule.
The global truncation error is computed for composite
integration rules when an interval between x_{}[x_{0},x_{n}] is divided into n
partial subintervals. The global truncation error is obtained by
summation of local truncation errors and all local rounding errors:
·
Composite trapezoidal rule:
e_{trapezoidal }=  _{}  I_{trapezoidal}(f;x_{0},x_{1},…,x_{n})
 < M_{2} _{}(x_{n} – x_{0}) + eps (x_{n –} x_{0}),
where M_{2} = max  f''(x). The
global truncation error is proportional to the length of the integration
interval and is order of O(h^{2}).
·
Composite Simpson rule:
e_{Simpson }=  _{}  I_{Simpson}(f;x_{0},x_{1},…,x_{n})
 < M_{4} _{}(x_{n} – x_{0}) + eps (x_{n –} x_{0}),
where M_{4} = max  f''''(x). The
global truncation error is proportional to the length of the integration
interval and is order of O(h^{4}).
h = 0.1; k = 1;
while (h > 0.0000001)
x = 0 : h : 1; y =
sqrt(1.x.^2); n = length(x)1;
yIexact = pi/4; % exact
integral, 1/4 of area of a unit disk
yItrap =
h*(y(1)+2*sum(y(2:n))+y(n+1))/2; % composite trapezoidal rule
yIsimp =
h*(y(1)+4*sum(y(2:2:n))+2*sum(y(3:2:n1))+y(n+1))/3;
yImid =
2*h*sum(y(2:2:n)); % composite
midpoint rule
eItrap(k) = abs(yItrapyIexact); eIsimp(k) = abs(yIsimpyIexact);
eImid(k) = abs(yImidyIexact); hst(k) = h; h = h/2; k = k+1;
end,
plot(hst,eItrap,'g:',hst,eIsimp,'b',hst,eImid,'r');
.
If
the step size h between two points becomes smaller, the
truncation error of the summation rule decreases. It decreases faster for
Simpson rule and slower for trapezoidal and midpoint rule. For example, if h
is halved, the global truncation error of the Simpson rule is
reduced by a factor of 16, while the global truncation errors of
the trapezoidal and midpoint rules are reduced only by a factor of 4.
Since
the rounding error is bounded by the integration interval multiplied by the
mashine precision eps, the error of numerical integration can be
reduced to that constant global number.
Example:
The
numerical approximations of the integral _{}, where I(t) is the current in a
resistorcapacitor network, are obtained with step size h = 10
(green pluses) and with step size h
= 5 (blue dots), versus the exact
integral (red solid curve). The composite trapezoidal rule's error
reduces with smaller step size h (blue dots are closer to the
exact red curve compared to the green pluses). The figure also shows that the
global truncation error for the integral grows with the length of the interval T.
MATLAB
numerical integration:
·
quad: evaluates
numerically the integral of a function by using adaptive Simpson quadrature
·
quadl: evaluates numerically the integral of a function by using adaptive
Lobatto quadrature
·
dblquad: evaluates the double integral of a function of two variables in a
rectangular domain
% The function for numerical integration
should be written as a MATLAB Mfile
% The functon should
accept a vector argument X and return a vector result Y
function [Y] = integrand(X)
% this Mfile sets up a
function y = f(x) = sqrt(1 + exp(x))
% which is the integrand
for the integral to be evaluated by function "quad"
Y = sqrt(1 + exp(X));
%
Compute the integral: I = int_0^2 sqrt(1 + exp(x)) dx
format long; I1 = quad(@integrand,0,2)
% the default tolerance is 10^(6) for absolute error of numerical
integration
tolerance = 10^(8); I2 = quad(@integrand,0,2,tolerance)
tolerance = 10^(12); I3 = quad(@integrand,0,2,tolerance)
I4 = quadl(@integrand,0,2)
h = 0.0001; x = 0 : h : 2; y = feval(@integrand,x); n =
length(y)1;
I5 = h*(y(1)+2*sum(y(2:n))+y(n+1))/2 % composite trapezoidal rule
I6 = h*(y(1)+4*sum(y(2:2:n))+2*sum(y(3:2:n1))+y(n+1))/3 %
composite Simpson rule
I7 = 2*h*sum(y(2:2:n)) % composite midpoint rule
I2 = 4.00699422326706
I3 = 4.00699422325470
I4 = 4.00699422322700
I5 = 4.00699422402304
I6 = 4.00699422325470
I7 = 4.00699422171802
W1 =
quad('sin(pi*x.^2)',0,1), W2 = quad(inline('sin(pi*x.^2)'),0,1)
% standard MATLAB functions can be used for numerical integration as strings
W2 = 0.5049
Romberg integration for higherorder
NewtonCotes integration formulas:
More
accurate integration formulas with smaller truncation error can be obtained by
interpolating several data points with higherorder interpolating polynomials.
For example, the thirdorder interpolating polynomial P_{3}(x) between
four data points leads to the Simpson 3/8 rule:
_{} _{} I_{Simpson 3/8}(f;x_{0},x_{3})
= _{}( y_{0 }+ 3y_{1 }+ 3y_{2 }+ y_{3 })
while
the fourthorder interpolating polynomial P_{4}(x) between
five data points leads to the Booles rule:
_{} _{} I_{Booles}(f;x_{0},x_{4})
= _{}( 7 y_{0 }+ 32 y_{1 }+ 12 y_{2 }+ 32
y_{3 } + 7 y_{4} )
The
higherorder integration formulas can be recovered with the use of the Romberg
integration algorithm based on the Richardson extrapolation algorithm.
·
Recursive integration formulas:
The
composite trapezoidal rule has the global truncation error of order O(h^{2}).
Denote the composite trapezoidal rule for the integral of f(x) between
[x_{0},x_{n}] as^{ }R_{1}(h).
Compute two numerical approximations of the integral with two step
sizes h and 2h:
_{}f(x) dx = R_{1}(h) + _{}h^{2 };^{ }_{}f(x) dx = R_{1}(2h) + 4_{}h^{4};
where
_{} is unknown
coefficient for the global truncation error. The number of trapezoids must be
even in order the numerical approximation with double stepsize (2h)
could be computed. By cancelling the truncation error of order O(h^{2}),
we define a new integration rule for the same integral:
_{}f(x) dx = _{} = R_{2}(h)
The
new integration rule R_{2}(h) for the same integral is
more accurate since the truncation error is of order O(h^{4}). It
is in fact the composite Simpson's rule as it can be checked directly. If the
step size is sufficiently small, the composite Simpson's rule gives a much
better numerical approximation for the integral, compared to the composite
trapezoidal rule.
The
process can be continued to find a higherorder integration formulas R_{m}(h)
with the truncation error O(h^{2m}). The recursive
formula for Romberg integration formulas:
R_{m+1}(h) =R_{m}(h)
+ _{}
·
Numerical algorithm
There
are two modifications of the Romberg integration algorithm: with doubling the
step size h and with halving the step size. The first
modification is used when a limited number of data values is available. The
second modification is more preferable if the function y = f(x) is
available analytically and the data values (x_{k},y_{k}) can
be computed for any step size h. In order to compute the
higherorder integration rules of the integral of f(x) between x_{0}
_{}x _{}x_{n }up to the order m, the
number n must be matched with m as: n = 2^{m1}.
In this case, there are sufficient number of points to compute
integration rules of lower order with larger step sizes: h, 2h, 4h, 8h,
…, (m1)h. The numerical approximations for the integrals can be
arranged in a table of recursive integration formulas starting with the
simplest approximation R_{1}(h) (composite trapezoidal
rule):
step size 
R_{1} 
R_{2} 
R_{3} 
R_{4} 
R_{5} 
h 
R_{1}(h) 




2h 
R_{1}(2h) 
R_{2}(h) 



4h 
R_{1}(4h) 
R_{2}(2h) 
R_{3}(h) 


8h_{} 
R_{1}(8h) 
R_{2}(4h) 
R_{3}(2h) 
R_{4}(h) 

16h 
R_{1}(16h) 
R_{2}(8h) 
R_{3}(4h) 
R_{4}(2h) 
R_{5}(h) 
The
diagonal entries are values of higherorder integration rules for the integral
of f(x) between x_{0} _{}x _{}x_{n }. The higherorder approximation R_{k}(h)
has the truncation error O(h^{2k}). If h is
small, the truncation error rapidly decreases with larger k.
However, the rounding error is constant with larger value of k. As
a result, no further increase in accuracy of numerical integration can be
obtained after some number m _{}M. The algorithm can be terminated when the relative
error falls below an error tolerance.
Example:
The figure
below presents the numerical approximations R_{1}(h) and
R_{2}(h) for the integral of the current I(t) with
the step size h = 10. Blue circles are found by the composite
Simpson rule R_{2}(h), green pluses are obtained by
composite trapezoidal rule R_{1}(h), and the exact
integral of I(t) is shown by a red solid curve. The composite
Simpson rule R_{2}(h) is clearly more accurate than the
composite trapezoidal rule R_{1}(h).