ex ’ 1

0 0

where x = hν/KB T and ± = (8πKB )/(ch)3 1.16 · 10’16 [J][K ’4 ][m’3 ].

4

∞

We also let f (x) = x3 /(ex ’ 1) and I(f ) = 0 f (x)dx.

To approximate I(f ) up to a previously ¬xed absolute error ¤ δ, we com-

pare method 1. introduced in Section 9.8.3 with Gauss-Laguerre quadra-

tures.

In the case of method 1. we proceed as follows. For any a > 0 we let

∞

a

I(f ) = 0 f (x)dx + a f (x)dx and try to ¬nd a function φ such that

∞ ∞

δ

f (x)dx ¤ φ(x)dx ¤ , (10.89)

2

a a

∞

the integral φ(x)dx being “easy” to compute. Once the value of a

a

has been found such that (10.89) is ful¬lled, we compute the integral

a

I1 (f ) = 0 f (x)dx using for instance the adaptive Cavalieri-Simpson for-

mula introduced in Section 9.7.2 and denoted in the following by AD.

A natural choice of a bounding function for f is φ(x) = Kx3 e’x , for a

suitable constant K > 1. Thus, we have K ≥ ex /(ex ’ 1), for any x > 0,

that is, letting x = a, K = ea /(ea ’1). Substituting back into (10.89) yields

∞

a3 + 3a2 + 6a + 6 δ

f (x)dx ¤ = ·(a) ¤ .

ea ’ 1 2

a

464 10. Orthogonal Polynomials in Approximation Theory

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

’0.2

0 2 4 6 8 10 12 14 16

FIGURE 10.10. Distribution of quadrature nodes and graph of the integrand

function

Letting δ = 10’3 , we see that (10.89) is satis¬ed by taking a 16. Pro-

gram 77 for computing I1 (f ) with the AD method, setting hmin=10’3 and

tol=5 · 10’4 , yields the approximate value I1 6.4934 with a number of

(nonuniform) partitions equal to 25.

The distribution of the quadrature nodes produced by the adaptive algo-

rithm is plotted in Figure 10.10. Globally, using method 1. yields an approx-

imation of I(f ) equal to 6.4984. Table 10.1 shows, for sake of comparison,

some approximate values of I(f ) obtained using the Gauss-Laguerre formu-

lae with the number of nodes varying between 2 to 20. Notice that, taking

n = 4 nodes, the accuracy of the two computational procedures is roughly

equivalent.

n In (f )

2 6.413727469517582

3 6.481130171540022

4 6.494535639802632

5 6.494313365790864

10 6.493939967652101

15 6.493939402671590

20 6.493939402219742

∞

x3 /(ex ’ 1)dx with

TABLE 10.1. Approximate evaluation of I(f ) = 0

Gauss-Laguerre quadratures

10.13.2 Numerical Solution of Schr¨dinger Equation

o

Let us consider the following di¬erential equation arising in quantum me-

chanics known as the Schr¨dinger equation

o

‚2ψ

‚ψ

=’ x∈R

i , t > 0. (10.90)

2m ‚x2

‚t

The symbols i and denote the imaginary unit and the reduced Planck

constant, respectively. The complex-valued function ψ = ψ(x, t), the solu-

10.13 Applications 465

tion of (10.90), is called a wave function and the quantity |ψ(x, t)|2 de¬nes

the probability density in the space x of a free electron of mass m at time

t (see [FRL55]).

The corresponding Cauchy problem may represent a physical model for

describing the motion of an electron in a cell of an in¬nite lattice (for more

details see, e.g., [AF83]).

Consider the initial condition ψ(x, 0) = w(x), where w is the step func-

√

tion that takes the value 1/ 2b for |x| ¤ b and is zero for |x| > b, with

b = a/5, and where 2a represents the inter-ionic distance in the lattice.

Therefore, we are searching for periodic solutions, with period equal to 2a.

Solving problem (10.90) can be carried out using Fourier analysis as

follows. We ¬rst write the Fourier series of w and ψ (for any t > 0)

a

N/2’1

1

w(x)e’iπkx/a dx,

wk eiπkx/a ,

w(x) = wk =

2a

k=’N/2 ’a

(10.91)

a

N/2’1

1

ψ(x, t)e’iπkx/a dx.

ψk (t)eiπkx/a , ψk (t) =

ψ(x, t) =

2a

k=’N/2 ’a

Then, we substitute back (10.91) into (10.90), obtaining the following

Cauchy problem for the Fourier coe¬cients ψk , for k = ’N/2, . . . , N/2 ’ 1

± 2

ψ (t) = ’i kπ

ψk (t),

k

2m a (10.92)

ψk (0) = wk .

The coe¬cients {wk } have been computed by regularizing the coe¬cients

{wk } of the step function w using the Lanczos smoothing (10.56) in order

to avoid the Gibbs phenomenon arising around the discontinuities of w (see

Section 10.9.1).

After solving (10.92), we ¬nally get, recalling (10.91), the following expres-

sion for the wave function

N/2’1

wk e’iEk t/ eiπkx/a ,

ψN (x, t) = (10.93)

k=’N/2

where the coe¬cients Ek = (k 2 π 2 2 )/(2ma2 ) represent, from the physical

standpoint, the energy levels that the electron may assume in its motion

within the potential well.

To compute the coe¬cients wk (and, as a consequence, wk ), we have used

the MATLAB intrinsic function fft (see Section 10.9.2), employing N =

—¦

26 = 64 points and letting a = 10 A= 10’9 [m]. Time analysis has been

carried out up to T = 10 [s], with time steps of 1 [s]; in all the reported

466 10. Orthogonal Polynomials in Approximation Theory

0.35 0.35

0.3 0.3

0.25 0.25

0.2 0.2

0.15 0.15

0.1 0.1

0.05 0.05

0 0

’10 ’5 0 5 10 ’10 ’5 0 5 10

FIGURE 10.11. Probability density |ψ(x, t)|2 at t = 0, 2, 5 [s], corresponding to

a step function as initial datum: solution without ¬ltering (left), with Lanczos

¬ltering (right)

—¦

graphs, the x-axis is measured in [A], while the y-axes are respectively in

units of 105 [m’1/2 ] and 1010 [m’1 ].

In Figure 10.11 we draw the probability density |ψ(x, t)|2 at t = 0, 2

and 5 [s]. The result obtained without the regularizing procedure above is

shown on the left, while the same calculation with the “¬ltering” of the

Fourier coe¬cients is reported on the right. The second plot demonstrates

the smoothing e¬ect on the solution by the regularization, at the price of

a slight enlargement of the step-like initial probability distribution.

Finally, it is interesting to apply Fourier analysis to solve problem (10.90)

starting from a smooth initial datum. For this, we choose an initial probabil-

ity density w of Gaussian form such that w 2 = 1. The solution |ψ(x, t)|2 ,

this time computed without regularization, is shown in Figure 10.12, at

t = 0, 2, 5, 7, 9[s]. Notice the absence of spurious oscillations with respect

to the previous case.

0.2

0.18

0.16

0.14

0.12

0.1

0.08

0.06

0.04

0.02

0

’10 ’5 0 5 10

FIGURE 10.12. Probability density |ψ(x, t)|2 at t = 0, 2, 5, 7, 9[s], corresponding

to an initial datum with Gaussian form

10.14 Exercises 467

10.14 Exercises

1. Prove the three-term relation (10.11).

[Hint: set x = cos(θ), for 0 ¤ θ ¤ π.]

2. Prove (10.31).

[Hint: ¬rst prove that vn n = (vn , vn )1/2 , Tk n = Tk w for k < n and

Tn 2 = 2 Tn 2 (see [QV94], formula (4.3.16)). Then, the thesis follows

n w

from (10.29) multiplying by Tl (l = k) and taking (·, ·)n .]

3. Prove (10.24) after showing that (f ’ ΠGL f ) ¤ Cn1’s f s,ω .

ω

n

[Hint: use the Gagliardo-Nirenberg inequality

max |f (x)| ¤ f 1/2 1/2

f

’1¤x¤1

valid for any f ∈ L2 with f ∈ L2 . Next, use the relation that has been just

shown to prove (10.24).]

1/2

is a norm for Pn .

4. Prove that the discrete seminorm f = (f, f )n

n

5. Compute weights and nodes of the following quadrature formulae

b n

w(x)f (x)dx = ωi f (xi ),

i=0

a

in such a way that the order is maximum, setting

√

ω(x) = x, a = 0, b = 1, n = 1;

a = ’1,

ω(x) = 2x2 + 1, b = 1, n = 0;

2 if 0 < x ¤ 1,

a = ’1,

ω(x) = b = 1, n = 1.

1 if ’ 1 ¤ x ¤ 0

√

[Solution: for ω(x) = x, the nodes x1 = 5 + 2 10/7, x2 = 5 ’ 2 10/7

9 9 9 9

are obtained, from which the weights can be computed (order 3); for ω(x) =

2x2 + 1, we get x1√ 3/5 and ω1 = 5/3 (order 1); for ω(x) = 2x2 + 1, we

= √

have x1 = 22 + 22 155, x2 = 22 ’ 22 155 (order 3).]

1 1 1 1

6. Prove (10.40).

—

[Hint: notice that (ΠGL f, Lj )n = fk (Lk , Lj )n = . . . , distinguishing the

n k

case j < n from the case j = n.]

7. Show that ||| · |||, de¬ned in (10.45), is an essentially strict seminorm.

[Solution : use the Cauchy-Schwarz inequality (1.14) to check that the

triangular inequality is satis¬ed. This proves that ||| · ||| is a seminorm. The

second part of the exercise follows after a direct computation.]

8. Consider in an interval [a, b] the nodes

b’a

1

xj = a + j ’ j = 1, 2, . . . , m

2 m

468 10. Orthogonal Polynomials in Approximation Theory

for m ≥ 1. They are the midpoints of m equally spaced intervals in [a, b].

Let f be a given function; prove that the least-squares polynomial rn with

respect to the weight w(x) = 1 minimizes the error average, de¬ned as

1/2

m

1

[f (xj ) ’ rn (xj )]2

E = lim .

m j=1

m’∞

9. Consider the function

1 2

n

f (x) ’ j

F (a0 , a1 , . . . , an ) = aj x dx

j=0

0

and determine the coe¬cients a0 , a1 , . . . , an in such a way that F is mini-

mized. Which kind of linear system is obtained?

[Hint: enforce the conditions ‚F/‚ai = 0 with i = 0, 1, . . . , n. The matrix

of the ¬nal linear system is the Hilbert matrix (see Example 3.2, Chapter

3) which is strongly ill-conditioned.]

11

Numerical Solution of Ordinary

Di¬erential Equations

In this chapter we deal with the numerical solutions of the Cauchy problem

for ordinary di¬erential equations (henceforth abbreviated by ODEs). After

a brief review of basic notions about ODEs, we introduce the most widely

used techniques for the numerical approximation of scalar equations. The

concepts of consistency, convergence, zero-stability and absolute stability

will be addressed. Then, we extend our analysis to systems of ODEs, with

emphasis on sti¬ problems.

11.1 The Cauchy Problem

The Cauchy problem (also known as the initial-value problem) consists of

¬nding the solution of an ODE, in the scalar or vector case, given suitable

initial conditions. In particular, in the scalar case, denoting by I an interval

of R containing the point t0 , the Cauchy problem associated with a ¬rst

order ODE reads:

¬nd a real-valued function y ∈ C 1 (I), such that

t ∈ I,

y (t) = f (t, y(t)),

(11.1)

y(t0 ) = y0 ,

where f (t, y) is a given real-valued function in the strip S = I —(’∞, +∞),

which is continuous with respect to both variables. Should f depend on t

only through y, the di¬erential equation is called autonomous.

470 11. Numerical Solution of Ordinary Di¬erential Equations

Most of our analysis will be concerned with one single di¬erential equa-

tion (scalar case). The extension to the case of systems of ¬rst-order ODEs

will be addressed in Section 11.9.

If f is continuous with respect to t, then the solution to (11.1) satis¬es

t

y(t) ’ y0 = f („, y(„ ))d„. (11.2)

t0

Conversely, if y is de¬ned by (11.2), then it is continuous in I and y(t0 ) =

y0 . Moreover, since y is a primitive of the continuous function f (·, y(·)),

y ∈ C 1 (I) and satis¬es the di¬erential equation y (t) = f (t, y(t)).

Thus, if f is continuous the Cauchy problem (11.1) is equivalent to the

integral equation (11.2). We shall see later on how to take advantage of

this equivalence in the numerical methods.

Let us now recall two existence and uniqueness results for (11.1).

1. Local existence and uniqueness.

Suppose that f (t, y) is locally Lipschitz continuous at (t0 , y0 ) with

respect to y, that is, there exist two neighborhoods, J ⊆ I of t0 of

width rJ , and Σ of y0 of width rΣ , and a constant L > 0, such that

|f (t, y1 ) ’ f (t, y2 )| ¤ L|y1 ’ y2 | ∀t ∈ J, ∀y1 , y2 ∈ Σ. (11.3)

Then, the Cauchy problem (11.1) admits a unique solution in a neigh-

borhood of t0 with radius r0 with 0 < r0 < min(rJ , rΣ /M, 1/L),

where M is the maximum of |f (t, y)| on J — Σ. This solution is called

the local solution.

Notice that condition (11.3) is automatically satis¬ed if f has con-

tinuous derivative with respect to y: indeed, in such a case it su¬ces

to choose L as the maximum of |‚f (t, y)/‚y| in J — Σ.

2. Global existence and uniqueness. The problem admits a unique