A = tridiagK’1 (’12, 24, ’12),

B = tridiagK’1 (’6, 0, 6),

C = tridiagK’1 (2, 8, 2).

System (8.68) has size equal to 2(K ’ 1); eliminating the unknown p from

the second equation, we get the reduced system (of size K ’ 1)

A ’ BC’1 BT u = b1 . (8.69)

366 8. Polynomial Interpolation

Since B is skew-symmetric and A is symmetric and positive de¬nite (s.p.d.),

the matrix M = A ’ BC’1 BT is s.p.d. too. Using Cholesky factorization for

solving system (8.69) is impractical as C’1 is full. An alternative is thus the

conjugate gradient method (CG) supplied with a suitable preconditioner

as the spectral condition number of M is of the order of h’4 = K 4 .

We notice that computing the residual at each step k ≥ 0 requires solv-

ing a linear system whose right side is the vector BT u(k) , u(k) being the

current iterate of CG method, and whose coe¬cient matrix is matrix C.

This system can be solved using the Thomas algorithm (3.53) with a cost

of the order of K ¬‚ops.

The CG algorithm terminates in correspondence to the lowest value of k

for which r(k) 2 ¤ u b1 2 , where r(k) is the residual of system (8.69) and

u is the roundo¬ unit.

The results obtained running the CG method in the case of a uniform

partition of [0, 1] with K = 50 elements and setting ± = f = 1 are sum-

marized in Figure 8.17 (right), which shows the convergence histories of

the method in both nonpreconditioned form (denoted by “Non Prec.”) and

with SSOR preconditioner (denoted by “Prec.”), having set the relaxation

parameter ω = 1.95.

We notice that the CG method does not converge within K ’ 1 steps

due to the e¬ect of the rounding errors. Notice also the e¬ectiveness of the

SSOR preconditioner in terms of the reduction of the number of iterations.

However, the high computational cost of this preconditioner prompts us to

devise another choice. Looking at the structure of the matrix M a natural

preconditioner is M = A ’ BC’1 BT , where C is the diagonal matrix whose

K’1

entries are cii = j=1 |cij |. The matrix M is banded so that its inversion

requires a strongly reduced cost than for the SSOR preconditioner. More-

over, as shown in Table 8.6, using M provides a dramatic decrease of the

number of iterations to converge.

M

K Without Precond. SSOR

25 51 27 12

50 178 61 25

100 685 118 33

200 2849 237 34

TABLE 8.6. Number of iterations as a function of K

8.8.2 Geometric Reconstruction Based on Computer

Tomographies

A typical application of the algorithms presented in Section 8.7 deals with

the reconstruction of the three-dimensional structure of internal organs of

human body based on computer tomographies (CT).

8.8 Applications 367

FIGURE 8.18. Cross-section of a blood vessel (left) and an associated character-

istic polygon using 16 points Pi (right)

The CT usually provides a sequence of images which represent the sections

of an organ at several horizontal planes; as a convention, we say that the

CT produces sections of the x, y plane in correspondance of several values of

z. The result is analogous to what we would get by sectioning the organ at

di¬erent values of z and taking the picture of the corresponding sections.

Obviously, the great advantage in using the CT is that the organ under

investigation can be visualized without being hidden by the neighboring

ones, as happens in other kinds of medical images, e.g., angiographies.

The image that is obtained for each section is coded into a matrix of

pixels (abbreviation of pictures elements) in the x, y plane; a certain value

is associated with each pixel expressing the level of grey of the image at

that point. This level is determined by the density of X rays which are

collected by a detector after passing through the human body. In practice,

the information contained in a CT at a given value of z is expressed by a

set of points (xi , yi ) which identify the boundary of the organ at z.

To improve the diagnostics it is often useful to reconstruct the three-

dimensional structure of the organ under examination starting from the

sections provided by the CT. With this aim, it is necessary to convert the

information coded by pixels into a parametric representation which can be

expressed by suitable functions interpolating the image at some signi¬cant

points on its boundary. This reconstruction can be carried out by using the

methods described in Section 8.7 as shown in Figure 8.19.

A set of curves like those shown in Figure 8.19 can be suitably stacked to

provide an overall three-dimensional view of the organ under examination.

368 8. Polynomial Interpolation

(c)

(a)

(b)

FIGURE 8.19. Reconstruction of the internal vessel of Figure 8.18 using di¬erent

interpolating splines with the same characteristic polygon: (a) B´zier curves, (b)

e

parametric splines and (c) parametric B-splines

8.9 Exercises

1. Prove that the characteristic polynomials li ∈ Pn de¬ned in (8.3) form a

basis for Pn .

2. An alternative approach to the method in Theorem 8.1, for constructing

the interpolating polynomial, consists of directly enforcing the n + 1 in-

terpolation constraints on Πn and then computing the coe¬cients ai . By

doing so, we end up with a linear system Xa= y, with a = (a0 , . . . , an )T ,

y = (y0 , . . . , yn )T and X = [xj ]. X is called Vandermonde matrix. Prove

i

that X is nonsingular if the nodes xi are distinct.

(xi ’ xj ) by recursion on n.]

[Hint: show that det(X)=

0¤j<i¤n

n

(xi ’ xj ) where ωn+1 is the nodal polynomial

3. Prove that ωn+1 (xi ) =

j=0

j=i

(8.6). Then, check (8.5).

4. Provide an estimate of ωn+1 ∞ , in the cases n = 1 and n = 2, for a

distribution of equally spaced nodes.

5. Prove that

(n ’ 1)!hn’1 |(x ’ xn’1 )(x ’ xn )| ¤ |ωn+1 (x)| ¤ n!hn’1 |(x ’ xn’1 )(x ’ xn )|,

where n is even, ’1 = x0 < x1 < . . . < xn’1 < xn = 1, x ∈ (xn’1 , xn ) and

h = 2/n.

8.9 Exercises 369

[Hint : let N = n/2 and show ¬rst that

ωn+1 (x) = (x + N h)(x + (N ’ 1)h) . . . (x + h)x

(8.70)

(x ’ h) . . . (x ’ (N ’ 1)h)(x ’ N h).

Then, take x = rh with N ’ 1 < r < N .]

6. Under the assumptions of Exercise 5, show that |ωn+1 | is maximum if

x ∈ (xn’1 , xn ) (notice that |ωn+1 | is an even function).

[Hint : use (8.70) to prove that |ωn+1 (x + h)/ωn+1 (x)| > 1 for any x ∈

(0, xn’1 ) with x not coinciding with any interpolation node.]

7. Prove the recursive relation (8.19) for Newton divided di¬erences.

8. Determine an interpolating polynomial Hf ∈ Pn such that

(Hf )(k) (x0 ) = f (k) (x0 ), k = 0, . . . , n,

and check that

n

f (j) (x0 )

(x ’ x0 )j ,

Hf (x) =

j!

j=0

that is, the Hermite interpolating polynomial on one node coincides with

the Taylor polynomial.

9. Given the following set of data

f0 = f (’1) = 1, f1 = f (’1) = 1, f2 = f (1) = 2, f3 = f (2) = 1 ,

prove that the Hermite-Birko¬ interpolating polynomial H3 does not exist

for them.

[Solution : letting H3 (x) = a3 x3 + a2 x2 + a1 x + a0 , one must check that

the matrix of the linear system H3 (xi ) = fi for i = 0, . . . , 3 is singular.]

10. Check that any sk ∈ Sk [a, b] admits a representation of the form

g

k

ci (x ’ xi )k ,

i

sk (x) = bi x + +

i=0 i=1

that is, 1, x, x2 , . . . , xk , (x ’ x1 )k , . . . , (x ’ xg )k form a basis for Sk [a, b].

+ +

11. Prove Property 8.2 and check its validity even in the case where the spline

s satis¬es conditions of the form s (a) = f (a), s (b) = f (b).

[Hint: start from

xi

b n

f (x) ’ s (x) s (x)dx = f (x) ’ s (x) s dx

i=1

a xi’1

and integrate by parts twice.]

370 8. Polynomial Interpolation

x2 x4 x6

12. Let f (x) = cos(x) = 1 ’ ’

+ + . . . ; then, consider the following

2! 4! 6!

rational approximation

a0 + a2 x2 + a4 x4

r(x) = , (8.71)

1 + b2 x2

called the Pad´ approximation. Determine the coe¬cients of r in such a

e

way that

f (x) ’ r(x) = γ8 x8 + γ10 x10 + . . .

[Solution: a0 = 1, a2 = ’7/15, a4 = 1/40, b2 = 1/30.]

13. Assume that the function f of the previous exercise is known at a set of n

equally spaced points xi ∈ (’π/2, π/2) with i = 0, . . . , n. Repeat Exercise

12, determining, by using MATLAB, the coe¬cients of r in such a way

i=0 |f (xi ) ’ r(xi )| is minimized. Consider the cases

n 2

that the quantity

n = 5 and n = 10.

9

Numerical Integration

In this chapter we present the most commonly used methods for numer-

ical integration. We will mainly consider one-dimensional integrals over

bounded intervals, although in Sections 9.8 and 9.9 an extension of the tech-

niques to integration over unbounded intervals (or integration of functions

with singularities) and to the multidimensional case will be considered.

9.1 Quadrature Formulae

Let f be a real integrable function over the interval [a, b]. Computing ex-

b

plicitly the de¬nite integral I(f ) = a f (x)dx may be di¬cult or even

impossible. Any explicit formula that is suitable for providing an approxi-

mation of I(f ) is said to be a quadrature formula or numerical integration

formula.

An example can be obtained by replacing f with an approximation fn ,

depending on the integer n ≥ 0, then computing I(fn ) instead of I(f ).

Letting In (f ) = I(fn ), we have

b

n ≥ 0.

In (f ) = fn (x)dx, (9.1)

a

The dependence on the end points a, b is always understood, so we write

In (f ) instead of In (f ; a, b).

372 9. Numerical Integration

If f ∈ C 0 ([a, b]), the quadrature error En (f ) = I(f ) ’ In (f ) satis¬es

b

|En (f )| ¤ |f (x) ’ fn (x)|dx ¤ (b ’ a) f ’ fn ∞.

a

Therefore, if for some n, f ’ fn ∞ < µ, then |En (f )| ¤ µ(b ’ a).

The approximant fn must be easily integrable, which is the case if, for

example, fn ∈ Pn . In this respect, a natural approach consists of using

fn = Πn f , the interpolating Lagrange polynomial of f over a set of n + 1

distinct nodes {xi }, with i = 0, . . . , n. By doing so, from (9.1) it follows

that

b

n

In (f ) = f (xi ) li (x)dx, (9.2)

i=0 a

where li is the characteristic Lagrange polynomial of degree n associated

with node xi (see Section 8.1). We notice that (9.2) is a special instance of

the following quadrature formula

n

In (f ) = ±i f (xi ), (9.3)

i=0

b

where the coe¬cients ±i of the linear combination are given by a li (x)dx.

Formula (9.3) is a weighted sum of the values of f at the points xi , for

i = 0, . . . , n. These points are said to be the nodes of the quadrature

formula, while the numbers ±i ∈ R are its coe¬cients or weights. Both

weights and nodes depend in general on n; again, for notational simplicity,

this dependence is always understood.

Formula (9.2), called the Lagrange quadrature formula, can be generalized

to the case where also the values of the derivative of f are available. This

leads to the Hermite quadrature formula (see Section 9.5)

1 n

±ik f (k) (xi )

In (f ) = (9.4)

k=0 i=0

where the weights are now denoted by ±ik .

Both (9.2) and (9.4) are interpolatory quadrature formulae, since the

function f has been replaced by its interpolating polynomial (Lagrange

and Hermite polynomials, respectively). We de¬ne the degree of exactness

of a quadrature formula as the maximum integer r ≥ 0 for which

∀f ∈ Pr .

In (f ) = I(f ),

Any interpolatory quadrature formula that makes use of n + 1 distinct

nodes has degree of exactness equal to at least n. Indeed, if f ∈ Pn , then

9.2 Interpolatory Quadratures 373

Πn f = f and thus In (Πn f ) = I(Πn f ). The converse statement is also true,

that is, a quadrature formula using n + 1 distinct nodes and having degree

of exactness equal at least to n is necessarily of interpolatory type (for the

proof see [IK66], p. 316).

As we will see in Section 10.2, the degree of exactness of a Lagrange

quadrature formula can be as large as 2n + 1 in the case of the so-called

Gaussian quadrature formulae.

9.2 Interpolatory Quadratures

We consider three remarkable instances of formula (9.2), corresponding to

n = 0, 1 and 2.

9.2.1 The Midpoint or Rectangle Formula

This formula is obtained by replacing f over [a, b] with the constant function

equal to the value attained by f at the midpoint of [a, b] (see Figure 9.1,

left). This yields

a+b

I0 (f ) = (b ’ a)f (9.5)

2

with weight ±0 = b ’ a and node x0 = (a + b)/2. If f ∈ C 2 ([a, b]), the

quadrature error is

b’a

h3

(9.6)

E0 (f ) = f (ξ), h = ,

3 2

where ξ lies within the interval (a, b).

f (x) f (x)

x x

xm’1

x0

a x0 b xk

FIGURE 9.1. The midpoint formula (left); the composite midpoint formula

(right)

Indeed, expanding f in a Taylor™s series around c = (a + b)/2 and trun-

cating at the second-order, we get

f (x) = f (c) + f (c)(x ’ c) + f (·(x))(x ’ c)2 /2,

374 9. Numerical Integration

from which, integrating on (a, b) and using the mean-value theorem, (9.6)

follows. From this, it turns out that (9.5) is exact for constant and a¬ne

functions (since in both cases f (ξ) = 0 for any ξ ∈ (a, b)), so that the

midpoint rule has degree of exactness equal to 1.

It is worth noting that if the width of the integration interval [a, b] is

not su¬ciently small, the quadrature error (9.6) can be quite large. This

drawback is common to all the numerical integration formulae that will

be described in the three forthcoming sections and can be overcome by

resorting to their composite counterparts as discussed in Section 9.4.

Suppose now that we approximate the integral I(f ) by replacing f over

[a, b] with its composite interpolating polynomial of degree zero, constructed

on m subintervals of width H = (b ’ a)/m, for m ≥ 1 (see Figure 9.1,

right). Introducing the quadrature nodes xk = a + (2k + 1)H/2, for k =

0, . . . , m ’ 1, we get the composite midpoint formula

m’1

m ≥ 1.

I0,m (f ) = H f (xk ), (9.7)

k=0

The quadrature error E0,m (f ) = I(f ) ’ I0,m (f ) is given by

b’a 2 b’a

(9.8)

E0,m (f ) = H f (ξ), H =

24 m

provided that f ∈ C 2 ([a, b]) and where ξ ∈ (a, b). From (9.8) we conclude

that (9.7) has degree of exactness equal to 1; (9.8) can be proved by recalling

(9.6) and using the additivity of integrals. Indeed, for k = 0, . . . , m ’ 1 and

ξk ∈ (a + kH, a + (k + 1)H),

m’1 m’1

H2 b ’ a b’a 2

3

E0,m (f ) = f (ξk )(H/2) /3 = f (ξk ) = H f (ξ).

24 m 24

k=0 k=0

The last equality is a consequence of the following theorem, that is applied

letting u = f and δj = 1 for j = 0, . . . , m ’ 1.

Theorem 9.1 (discrete mean-value theorem) Let u ∈ C 0 ([a, b]) and

let xj be s + 1 points in [a, b] and δj be s + 1 constants, all having the same

sign. Then there exists · ∈ [a, b] such that

s s

δj u(xj ) = u(·) δj . (9.9)

j=0 j=0

¯

Proof. Let um = minx∈[a,b] u(x) = u(¯) and uM = maxx∈[a,b] u(x) = u(x),

x

¯