<< . .

. 57
( : 95)

. . >>

a ± b

S ± A± N
a b

FIGURE 9.3. Distribution of the integration intervals at the generic step of the
adaptive algorithm and updating of the integration grid

2. If the local error test (9.45) fails, then:
(j) A is halved, and the new active interval is set to A = [±, ± ] with
± = (± + β)/2 (corresponding to the path (II) in Figure 9.3);
(jj) we let N ← N ∪ [± , β], β ← ± ;
(jjj) a new error estimate is provided.

In order to prevent the algorithm from generating too small stepsizes, it
is convenient to monitor the width of A and warn the user, in case of
an excessive reduction of the steplength, about the presence of a possible
singularity in the integrand function (see Section 9.8).

Example 9.9 Let us employ Cavalieri-Simpson adaptive integration for com-
puting the integral

tan’1 (10x)dx
I(f ) =
= 4tan’1 (40) + 3tan’1 (’30) ’ (1/20) log(16/9) 1.54201193.

Running Program 77 with tol = 10’4 and hmin = 10’3 yields an approximation
of the integral with an absolute error of 2.104 · 10’5 . The algorithm performs
77 functional evaluations, corresponding to partitioning the interval [a, b] into
38 nonuniform subintervals. We notice that the corresponding composite formula
with uniform stepsize would have required 128 subintervals with an absolute error
of 2.413 · 10’5 .
In Figure 9.4 (left) we show, together with the plot of the integrand function,
the distribution of the quadrature nodes as a function of x, while on the right
the integration step density (piecewise constant) ∆h (x) is shown, de¬ned as the
inverse of the step size h over each active interval A. Notice the high value attained
by ∆h at x = 0, where the derivative of the integrand function is maximum. •

The adaptive algorithm described above is implemented in Program 77.
Among the input parameters, hmin is the minimum admissible value of
the integration steplength. In output the program returns the approximate
9.7 Automatic Integration 397









’3 ’2 ’1 0 1 2 3 4
’3 ’2 ’1 0 1 2 3 4

FIGURE 9.4. Distribution of quadrature nodes (left); density of the integration
stepsize in the approximation of the integral of Example 9.9 (right)

value of the integral integ, the total number of functional evaluations nfv
and the set of integration points xfv.

Program 77 - simpadpt : Adaptive Cavalieri-Simpson formula

function [integ,xfv,nfv]=simpadpt(a,b,tol,fun,hmin);
integ=0; level=0; i=1; alfa(i)=a; beta(i)=b;
step=(beta(i)-alfa(i))/4; nfv=0;
for k=1:5, x=a+(k-1)*step; f(i,k)=eval(fun); nfv=nfv+1; end
while (i > 0),
S=0; S2=0; h=(beta(i)-alfa(i))/2; S=(h/3)*(f(i,1)+4*f(i,3)+f(i,5));
h=h/2; S2=(h/3)*(f(i,1)+4*f(i,2)+f(i,3));
tolrv=tol*(beta(i)-alfa(i))/(b-a); errrv=abs(S-S2)/10;
if (errrv > tolrv)
i=i+1; alfa(i)=alfa(i-1); beta(i)=(alfa(i-1)+beta(i-1))/2;
if (len >= hmin),
if (len <= 11*hmin)
disp(™ Steplength close to hmin ™),
str=sprintf(™The approximate integral is %12.7e™,integ);disp(str),end;
step=len/4; x=alfa(i)+step; f(i,2)=eval(fun);
nfv=nfv+1; x=beta(i)-step; f(i,4)=eval(fun); nfv=nfv+1;
else, xfv=xfv™; disp(™ Too small steplength ™)
str=sprintf(™The approximate integral is %12.7e™,integ);
disp(str), return
end, else
integ=integ+S2; level=level+1; if (level==1),
for k=1:5, xfv(k)=alfa(i)+(k-1)*h; end; ist=5;
else, for k=1:4, xfv(ist+k)=alfa(i)+k*h; end; ist=ist+4; end;
if (beta(i)==b), xfv=xfv™;
str=sprintf(™The approximate integral is %12.7e™,integ);
disp(str), return, end; i=i-1; alfa(i)=beta(i+1);
f(i,1)=f(i+1,5); f(i,3)=f(i,4); step=abs(beta(i)-alfa(i))/4;
x=alfa(i)+step; f(i,2)=eval(fun); nfv=nfv+1; x=beta(i)-step;
f(i,4)=eval(fun); nfv=nfv+1;
398 9. Numerical Integration


9.8 Singular Integrals
In this section we extend our analysis to deal with singular integrals, arising
when f has ¬nite jumps or is even in¬nite at some point. Besides, we
will consider the case of integrals of bounded functions over unbounded
intervals. We brie¬‚y address the most relevant numerical techniques for
properly handling these integrals.

9.8.1 Integrals of Functions with Finite Jump Discontinuities
Let c be a known point within [a, b] and assume that f is a continuous and
bounded function in [a, c) and (c, b], with ¬nite jump f (c+ ) ’ f (c’ ). Since
b c b
I(f ) = f (x)dx = f (x)dx + f (x)dx, (9.46)
a a c

any integration formula of the previous sections can be used on [a, c’ ] and
[c+ , b] to furnish an approximation of I(f ). We proceed similarly if f admits
a ¬nite number of jump discontinuities within [a, b].
When the position of the discontinuity points of f is not known a priori,
a preliminary analysis of the graph of the function should be carried out.
Alternatively, one can resort to an adaptive integrator that is able to detect
the presence of discontinuities when the integration steplength falls below
a given tolerance (see Section 9.7.2).

9.8.2 Integrals of In¬nite Functions
Let us deal with the case in which limx’a+ f (x) = ∞; similar consider-
ations hold when f is in¬nite as x ’ b’ , while the case of a point of
singularity c internal to the interval [a, b] can be recast to one of the pre-
vious two cases owing to (9.46). Assume that the integrand function is of
the form
0 ¤ µ < 1,
f (x) = ,
(x ’ a)µ
where φ is a function whose absolute value is bounded by M . Then
(b ’ a)1’µ
|I(f )| ¤ M lim dx = M .
(x ’ a)µ 1’µ
9.8 Singular Integrals 399

Suppose we wish to approximate I(f ) up to a ¬xed tolerance δ. For this, let
us describe the following two methods (for further details, see also [IK66],
Section 7.6, and [DR75], Section 2.12 and Appendix 1).

Method 1. For any µ such that 0 < µ < (b ’ a), we write the singular
integral as I(f ) = I1 + I2 , where
a+µ b
φ(x) φ(x)
I1 = dx, I2 = dx.
(x ’ a)µ (x ’ a)µ
a a+µ

The computation of I2 is not troublesome. After replacing φ by its p-th
order Taylor™s expansion around x = a, we obtain

(x ’ a)p+1 (p+1)
φ(x) = ¦p (x) + φ (ξ(x)), (9.47)
(p + 1)!
φ(k) (a)(x ’ a)k /k!. Then
where ¦p (x) = k=0

µk φ(k) (a) 1
(x ’ a)p+1’µ φ(p+1) (ξ(x))dx.
I1 = µ1’µ +
k!(k + 1 ’ µ) (p + 1)!
k=0 a

Replacing I1 by the ¬nite sum, the corresponding error E1 can be bounded
|E1 | ¤ max |φ(p+1) (x)|, p ≥ 0. (9.48)
(p + 1)!(p + 2 ’ µ) a¤x¤a+µ

For ¬xed p, the right side of (9.48) is an increasing function of µ. On the
other hand, taking µ < 1 and assuming that the successive derivatives of φ
do not grow too much as p increases, the same function is decreasing as p
Let us next approximate I2 using a composite Newton-Cotes formula
with m subintervals and n quadrature nodes for each subinterval, n being
an even integer. Recalling (9.26) and aiming at equidistributing the error
δ between I1 and I2 , it turns out that
b ’ a ’ µ |Mn | b’a’µ
|E2 | ¤ M (n+2)
(µ) = δ/2, (9.49)
(n + 2)! nn+3 m
dn+2 φ(x)
M (n+2)
(µ) = max .
(x ’ a)µ

The value of the constant M(n+2) (µ) grows rapidly as µ tends to zero; as
a consequence, (9.49) might require such a large number of subintervals
mµ = m(µ) to make the method at hand of little practical use.
400 9. Numerical Integration

Example 9.10 Consider the singular integral (known as the Fresnel integral)
√ dx.
I(f ) = (9.50)

Expanding the integrand function in a Taylor™s series around the origin and
applying the theorem of integration by series, we get

(’1)k 1
(π/2)2k+1/2 .
I(f ) =
(2k)! (2k + 1/2)

Truncating the series at the ¬rst 10 terms, we obtain an approximate value of
the integral equal to 1.9549.
Using the composite Cavalieri-Simpson formula, the a priori estimate (9.49)
yields, as µ tends to zero and letting n = 2, |M2 | = 4/15,
0.018 π 5
mµ .
δ 2
For δ = 10’4 , taking µ = 10’2 , it turns out that 1140 (uniform) subintervals are
needed, while for µ = 10’4 and µ = 10’6 the number of subintervals is 2 · 105
and 3.6 · 107 , respectively.
As a comparison, running Program 77 (adaptive integration with Cavalieri-
Simpson formula) with a = µ = 10’10 , hmin = 10’12 and tol = 10’4 , we
get the approximate value 1.955 for the integral at the price of 1057 functional
evaluations, which correspond to 528 nonuniform subdivisions of the interval

[0, π/2].

Method 2. Using the Taylor expansion (9.47) we obtain
b b
φ(x) ’ ¦p (x) ¦p (x)
I(f ) = dx + dx = I1 + I2 .
(x ’ a)µ (x ’ a)µ
a a

Exact computation of I2 yields
(b ’ a)k φ(k) (a)
I2 = (b ’ a) 1’µ
. (9.51)
k!(k + 1 ’ µ)

The integral I1 is, for p ≥ 0
b b
φ(p+1) (ξ(x))
(x ’ a)p+1’µ
I1 = dx = g(x)dx. (9.52)
(p + 1)!
a a

Unlike the case of method 1, the integrand function g does not blow up
at x = a, since its ¬rst p derivatives are ¬nite at x = a. As a consequence,
assuming we approximate I1 using a composite Newton-Cotes formula, it is
possible to give an estimate of the quadrature error, provided that p ≥ n+2,
for n ≥ 0 even, or p ≥ n + 1, for n odd.
9.8 Singular Integrals 401

Example 9.11 Consider again the singular Fresnel integral (9.50), and assume
we use the composite Cavalieri-Simpson formula for approximating I1 . We will
take p = 4 in (9.51) and (9.52). Computing I2 yields the value (π/2)1/2 (2 ’
(1/5)(π/2)2 + (1/108)(π/2)4 ) 1.9588. Applying the error estimate (9.26) with
n = 2 shows that only 2 subdivisions of [0, π/2] su¬ce for approximating I1 up
to an error δ = 10’4 , obtaining the value I1 ’0.0173. As a whole, method 2

returns for (9.50) the approximate value 1.9415.

9.8.3 Integrals over Unbounded Intervals
Let f ∈ C 0 ([a, +∞)); should it exist and be ¬nite, the following limit

lim f (x)dx

is taken as being the value of the singular integral

I(f ) = f (x)dx = lim f (x)dx. (9.53)

An analogous de¬nition holds if f is continuous over (’∞, b], while for a
function f : R ’ R, integrable over any bounded interval, we let
∞ c +∞
f (x)dx = f (x)dx + f (x)dx (9.54)
’∞ ’∞ c

if c is any real number and the two singular integrals on the right hand side
of (9.54) are convergent. This de¬nition is correct since the value of I(f )
does not depend on the choice of c.
A su¬cient condition for f to be integrable over [a, +∞) is that
∃ρ > 0, such that lim x1+ρ f (x) = 0,

that is, we require f to be in¬nitesimal of order > 1 with respect to 1/x
as x ’ ∞. For the numerical approximation of (9.53) up to a tolerance δ,
we consider the following methods, referring for further details to [DR75],
Chapter 3.

Method 1. To compute (9.53), we can split I(f ) as I(f ) = I1 + I2 , where

I1 = a f (x)dx and I2 = c f (x)dx.
The end-point c, which can be taken arbitrarily, is chosen in such a way
that the contribution of I2 is negligible. Precisely, taking advantage of the
asymptotic behavior of f , c is selected to guarantee that I2 equals a fraction
of the ¬xed tolerance, say, I2 = δ/2.
Then, I1 will be computed up to an absolute error equal to δ/2. This
ensures that the global error in the computation of I1 + I2 is below the
tolerance δ.
402 9. Numerical Integration

Example 9.12 Compute up to an error δ = 10’3 the integral

cos2 (x)e’x dx = 3/5.
I(f ) =

cos2 (x)e’x dx ¤ e’x dx = e’c ; re-
For any given c > 0, we have I2 =
quiring that e = δ/2, one gets c 7.6. Then, assuming we use the compos-
ite trapezoidal formula for approximating I1 , thanks to (9.27) with n = 1 and
M = max0¤x¤c |f (x)| 1.04, we obtain m ≥ M c3 /(6δ) = 277.
Program 72 returns the value I1 0.599905, instead of the exact value I1 =
3/5 ’ e (cos (c) ’ (sin(2c) + 2 cos(2c))/5) 0.599842, with an absolute error of

about 6.27 · 10’5 . The global numerical outcome is thus I1 + I2 0.600405, with
an absolute error with respect to I(f ) equal to 4.05 · 10’4 . •

Method 2. For any real number c, we let I(f ) = I1 + I2 , as for method 1,
then we introduce the change of variable x = 1/t in order to transform I2
into an integral over the bounded interval [0, 1/c]
1/c 1/c

f (t)t’2 dt =
I2 = g(t)dt. (9.55)
0 0

If g(t) is not singular at t = 0, (9.55) can be treated by any quadrature
formula introduced in this chapter. Otherwise, one can resort to the inte-
gration methods considered in Section 9.8.2.

Method 3. Gaussian interpolatory formulae are used, where the integra-
tion nodes are the zeros of Laguerre and Hermite orthogonal polynomials
(see Section 10.5).

9.9 Multidimensional Numerical Integration
Let „¦ be a bounded domain in R2 with a su¬ciently smooth boundary. We
consider the problem of approximating the integral I(f ) = „¦ f (x, y)dxdy,

<< . .

. 57
( : 95)

. . >>