<< . .

. 53
( : 95)



. . >>

b1 ∈ RK’1 is the vector of components equal to h4 f /±, while
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
¯

<< . .

. 53
( : 95)



. . >>