<< . .

. 67
( : 95)



. . >>

E= dx,
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

<< . .

. 67
( : 95)



. . >>