m 4000 5600 250 50

3.24e ’ 10 3.30e ’ 10 2.98e ’ 10 3.58e ’ 11 3.18e ’ 10

Err

20007 29013 2519 5772 3540

flops

√

1/β

1 ’ K 2 x2 dx, with

TABLE 9.12. Methods for approximating I(f ) = 4π±

√ 0

±2 = (3 ’ 2 2)/100, β = 10 and K 2 = β 2 (1 ’ ±2 β 2 )

5

10

1

0.8

4

10

0.6

0.4

3

10

0.2

0

2

10

0.1

0.06 0.08

0.04

0 0.02

’0.2

0.02 0.04 0.06 0.1

0.08

0

FIGURE 9.7. Distribution of quadrature nodes (left); integration stepsize density

in the approximation of integral (9.62) (right)

9.10.2 Computation of the Wind Action on a Sailboat Mast

Let us consider the sailboat schematically drawn in Figure 9.8 (left) and

subject to the action of the wind force. The mast, of length L, is denoted by

the straight line AB, while one of the two shrouds (strings for the side sti¬-

ening of the mast) is represented by the straight line BO. Any in¬nitesimal

element of the sail transmits to the corresponding element of length dx of

the mast a force of magnitude equal to f (x)dx. The change of f along with

the height x, measured from the point A (basis of the mast), is expressed

by the following law

±x ’γx

f (x) = e ,

x+β

where ±, β and γ are given constants.

9.10 Applications 411

The resultant R of the force f is de¬ned as

L

f (x)dx ≡ I(f ),

R= (9.63)

0

and is applied at a point at distance equal to b (to be determined) from

the basis of the mast.

B

B

mast

shroud

wind

1111111111

0000000000

direction

0000000000

1111111111 L

R

1111111111

0000000000 f dx

T

dx

T

b

M

A

O

O A

H

V

FIGURE 9.8. Schematic representation of a sailboat (left); forces acting on the

mast (right)

Computing R and the distance b, given by b = I(xf )/I(f ), is crucial for the

structural design of the mast and shroud sections. Indeed, once the values

of R and b are known, it is possible to analyze the hyperstatic structure

mast-shroud (using for instance the method of forces), thus allowing for the

computation of the reactions V , H and M at the basis of the mast and the

traction T that is transmitted by the shroud, and are drawn in Figure 9.8

(right). Then, the internal actions in the structure can be found, as well as

the maximum stresses arising in the mast AB and in the shroud BO, from

which, assuming that the safety veri¬cations are satis¬ed, one can ¬nally

design the geometrical parameters of the sections of AB and BO.

For the approximate computation of R we have considered the compos-

ite midpoint, trapezoidal and Cavalieri-Simpson rules, denoted henceforth

by (MP), (TR) and (CS), and, for a comparison, the adaptive Cavalieri-

Simpson quadrature formula introduced in Section 9.7.2 and denoted by

(AD). Since a closed-form expression for the integral (9.63) is not available,

the composite rules have been applied taking mk = 2k uniform partitions

of [0, L], with k = 0, . . . , 15.

We have assumed in the numerical experiments ± = 50, β = 5/3 and

γ = 1/4 and we have run Program 77 taking tol=10’4 and hmin=10’3 .

The sequence of integrals computed using the composite formulae has been

stopped at k = 12 (corresponding to mk = 212 = 4096) since the remaining

412 9. Numerical Integration

0

10

’1

10

’2

10

’3

(TR)

10

’4

10 (PM)

’5

10

’6

10

(CS)

’7

10

’8

10 (AD)

’9

10

0 20 40 60 80 100 120

FIGURE 9.9. Relative errors in the approximate computation of the integral

(±xe’γx )/(x + β)dx

L

0

elements, in the case of formula CS, di¬er among them only up to the last

signi¬cant ¬gure. Therefore, we have assumed as the exact value of I(f )

(CS)

the outcome I12 = 100.0613683179612 provided by formula CS.

(CS)

We report in Figure 9.9 the logarithmic plots of the relative error |I12 ’

Ik |/I12 , for k = 0, . . . , 7, Ik being the generic element of the sequence for

the three considered formulae. As a comparison, we also display the graph

of the relative error in the case of formula AD, applied on a number of

(nonuniform) partitions equivalent to that of the composite rules.

Notice how, for the same number of partitions, formula AD is more accu-

rate, with a relative error of 2.06 · 10’7 obtained using 37 (nonuniform)

partitions of [0, L]. Methods PM and TR achieve a comparable accuracy

employing 2048 and 4096 uniform subintervals, respectively, while formula

CS requires about 64 partitions. The e¬ectiveness of the adaptivity pro-

cedure is demonstrated by the plots in Figure 9.10, which show, together

with the graph of f , the distribution of the quadrature nodes (left) and the

function ∆h (x) (right) that expresses the (piecewise constant) density of

the integration stepsize h, de¬ned as the inverse of the value of h over each

active interval A (see Section 9.7.2).

Notice also the high value of ∆h at x = 0, where the derivatives of f are

maximum.

9.11 Exercises

1. Let E0 (f ) and E1 (f ) be the quadrature errors in (9.6) and (9.12). Prove

that |E1 (f )| 2|E0 (f )|.

2. Check that the error estimates for the midpoint, trapezoidal and Cavalieri-

Simpson formulae, given respectively by (9.6), (9.12) and (9.16), are special

instances of (9.19) or (9.20). In particular, show that M0 = 2/3, K1 =

9.11 Exercises 413

20 30

25

15

20

10

15

5

10

0

5

’5 0

0 2 4 6 8 10 0 2 4 6 8 10

FIGURE 9.10. Distribution of quadrature nodes (left); integration step density

in the approximation of the integral 0 (±xe’γx )/(x + β)dx (right)

L

’1/6 and M2 = ’4/15 and determine, using the de¬nition, the degree of

exactness r of each formula.

b

[Hint: ¬nd r such that In (xk ) = a xk dx, for k = 0, . . . , r, and In (xj ) =

bj

x dx, for j > r.]

a

n

3. Let In (f ) = k=0 ±k f (xk ) be a Lagrange quadrature formula on n + 1

nodes. Compute the degree of exactness r of the formulae:

(a) I2 (f ) = (2/3)[2f (’1/2) ’ f (0) + 2f (1/2)],

(b) I4 (f ) = (1/4)[f (’1) + 3f (’1/3) + 3f (1/3) + f (1)].

Which is the order of in¬nitesimal p for (a) and (b)?

[Solution: r = 3 and p = 5 for both I2 (f ) and I4 (f ).]

4. Compute df [x0 , . . . , xn , x]/dx by checking (9.22).

[Hint: proceed by computing directly the derivative at x as an incremental

ratio, in the case where only one node x0 exists, then upgrade progressively

the order of the divided di¬erence.]

√

1

5. Let Iw (f ) = 0 w(x)f (x)dx with w(x) = x, and consider the quadrature

formula Q(f ) = af (x1 ). Find a and x1 in such a way that Q has maximum

degree of exactness r.

[Solution: a = 2/3, x1 = 3/5 and r = 1.]

6. Let us consider the quadrature formula Q(f ) = ±1 f (0) + ±2 f (1) + ±3 f (0)

1

for the approximation of I(f ) = 0 f (x)dx, where f ∈ C 1 ([0, 1]). Determine

the coe¬cients ±j , for j = 1, 2, 3 in such a way that Q has degree of

exactness r = 2.

[Solution: ±1 = 2/3, ±2 = 1/3 and ±3 = 1/6.]

7. Apply the midpoint, trapezoidal and Cavalieri-Simpson composite rules to

approximate the integral

1

|x|ex dx,

’1

and discuss their convergence as a function of the size H of the subintervals.

414 9. Numerical Integration

1

8. Consider the integral I(f ) = 0 ex dx and estimate the minimum number m

of subintervals that is needed for computing I(f ) up to an absolute error

¤ 5 · 10’4 using the composite trapezoidal (TR) and Cavalieri-Simpson

(CS) rules. Evaluate in both cases the absolute error Err that is actually

made.

[Solution: for TR, we have m = 17 and Err = 4.95 · 10’4 , while for CS,

m = 2 and Err = 3.70 · 10’5 .]

9. Consider the corrected trapezoidal formula (9.30) and check that |E1 (f )|

corr

corr

4|E2 (f )|, where E1 (f ) and E2 (f ) are de¬ned in (9.31) and (9.16), re-

spectively.

10. Compute, with an error less than 10’4 , the following integrals:

∞

sin(x)/(1 + x4 )dx;

(a) 0

∞

e’x (1 + x)’5 dx;

(b) 0

2

∞

cos(x)e’x dx.

(c) ’∞

11. Use the reduction midpoint and trapezoidal formulae for computing the

y

dxdy over the domain „¦ = (0, 1)2 . Run

double integral I(f ) = „¦

(1 + xy)

Programs 78 and 79 with M = 2i , for i = 0, . . . , 10 and plot in log-scale

the absolute error in the two cases as a function of M . Which method is

the most accurate? How many functional evaluations are needed to get an

(absolute) accuracy of the order of 10’6 ?

[Solution: the exact integral is I(f ) = log(4) ’ 1, and almost 2002 = 40000

functional evaluations are needed.]

10

Orthogonal Polynomials in

Approximation Theory

Trigonometric polynomials, as well as other orthogonal polynomials like

Legendre™s and Chebyshev™s, are widely employed in approximation theory.

This chapter addresses the most relevant properties of orthogonal poly-

nomials, and introduces the transforms associated with them, in particular

the discrete Fourier transform and the FFT, but also the Zeta and Wavelet

transforms.

Application to interpolation, least-squares approximation, numerical dif-

ferentiation and Gaussian integration are addressed.

10.1 Approximation of Functions by Generalized

Fourier Series

Let w = w(x) be a weight function on the interval (’1, 1), i.e., a nonneg-

ative integrable function in (’1, 1). Let us denote by {pk , k = 0, 1, . . . } a

system of algebraic polynomials, with pk of degree equal to k for each k,

mutually orthogonal on the interval (’1, 1) with respect to w. This means

that

1

pk (x)pm (x)w(x)dx = 0 if k = m.

’1

416 10. Orthogonal Polynomials in Approximation Theory

1 1/2

Set (f, g)w = ’1 f (x)g(x)w(x)dx and f w = (f, f )w ; (·, ·)w and · w

are respectively the scalar product and the norm for the function space

1

f : (’1, 1) ’ R, f 2 (x)w(x)dx < ∞ .

L2 L2 (’1, 1)

= = (10.1)

w w

’1

For any function f ∈ L2 the series

w

+∞

(f, pk )w

Sf = fk p k , with fk = ,

pk 2 w

k=0

is called the generalized Fourier series of f , and fk is the k-th Fourier

coe¬cient. As is well-known, Sf converges in average (or in the sense of

L2 ) to f . This means that, letting for any integer n

w

n

fn (x) = fk pk (x) (10.2)

k=0

(fn ∈ Pn is the truncation of order n of the generalized Fourier series of

f ), the following convergence result holds

f ’ fn

lim = 0.

w

n’+∞

Thanks to Parseval™s equality, we have

+∞

2 2 2

f = fk pk

w w

k=0

+∞

and, for any n, f ’fn 2 = k=n+1 fk pk 2 is the square of the remainder

2

w w

of the generalized Fourier series.

The polynomial fn ∈ Pn satis¬es the following minimization property

f ’ fn = min f ’ q w. (10.3)

w

q∈Pn

+∞

Indeed, since f ’ fn = k=n+1 fk pk , the property of orthogonality of

polynomials {pk } implies (f ’ fn , q)w = 0 ∀q ∈ Pn . Then, the Cauchy-

Schwarz inequality (8.29) yields

f ’ fn = (f ’ fn , f ’ fn )w = (f ’ fn , f ’ q)w + (f ’ fn , q ’ fn )w

2

w

= (f ’ fn , f ’ q)w ¤ f ’ fn f ’q ∀q ∈ Pn ,

w,

w

and (10.3) follows since q is arbitrary in Pn . In such a case, we say that fn

is the orthogonal projection of f over Pn in the sense of L2 . It is therefore

w

interesting to compute the coe¬cients fk of fn . As will be seen in later