(I)

S A

a ± b

(II)

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

4

tan’1 (10x)dx

I(f ) =

’3

= 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

80

2

70

1.5

60

1

50

0.5

40

0

30

’0.5

20

’1

10

’1.5

0

’2

’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));

S2=S2+(h/3)*(f(i,3)+4*f(i,4)+f(i,5));

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;

f(i,1)=f(i-1,1);f(i,3)=f(i-1,2);f(i,5)=f(i-1,3);len=abs(beta(i)-alfa(i));

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

end

end

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

φ(x)

0 ¤ µ < 1,

f (x) = ,

(x ’ a)µ

where φ is a function whose absolute value is bounded by M . Then

b

(b ’ a)1’µ

1

|I(f )| ¤ M lim dx = M .

(x ’ a)µ 1’µ

+

t’a

t

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)

p≥0

φ(x) = ¦p (x) + φ (ξ(x)), (9.47)

(p + 1)!

p

φ(k) (a)(x ’ a)k /k!. Then

where ¦p (x) = k=0

a+µ

p

µ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

as

µp+2’µ

|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

grows.

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

n+2

b ’ a ’ µ |Mn | b’a’µ

|E2 | ¤ M (n+2)

(µ) = δ/2, (9.49)

(n + 2)! nn+3 m

where

dn+2 φ(x)

M (n+2)

(µ) = max .

(x ’ a)µ

dxn+2

a+µ¤x¤b

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)

π/2

cos(x)

√ dx.

I(f ) = (9.50)

x

0

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)

k=0

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,

1/4

0.018 π 5

µ’9/2

’µ

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

p

(b ’ a)k φ(k) (a)

I2 = (b ’ a) 1’µ

. (9.51)

k!(k + 1 ’ µ)

k=0

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

t

lim f (x)dx

t’+∞

a

is taken as being the value of the singular integral

t

∞

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

t’+∞

a

a

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,

x’+∞

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

∞

c

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 ) =

0

∞

∞

cos2 (x)e’x dx ¤ e’x dx = e’c ; re-

For any given c > 0, we have I2 =

c

c

’c

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

1/2

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 =

’c

3/5 ’ e (cos (c) ’ (sin(2c) + 2 cos(2c))/5) 0.599842, with an absolute error of

2

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,