<<

. 23
( 26)



>>

® 
xi+1 xi+1
xi

µ °ui’1 •i •i+1 dx»
2
•i’1 •i dx + ui (•i ) dx + ui+1
® 
xi’1 xi’1 xi
xi+1 xi+1
xi

+β °ui’1 •i+1 •i dx» = 0.
•i’1 •i dx + ui •i •i dx + ui+1
xi’1 xi’1 xi

Assuming a uniform partition of [0, 1] with xi = xi’1 + h for i = 1, . . . , n,
h = 1/n, and noting that •j (x) = h if xj’1 ¤ x ¤ xj , •j (x) = ’ h if
1 1

xj ¤ x ¤ xj+1 , we deduce that
µ 1
(’ui’1 + 2ui ’ ui+1 ) + β (ui+1 ’ ui’1 ) = 0, i = 1, . . . , n ’ 1.
h 2
(12.75)
12.5 Advection-Di¬usion Equations 563

Multiplying by h/µ and de¬ning the local P´clet number to be
e

|β|h
Pe =

we ¬nally obtain

(Pe ’ 1) ui+1 + 2ui ’ (Pe + 1) ui’1 = 0. i = 1, . . . , n ’ 1. (12.76)

This is a linear di¬erence equation which admits a solution of the form
ui = A1 ρi + A2 ρi for suitable constants A1 , A2 (see Section 11.4), where
1 2
ρ1 and ρ2 are the two roots of the following characteristic equation

(Pe ’ 1) ρ2 + 2ρ ’ (Pe + 1) = 0.

±
Thus
√ 1 + Pe
’1 ± 1 + Pe2 ’ 1  ,
1 ’ Pe
ρ1,2 = =
 1.
Pe ’ 1

Imposing the boundary conditions at x = 0 and x = 1 gives
n
A1 = 1/(1 ’ ), A2 = ’A1 ,
1+Pe
1’Pe

so that the solution of (12.76) is
i n
1 + Pe 1 + Pe
1’ / 1’
ui = i = 0, . . . , n.
1 ’ Pe 1 ’ Pe

We notice that if Pe > 1, a power with a negative base appears at the
numerator which gives rise to an oscillating solution. This is clearly visible
in Figure 12.7 where for several values of the local P´clet number, the
e
solution of (12.76) is compared with the exact solution (sampled at the
mesh nodes) corresponding to a value of the global P´clet equal to 50.
e
The simplest remedy for preventing the oscillations consists of choosing
a su¬ciently small mesh stepsize h in such a way that Pe < 1. However this
approach is often impractical: for example, if β = 1 and µ = 5 · 10’5 one
should take h < 10’4 which amounts to dividing [0, 1] into 10000 subin-
tervals, a strategy that becomes unfeasible when dealing with multidimen-
sional problems. Other strategies can be pursued, as will be addressed in
the next sections.


12.5.2 The Relationship between Finite Elements and Finite
Di¬erences; the Numerical Viscosity
To examine the behaviour of the ¬nite di¬erence method (FD) when applied
to the solution of advection-di¬usion problems and its relationship with the
564 12. Two-Point Boundary Value Problems

1




0.5




0




’0.5


0.75 0.8 0.85 0.9 0.95 1



FIGURE 12.7. Finite di¬erence solution of the advection-di¬usion problem
(12.70) (with Pegl = 50) for several values of the local P´clet number. Solid
e
line: exact solution, dot-dashed line: Pe = 2.63, dotted line: Pe = 1.28, dashed
line: Pe = 0.63

¬nite element method (FE), we again consider the one-dimensional problem
(12.70) with a uniform meshsize h.
To ensure that the local discretization error is of second order we approx-
imate u (xi ) and u (xi ), i = 1, . . . , n ’ 1, by the centred ¬nite di¬erences
(10.61) and (10.65) respectively (see Section 10.10.1). We obtain the fol-
lowing FD problem
±
 ’µ ui+1 ’ 2ui + ui’1 + β ui+1 ’ ui’1 = 0, i = 1, . . . , n ’ 1
h2 (12.77)
2h
 u = 0, u = 1.
0 n

If we multiply by h for every i = 1, . . . , n ’ 1, we obtain exactly the same
equation (12.75) that was obtained using piecewise linear ¬nite elements.
The equivalence between FD and FE can be pro¬tably employed to de-
vise a cure for the oscillations arising in the approximate solution of (12.75)
when the local P´clet number is larger than 1. The important observation
e
here is that the instability in the FD solution is due to the fact that the
discretization scheme is a centred one. A possible remedy consists of ap-
proximating the ¬rst derivative by a one-sided ¬nite di¬erence according
to the direction of the transport ¬eld. Precisely, we use the backward dif-
ference if the convective coe¬cient β is positive and the forward di¬erence
otherwise. The resulting scheme when β > 0 is
ui+1 ’ 2ui + ui’1 ui ’ ui’1
’µ i = 1, . . . , n ’ 1
+β =0 (12.78)
h2 h
which, for µ = 0 reduces to ui = ui’1 and therefore yields the desired
constant solution of the limit problem βu = 0. This one-side discretization
12.5 Advection-Di¬usion Equations 565

of the ¬rst derivative is called upwind di¬erencing: the price to be paid for
the enhanced stability is a loss of accuracy since the upwind ¬nite di¬erence
introduces a local discretization error of O(h) and not of O(h2 ) as happens
using centred ¬nite di¬erences.
Noting that
ui ’ ui’1 ui+1 ’ ui’1 h ui+1 ’ 2ui + ui’1

= ,
h2
h 2h 2
the upwind ¬nite di¬erence can be interpreted as the sum of a centred ¬nite
di¬erence approximating the ¬rst derivative and of a term proportional to
the discretization of the second-order derivative. Consequently, (12.78) is
equivalent to
ui+1 ’ 2ui + ui’1 ui+1 ’ ui’1
’µh i = 1, . . . , n ’ 1 (12.79)
+β =0
h2 2h
where µh = µ(1 + Pe). This amounts to having replaced the di¬erential
equation (12.72) with the perturbed one

’µh u + βu = 0, (12.80)

then using centred ¬nite di¬erences to approximate both u and u . The
perturbation
βh
’µ Pe u = ’ u (12.81)
2
is called the numerical viscosity (or arti¬cial di¬usion). In Figure 12.8
a comparison between centred and upwinded discretizations of problem
(12.72) is shown.
More generally, we can resort to a centred scheme of the form (12.80)
with the following viscosity

µh = µ(1 + φ(Pe)) (12.82)

where φ is a suitable function of the local P´clet number satisfying
e

lim φ(t) = 0.
t’0+

Notice that when φ(t) = 0 for all t, one recovers the centred ¬nite di¬erence
method (12.77), while if φ(t) = t the upwind ¬nite di¬erence scheme (12.78)
(or, equivalently, (12.79)) is obtained. Other choices are admissible as well.
For instance, taking

φ(t) = t ’ 1 + B(2t) (12.83)

where B(t) is the Bernoulli function de¬ned as B(t) = t/(et ’ 1) for t = 0
and B(0) = 1, yields the so called exponential ¬tting ¬nite di¬erence scheme
which is also well known as the Scharfetter-Gummel (SG) method [SG69].
566 12. Two-Point Boundary Value Problems

1.2



1



0.8



0.6



0.4



0.2



0



’0.2



’0.4

0.7 0.75 0.8 0.85 0.9 0.95 1



FIGURE 12.8. Finite di¬erence solution of (12.72) (with µ = 1/100) using a
centred discretization (dashed and dotted lines) and the arti¬cial viscosity (12.81)
(dot-dashed and starred lines). The solid line denotes the exact solution. Notice
the e¬ect of eliminating the oscilations when the local P´clet number is large;
e
conversely, notice also the corresponding loss of accuracy for low values of the
local P´clet number
e

2.5



2



1.5



1



0.5



0
0 0.5 1 1.5 2 2.5




FIGURE 12.9. The functions φU P (dash-dotted line) and φSG (solid line)


Remark 12.6 Denoting respectively by φC , φU P and φSG the previous
three functions, i.e. φC = 0, φU P (t) = t and φSG (t) = t ’ 1 + B(2t) (see
Figure 12.9), we notice that φSG φU P as Pe ’ +∞ while φSG = O(h2 )
and φU P = O(h) if Pe ’ 0+ . Therefore, the SG method is second-order
accurate with respect to h and for this reason it is an optimal viscosity
upwind method. Actually, one can show (see [HGR96], pp. 44-45) that
if f is piecewise constant over the grid partition the SG scheme yields a
numerical solution uSG which is nodally exact, i.e., uSG (xi ) = u(xi ) for
h h
each node xi , irrespectively of h (and, thus, of the size of the local P´clet
e
number). This is demonstrated in Figure 12.10.
12.5 Advection-Di¬usion Equations 567

1




0.5




0

0.9 0.95 1


FIGURE 12.10. Comparison between the numerical solutions of problem (12.72)
(with µ = 1/200) obtained using the arti¬cial viscosity (12.81) (dashed line where
the symbol denotes the nodal values) and with the optimal viscosity (12.83)
(dotted line where the symbol —¦ denotes the nodal values) in the case where
Pe = 1.25. The solid line denotes the exact solution

The new local P´clet number associated with the scheme (12.79)-(12.82) is
e
de¬ned as
|β|h Pe
Pe— = = .
2µh (1 + φ(Pe))
For both the upwind and the SG schemes we have Pe— < 1 for any value
of h. This implies that the matrix associated with these methods is a M-
matrix for any h (see De¬nition 1.25 and Exercise 13), and, in turn, that the
numerical solution uh satis¬es a discrete maximum principle (see Section
12.2.2).


12.5.3 Stabilized Finite Element Methods
In this section we extend the use of numerical viscosity introduced in
the previous section for ¬nite di¬erences to the Galerkin method using
¬nite elements of arbitrary degree k ≥ 1. For this purpose we consider
the advection-di¬usion problem (12.70) where the viscosity coe¬cient µ is
replaced by (12.82). This yields the following modi¬cation of the original
Galerkin problem (12.72):
0 k,0
¬nd uh ∈ Xh = vh ∈ Xh : vh (0) = vh (1) = 0
k
such that

1
0 k,0
= ’ βvh dx ∀vh ∈ Xh ,
ah (uh , vh ) (12.84)
0

where
ah (u, v) = a(u, v) + b(u, v),
568 12. Two-Point Boundary Value Problems

and
1

b(u, v) = µ φ(Pe) u v dx
0
is called the stabilization term. Since ah (v, v) = µh |v|2 for all v ∈ H1 (0, 1)
1 0
and µh /µ = (1 + φ(Pe)) ≥ 1, the modi¬ed problem (12.84) enjoys more
favorable monotonicity properties than the corresponding non stabilized
Galerkin formulation (12.75).
0
0
To prove convergence, it is su¬cient to show that uh tends to u as h ’ 0,
0
where u (x) = u(x) ’ x. This is done in the following theorem, where we
0
assume that u (and henceforth u) has the required regularity.

Theorem 12.4 If k = 1 then
0 0
0
| u ’ uh |H1 (0,1) ¤ Ch G(u) (12.85)
0
where C > 0 is a suitable constant independent of h and u, and
±
 |u| 1
0 0
 H (0,1) + | u |H2 (0,1) for the upwind method
0
u) =
G(
0
 |u| 2 for the SG method.
H (0,1)

Further, if k = 2 the SG method yields the improved error estimate
0 0 0
0
| u ’ uh |H1 (0,1) ¤ Ch2 (| u |H1 (0,1) + | u |H3 (0,1) ). (12.86)
Proof. From (12.70) we obtain
1
0
a(u, vh ) = ’ ∀vh ∈ Xh .
k,0
βvh dx,
0

By comparison with (12.84) we get
0 0
0
ah (u ’ uh , vh ) = b(u, vh ), ∀vh ∈ Xh .
k,0
(12.87)
0 0
Denote by Eh =u ’ uh the discretization error and recall that the space H1 (0, 1)
0
is endowed with the norm (12.49). Then,
0 0 0 0
µh |Eh |2 1 (0,1) = ah (Eh , Eh ) = ah (Eh , u ’Πk u) + ah (Eh , Πk u ’ uh )
h h
H
0 0 0 0 0
= ah (Eh , u ’Πk u) + b(u, Πk u ’ uh )
h h

0 0
where we have applied (12.87) with vh = Πk u ’ uh . Using the Cauchy-Schwarz
h
inequality we get
0 0
µh |Eh |2 1 (0,1) ¤ Mh |Eh |H1 (0,1) | u ’Πk u |H1 (0,1)
h
H
1
(12.88)
0
0 0

(Πk u
+µφ(Pe) u uh ) dx
h
0
12.5 Advection-Di¬usion Equations 569

where Mh = µh + |β|CP is the continuity constant of the bilinear form ah (·, ·)
and CP is the Poincar´ constant introduced in (12.16).
e
Notice that if k = 1 (corresponding to piecewise linear ¬nite elements) and
φ = φSG (SG optimal viscosity) the quantity in the second integral is identically
0
0
zero since uh = Π1 u, as pointed out in Remark 12.6. Then, from (12.88) we get
h

|β|CP 0 0
|Eh |H1 (0,1) ¤ | u ’Π1 u |H1 (0,1) .
1+ h
µh

Noting that µh > µ, using (12.71) and the interpolation estimate (8.27), we ¬nally
obtain the error bound
0
|Eh |H1 (0,1) ¤ C(1 + 2Pegl CP )h| u |H2 (0,1) .

In the general case the error inequality (12.88) can be further manipulated. Using
the Cauchy-Schwarz and triangular inequalities we obtain
1
0
0 0 0 0
0
u (Πk u ’ uh ) dx ¤ | u |H1 (0,1) (|Πk u ’ u |H1 (0,1) + |Eh |H1 (0,1) )
h h
0

from which
0 0
µh |Eh |2 1 (0,1) ¤ |Eh |H1 (0,1) Mh | u ’Πk u |H1 (0,1)
h
H
0 0 0 0
+µφ(Pe)| u |H1 (0,1) + µφ(Pe)| u |H1 (0,1) | u ’Πk u |H1 (0,1) .
h


Using again the interpolation estimate (8.27) yields
0 0
µh |Eh |2 1 (0,1) ¤ |Eh |H1 (0,1) Mh Chk | u |Hk+1 (0,1) + µφ(Pe)| u |H1 (0,1)
H
0 0
+Cµφ(Pe)| u |H1 (0,1) hk | u |Hk+1 (0,1) .

Using Young™s inequality (12.40) gives

µh |Eh |2 1 (0,1)
H
µh |Eh |2 1 (0,1) ¤
H
2
3 0 0
(Mh Chk | u |Hk+1 (0,1) )2 + (µφ(Pe)| u |H1 (0,1) )2
+
4µh
from which it follows that
2 2
3 Mh 3µ
0 0
|Eh |2 1 (0,1) ¤ C 2 h2k | u |2 k+1 (0,1) + φ(Pe)2 | u |2 1 (0,1)
H
H
H
2 µh 2 µh
2µ 0 k0
+ φ(Pe)| u |H1 (0,1) Ch | u |Hk+1 (0,1) .
µh
Again using the fact that µh > µ and the de¬nition (12.71) we get (Mh /µh ) ¤
(1 + 2CP Pegl ) and then
32 0
|Eh |2 1 (0,1) ¤ C (1 + 2CP Pegl )2 h2k | u |2 k+1 (0,1)
H
H
2
3
0 0 0
+2φ(Pe)Chk | u |H1 (0,1) | u |Hk+1 (0,1) + φ(Pe)2 | u |2 1 (0,1) ,
H
2
570 12. Two-Point Boundary Value Problems

which can be bounded further as
0
|Eh |2 1 (0,1) ¤ M h2k | u |2 k+1 (0,1)
H H
(12.89)
0 0 0
φ(Pe)hk | u |H1 (0,1) | u |Hk+1 (0,1) + φ(Pe)2 | u |2 1 (0,1)
+ H


for a suitable positive constant M.
If φU P = Cµ h, where Cµ = β/µ, we obtain
0
|Eh |2 1 (0,1) ¤ Ch2 h2k’2 | u |2 k+1 (0,1)
H H
0 0 0
+hk’1 | u |H1 (0,1) | u |Hk+1 (0,1) + | u |2 1 (0,1)
H


which shows that using piecewise linear ¬nite elements (i.e., k = 1) plus the
upwind arti¬cial viscosity gives the linear convergence estimate (12.85).
In the case φ = φSG , assuming that for h su¬ciently small φSG ¤ Kh2 , for a
suitable positive constant K, we get
0
|Eh |2 1 (0,1) ¤ Ch4 h2(k’2) | u |2 k+1 (0,1)
H H
0 0 0
+hk’2 | u |H1 (0,1) | u |Hk+1 (0,1) + | u |2 1 (0,1)
H


which shows that using quadratic ¬nite elements (i.e., k = 2) plus the optimal
3
arti¬cial viscosity gives the second-order convergence estimate (12.86).

Programs 97 and 98 implement the computation of the arti¬cial and opti-
mal arti¬cial viscosities (12.81) and (12.83), respectively. These viscosities
can be selected by the user setting the input parameter stabfun in Pro-
gram 94 equal to artvisc or sgvisc. The function sgvisc employs the
function bern to evaluate the Bernoulli function in (12.83).

Program 97 - artvisc : Arti¬cial viscosity
function [Kupw,rhsbc] = artvisc(Nx, h, nu, beta)
Peclet=0.5*h*abs(beta);
for i=2:Nx,
dd(i-1)=(Peclet(i-1)+Peclet(i))/h;
if i > 2, ld(i-2)=-Peclet(i-1)/h; end
if i < Nx, ud(i-1)=-Peclet(i)/h; end
end
Kupw=spdiags([[ld 0]™,dd™,[0 ud]™],-1:1,Nx-1,Nx-1);
rhsbc = - [Peclet(1)/h, Peclet(Nx)/h];


Program 98 - sgvisc : Optimal arti¬cial viscosity
function [Ksg,rhsbc] = sgvisc(Nx, h, nu, beta)
Peclet=0.5*h*abs(beta)./nu;
[bp, bn]=bern(2*Peclet);
Peclet=Peclet-1+bp;
12.5 Advection-Di¬usion Equations 571

for i=2:Nx,
dd(i-1)=(nu(i-1)*Peclet(i-1)+nu(i)*Peclet(i))/h;
if i > 2, ld(i-2)=-nu(i-1)*Peclet(i-1)/h; end
if i < Nx, ud(i-1)=-nu(i)*Peclet(i)/h; end
end
Ksg=spdiags([[ld 0]™,dd™,[0 ud]™],-1:1,Nx-1,Nx-1);
rhsbc = - [nu(1)*Peclet(1)/h, nu(Nx)*Peclet(Nx)/h];



Program 99 - bern : Evaluation of the Bernoulli function

function [bp,bn]=bern(x)
xlim=1e-2; ax=abs(x);
if (ax == 0), bp=1.; bn=1.; return; end;
if (ax > 80),
if (x > 0), bp=0.; bn=x; return;
else, bp=-x; bn=0.; return; end;
end;
if (ax > xlim),
bp=x/(exp(x)-1); bn=x+bp; return;
else
ii=1; fp=1.;fn=1.; df=1.; s=1.;
while (abs(df) > eps),
ii=ii+1; s=-s; df=df*x/ii;
fp=fp+df; fn=fn+s*df;
bp=1./fp; bn=1./fn;
end;
return;
end




Example 12.3 We use Program 94 supplied with Programs 97 and 98 for the
numerical approximation of problem (12.70) in the case µ = 10’2 . Figure 12.11
shows the convergence behavior as a function of log(h) of the Galerkin method
without (G) and with numerical viscosity (upwind (UP) and SG methods are
employed). The ¬gure shows the logarithmic plots of the L2 (0, 1) and H1 (0, 1)-
norms, where the solid line denotes the UP method and the dashed and dotted
lines denote the G and SG methods, respectively. It is interesting to notice that
the UP and SG schemes exhibit the same (linear) convergence rate as the pure
Galerkin method in the H1 -norm, while the accuracy of the UP scheme in the
L2 -norm deteriorates dramatically because of the e¬ect of the arti¬cial viscosity
which is O(h). Conversely, the SG converges quadratically since the introduced
numerical viscosity is in this case O(h2 ) as h tends to zero. •
572 12. Two-Point Boundary Value Problems

102


101


100
||u ’ u h || 1
H (0,1)
’1
10


10’2

||u ’ u h || 2
L (0,1)
10’3


10’4 ’3
10’2 10’1
10


FIGURE 12.11. Convergence analysis for an advection-di¬usion problem


12.6 A Quick Glance to the Two-Dimensional Case
The game that we want to play is to extend (in a few pages) the basic ideas
illustrated so far to the two-dimensional case. The obvious generalization of
problem (12.1)-(12.2) is the celebrated Poisson problem with homogeneous
Dirichlet boundary condition

’ u=f in „¦,
(12.90)
u=0 on ‚„¦,

where u = ‚ 2 u/‚x2 + ‚ 2 u/‚y 2 is the Laplace operator and „¦ is a two-
dimensional bounded domain whose boundary is ‚„¦. If we allow „¦ to be
the unit square „¦ = (0, 1)2 , the ¬nite di¬erence approximation of (12.90)
that mimics (12.10) is

Lh uh (xi,j ) = f (xi,j ) for i, j = 1, . . . , N ’ 1,
(12.91)
uh (xi,j ) = 0 if i = 0 or N, j = 0 or N,

where xi,j = (ih, jh) (h = 1/N > 0) are the grid points and uh is a grid
function. Finally, Lh denotes any consistent approximation of the operator
L = ’ . The classical choice is

1
(4ui,j ’ ui+1,j ’ ui’1,j ’ ui,j+1 ’ ui,j’1 ) ,
Lh uh (xi,j ) = (12.92)
h2

where ui,j = uh (xi,j ), which amounts to adopting the second-order centred
discretization of the second derivative (10.65) in both x and y directions
(see Figure 12.12, left).
The resulting scheme is known as the ¬ve-point discretization of the
Laplacian. It is an easy (and useful) exercise for the reader to check that
the associated matrix Afd has (N ’ 1)2 rows and is pentadiagonal, with the
12.6 A Quick Glance to the Two-Dimensional Case 573




xi+1,j
xij xi,j+1
xi,j’1

xi’1,j



FIGURE 12.12. Left: ¬nite di¬erence grid and stencil for a squared domain.
Right: ¬nite element triangulation of a region around a hole


i-th row given by
±
 4 if j = i,

1
’1 if j = i ’ N ’ 1, i ’ 1, i + 1, i + N + 1,
(afd )ij = 2 (12.93)
h 

0 otherwise.

Moreover, Afd is symmetric positive de¬nite and it is also an M-matrix (see
Exercise 14). As expected, the consistency error associated with (12.92) is
second-order with respect to h, and the same holds for the discretization
error u ’ uh h,∞ of the method. More general boundary conditions than
the one considered in (12.90) can be dealt with by properly extending to
the two-dimensional case the mirror imaging technique described in Section
12.2.3 and in Exercise 10 (for a thorough discussion of this subject, see, e.g.,
[Smi85]).

The extension of the Galerkin method is (formally speaking) even more
straightforward and actually is still readable as (12.43) with, however, the
implicit understanding that both the function space Vh and the bilinear
form a(·, ·) be adapted to the problem at hand. The ¬nite element method
corresponds to taking

Vh = vh ∈ C 0 („¦) : vh |T ∈ Pk (T ) ∀T ∈ Th , vh |‚„¦ = 0 , (12.94)

where Th denotes here a triangulation of the domain „¦ as previously intro-
duced in Section 8.5.2, while Pk (k ≥ 1) is the space of piecewise polyno-
mials de¬ned in (8.35). Note that „¦ needs not to be a rectangular domain
(see Figure 12.12, right).
574 12. Two-Point Boundary Value Problems

As of the bilinear form a(·, ·), the same kind of mathematical manipula-
tions performed in Section 12.4.1 lead to

∇uh · ∇vh dxdy,
a(uh , vh ) =
„¦

where we have used the following Green™s formula that generalizes the
formula of integration by parts

’ ∇u · ∇v dxdy ’ ∇u · n v dγ,
u v dxdy = (12.95)
„¦ „¦ ‚„¦

for any u, v smooth enough and where n is the outward normal unit vector
on ‚„¦ (see Exercise 15).
The error analysis for the two-dimensional ¬nite element approximation
of (12.90) can still be performed through the combined use of Ce`™s lemma
a
and interpolation error estimates as in Section 12.4.5 and is summarized in
the following result, which is the two-dimensional counterpart of Property
12.1 (for its proof we refer, e.g., to [QV94], Theorem 6.2.1).

Property 12.2 Let u ∈ H1 („¦) be the exact solution of (12.90) and uh ∈ Vh
0
be its ¬nite element approximation using continuous piecewise polynomials
of degree k ≥ 1. Assume also that u ∈ Hs („¦) for some s ≥ 2. Then the
following error estimate holds
M
u ’ uh ¤ Chl u (12.96)
H1 („¦) Hl+1 („¦)
0 ±0
where l = min(k, s ’ 1). Under the same assumptions, one can also prove
that

u ’ uh ¤ Chl+1 u Hl+1 („¦) . (12.97)
L2 („¦)

We notice that, for any integer s ≥ 0, the Sobolev space Hs („¦) introduced
above is de¬ned as the space of functions with the ¬rst s partial derivatives
(in the distributional sense) belonging to L2 („¦). Moreover, H1 („¦) is the
0
1
space of functions of H („¦) such that u = 0 on ‚„¦. The precise mathemat-
ical meaning of this latter statement has to be carefully addressed since,
for instance, a function belonging to H1 („¦) does not necessarily mean to
0
be continuous everywhere. For a comprehensive presentation and analysis
of Sobolev spaces we refer to [Ada75] and [LM68].
Following the same procedure as in Section 12.4.3, we can write the ¬nite
element solution uh as
N
uh (x, y) = uj •j (x, y)
j=1
12.7 Applications 575

N
where {•j }j=1 is a basis for Vh . An example of such a basis in the case
k = 1 is provided by the so-called hat functions introduced in Section 8.5.2
(see Figure 8.7, right). The Galerkin ¬nite element method leads to the
solution of the linear system Afe u = f , where (afe )ij = a(•j , •i ).
Exactly as happens in the one-dimensional case, the matrix Afe is sym-
metric positive de¬nite and, in general, sparse, the sparsity pattern being
strongly dependent on the topology of Th and the numbering of its nodes.
Moreover, the spectral condition number of Afe is still O(h’2 ), which im-
plies that solving iteratively the linear system demands for the use of a
suitable preconditioner (see Section 4.3.2). If, instead, a direct solver is
used, one should resort to a suitable renumbering procedure, as explained
in Section 3.9.



12.7 Applications
In this section we employ the ¬nite element method for the numerical
approximation of two problems arising in ¬‚uid mechanics and in the sim-
ulation of biological systems.


12.7.1 Lubrication of a Slider
Let us consider a rigid slider moving in the direction x along a physical
support from which it is separated by a thin layer of a viscous ¬‚uid (which
is the lubricant). Suppose that the slider, of length L, moves with a ve-
locity U with respect to the plane support which is supposed to have an
in¬nite length. The surface of the slider that is faced towards the support
is described by the function s (see Figure 12.13, left).
Denoting by µ the viscosity of the lubricant, the pressure p acting on the
slider can be modeled by the following Dirichlet problem
±
 s3

’ = ’(U s) x ∈ (0, L),
p


 p(0) = 0, p(L) = 0.

Assume in the numerical computations that we are working with a conver-
gent-divergent slider of unit length, whose surface is s(x) = 1’3/2x+9/8x2
with µ = 1.
Figure 12.13 (right) shows the solution obtained using linear and quadratic
¬nite elements with an uniform grid size h = 0.2. The linear system has
been solved by the nonpreconditioned conjugate gradient method. To re-
duce the Euclidean norm of the residual below 10’10 , 4 iterations are
needed in the case of linear ¬nite elements while 9 are required for quadratic
¬nite elements.
576 12. Two-Point Boundary Value Problems


L
0.09

U 0.08


0.07


0.06


0.05


0.04


s(x) 0.03


0.02



x 0.01


0


’0.01
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1




FIGURE 12.13. Left: geometrical parameters of the slider considered in Section
12.7.1; right: pressure on a converging-diverging slider. The solid line denotes
the solution obtained used linear ¬nite elements (the symbol —¦ denotes the nodal
values), while the dashed line denotes the solution obtained using quadratic ¬nite
elements

Table 12.2 reports a numerical study of the condition number K2 (Afe ) as
a function of h. In the case of linear ¬nite elements we have denoted the
matrix by A1 , while A2 is the corresponding matrix for quadratic elements.
Here we assume that the condition number approaches h’p as h tends to
zero; the numbers pi are the estimated values of p. As can be seen, in
both cases the condition number grows like h’2 , however, for every ¬xed
h, K2 (A2 ) is much bigger than K2 (A1 ).

h K2 (A1 ) p1 K2 (A2 ) p2
0.10000 63.951 “ 455.24
0.05000 348.21 2.444 2225.7 2.28
0.02500 1703.7 2.290 10173.3 2.19
0.01250 7744.6 2.184 44329.6 2.12
0.00625 33579 2.116 187195.2 2.07
TABLE 12.2. Condition number of the sti¬ness matrix for linear and quadratic
¬nite elements




12.7.2 Vertical Distribution of Spore Concentration over
Wide Regions
In this section we are concerned with di¬usion and transport processes of
spores in the air, such as the endspores of bacteria and the pollen of ¬‚ow-
ering plants. In particular, we study the vertical concentration distribution
of spores and pollen grains over a wide area. These spores, in addition to
settling under the in¬‚uence of gravity, di¬use passively in the atmosphere.
12.7 Applications 577

The basic model assumes the di¬usivity ν and the settling velocity β to
be given constants, the averaging procedure incorporating various physical
processes such as small-scale convection and horizontal advection-di¬usion
which can be neglected over a wide horizontal area. Denoting by x ≥ 0 the
vertical direction, the steady-state distribution u(x) of the spore concen-
tration is the solution of

’νu + βu = 0 0 < x < H,
(12.98)
’νu (H) + βu(H) = 0,
u(0) = u0 ,

where H is a ¬xed height at which we assume a vanishing Neumann condi-
tion for the advective-di¬usive ¬‚ux ’νu +βu (see Section 12.4.1). Realistic
values of the coe¬cients are ν = 10 m2 s’1 and β = ’0.03 ms’1 ; as for
u0 , a reference concentration of 1 pollen grain per m3 has been used in the
numerical experiments, while the height H has been set equal to 10 km.
The global P´clet number is therefore Pegl = 15.
e
The Galerkin ¬nite element method with piecewise linear ¬nite elements
has been used for the approximation of (12.98). Figure 12.14 (left) shows
the solution computed by running Program 94 on a uniform grid with step-
size h = H/10. The solution obtained using the (non stabilized) Galerkin
formulation (G) is denoted by the solid line, while the dash-dotted and
dashed lines refer to the Scharfetter-Gummel (SG) and upwind (UP) sta-
bilized methods. Spurious oscillations can be noticed in the G solution while
the one obtained using UP is clearly overdi¬used with respect to the SG
solution that is nodally exact. The local P´clet number is equal to 1.5 in
e
this case. Taking h = H/100 yields for the pure Galerkin scheme a stable
result as shown in Figure 12.14 (right) where the solutions furnished by G
(solid line) and UP (dashed line) are compared.

1
1.2
0.9

1 0.8

0.7
0.8

0.6
0.6
0.5

0.4 0.4

0.3
0.2
0.2
0
0.1

0
’0.2
0 200 400 600 800 1000 1200 1400 1600 1800
0 2000 4000 6000 8000 10000



FIGURE 12.14. Vertical concentration of spores: G, SG and UP solutions with
h = H/10 (left) and h = H/100 (right, where only the portion [0, 2000] is shown).
The x-axis represents the vertical coordinate
578 12. Two-Point Boundary Value Problems

12.8 Exercises
1. Consider the boundary value problem (12.1)-(12.2) with f (x) = 1/x. Using
(12.3) prove that u(x) = ’x log(x). This shows that u ∈ C 2 (0, 1) but u(0)
is not de¬ned and u , u do not exist at x = 0 (’: if f ∈ C 0 (0, 1), but not
f ∈ C 0 ([0, 1]), then u does not belong to C 0 ([0, 1])).
2. Prove that the matrix Afd introduced in (12.8) is an M-matrix.
[Hint: check that Afd x ≥ 0 ’ x ≥ 0. To do this, for any ± > 0 set Afd,± =
Afd +±In’1 . Then, compute w = Afd,± x and prove that min1¤i¤(n’1) wi ≥
0. Finally, since the matrix Afd,± is invertible, being symmetric and positive
de¬nite, and since the entries of A’1 are continuous functions of ± ≥ 0,
fd,±
’1
one concludes that Afd,± is a nonnegative matrix as ± ’ 0.]
0
3. Prove that (12.13) de¬nes a norm for Vh .
4. Prove (12.15) by induction on m.
5. Prove the estimate (12.23).
[Hint: for each internal node xj , j = 1, . . . , n ’ 1, integrate by parts (12.21)
to get

„h (xj )
® 
xj xj +h
1 
= ’u (xj ) ’ u (t)(xj ’ h ’ t)2 dt ’ u (t)(xj + h ’ t)2 dt» .
°
h2
xj ’h xj


Then, pass to the squares and sum „h (xj )2 for j = 1, . . . , n ’ 1. On noting
that (a + b + c)2 ¤ 3(a2 + b2 + c2 ), for any real numbers a, b, c, and applying
the Cauchy-Schwarz inequality yields the desired result.]
6. Prove that Gk (xj ) = hG(xj , xk ), where G is Green™s function introduced in
(12.4) and Gk is its corresponding discrete counterpart solution of (12.4).
[Solution: we prove the result by verifying that Lh G = hek . Indeed, for a
¬xed xk the function G(xk , s) is a straight line on the intervals [0, xk ] and
[xk , 1] so that Lh G = 0 at every node xl with l = 0, . . . , k ’ 1 and l =
k + 1, . . . , n + 1. Finally, a direct computation shows that (Lh G)(xk ) = 1/h
which concludes the proof.]
7. Let g = 1 and prove that Th g(xj ) = 1 xj (1 ’ xj ).
2
[Solution: use the de¬nition (12.25) with g(xk ) = 1, k = 1, . . . , n ’ 1 and
recall that Gk (xj ) = hG(xj , xk ) from the exercise above. Then
® 
j n’1
Th g(xj ) = h ° xj (1 ’ xk )»
xk (1 ’ xj ) +
k=1 k=j+1


from which, after straightforward computations, one gets the desired re-
sult.]
8. Prove Young™s inequality (12.40).
¤ vh ∀vh ∈ Vh .
9. Show that vh h h,∞
12.8 Exercises 579

10. Consider the two-point boundary value problem (12.29) with the following
boundary conditions

»0 u(0) + µ0 u(0) = g0 , »1 u(1) + µ1 u(1) = g1 ,

where »j , µj and gj are given data (j = 0, 1). Using the mirror imaging
technique described in Section 12.2.3 write the ¬nite di¬erence discretiza-
tion of the equations corresponding to the nodes x0 and xn .
[Solution:

±1/2 γ0 ± 0 »0 u1 ±0 g0 f0
u0 ’ ±1/2 2 =
node x0 : + + +,
2
h 2 µ0 h h µ0 h 2
±n’1/2 γn ±n »1 un’1 ±n g1 fn
un ’ ±n’1/2 2 =
node xn : + + + .]
2
h 2 µ1 h h µ1 h 2

11. Discretize the fourth-order di¬erential operator Lu(x) = ’u(iv) (x) using
centred ¬nite di¬erences.
[Solution : apply twice the second order centred ¬nite di¬erence operator
Lh de¬ned in (12.9).]
12. Consider problem (12.41) with non homogeneous Neumann boundary con-
ditions ±u (0) = w0 , ±u (1) = w1 . Show that the solution satis¬es prob-
lem (12.43) where V = H1 (0, 1) and the right-hand side is replaced by
(f, v)+w1 v(1)’w0 v(0). Derive the formulation in the case of mixed bound-
ary conditions ±u (0) = w0 , u(1) = u1 .
13. Using Property 1.19 prove that the matrices corresponding to the stabilized
¬nite element method (12.79) using the upwind and SG arti¬cial viscosities
φU P and φSG (see Section 12.5.2) are M-matrices irrespective of h.
[Hint: let us denote respectively by AU P and ASG the two sti¬ness matrices
corrsponding to φU P and φSG . Take v(x) = 1 + x and set vi = 1 + xi ,
i = 0, . . . , n, being xi = ih, h = 1/n. Then, by a direct computation check
that (AU P v)i ≥ β > 0. As for the matrix ASG the same result can be
proved by noting that B(’t) = t + B(t) for any t ∈ R.]
14. Prove that the matrix Afd with entries given by (12.93) is symmetric pos-
itive de¬nite and it is also an M-matrix.
[Solution: to show that Afd is positive de¬nite, proceed as in the corre-
sponding proof in Section 12.2, then proceed as in Exercise 2.]
15. Prove the Green™s formula (12.95).
[Solution: ¬rst, notice that for any u, v su¬ciently smooth, div(v∇u) =
v u + ∇u · ∇v. Then, integrate this relation over „¦ and use the divergence
‚u
theorem div(v∇u) dxdy = v dγ.]
‚n
„¦ ‚„¦
13
Parabolic and Hyperbolic Initial
Boundary Value Problems




The ¬nal chapter of this book is devoted to the approximation of time-
dependent partial di¬erential equations. Parabolic and hyperbolic initial-
boundary value problems will be addressed and either ¬nite di¬erences and
¬nite elements will be considered for their discretization.



13.1 The Heat Equation
The problem we are considering is how to ¬nd a function u = u(x, t) for
x ∈ [0, 1] and t > 0 that satis¬es the partial di¬erential equation

‚u
(13.1)
+ Lu = f 0 < x < 1, t > 0,
‚t
subject to the boundary conditions

u(0, t) = u(1, t) = 0, t>0 (13.2)

and the initial condition

0 ¤ x ¤ 1.
u(x, 0) = u0 (x) (13.3)

The di¬erential operator L is de¬ned as

‚2u
Lu = ’ν . (13.4)
‚x2
582 13. Parabolic and Hyperbolic Initial Boundary Value Problems

Equation (13.1) is called the heat equation. In fact, u(x, t) describes the
temperature at the point x and time t of a metallic bar of unit length
that occupies the interval [0, 1], under the following conditions. Its thermal
conductivity is constant and equal to ν > 0, its extrema are kept at a
constant temperature of zero degrees, at time t = 0 its temperature at
point x is described by u0 (x), and f (x, t) represents the heat production
per unit length supplied at point x at time t. Here we are supposing that
the volumetric density ρ and the speci¬c heat per unit mass cp are both
constant and unitary. Otherwise, the temporal derivative ‚u/‚t should be
multiplied by the product ρcp in (13.1).
A solution of problem (13.1)-(13.3) is provided by a Fourier series. For
instance, if ν = 1 and f ≡ 0, it is given by

2
cn e’(nπ) t sin(nπx)
u(x, t) = (13.5)
n=1

where the coe¬cients cn are the Fourier sine coe¬cients of the initial datum
u0 (x), i.e.
1

cn = 2 u0 (x) sin(nπx) dx, n = 1, 2 . . .
0

If instead of (13.2) we consider the Neumann conditions

ux (0, t) = ux (1, t) = 0, t > 0, (13.6)

the corresponding solution (still in the case where f ≡ 0 and ν = 1) would
be

d0 2
dn e’(nπ) t cos(nπx),
u(x, t) = +
2 n=1

where the coe¬cients dn are the Fourier cosine coe¬cients of u0 (x), i.e.
1

dn = 2 u0 (x) cos(nπx) dx, n = 1, 2 . . .
0

These expressions show that the solution decays exponentially fast in time.
A more general result can be stated concerning the behavior in time of the
energy
1

u2 (x, t) dx.
E(t) =
0
13.1 The Heat Equation 583

Indeed, if we multiply (13.1) by u and integrate with respect to x over the
interval [0, 1], we obtain
1 1 1
‚2u ‚u2
‚u 1
(x, t)u(x, t) dx ’ ν (x, t)u(x, t) dx = (x, t) dx
‚x2
‚t 2 ‚t
0 0 0
1
2 x=1
‚u ‚u
(x, t) dx ’ ν
+ν (x, t)u(x, t)
‚x ‚x x=0
0
1
2
1 ‚u
= E (t) + ν (x, t) dx,
2 ‚x
0

having used integration by parts, the boundary conditions (13.2) or (13.6),
and interchanged di¬erentiation and integration.
Using the Cauchy-Schwarz inequality (8.29) yields
1

f (x, t)u(x, t) dx ¤ (F (t))1/2 (E(t))1/2
0

1
f 2 (x, t) dx. Then
where F (t) = 0

1
2
‚u
dx ¤ 2(F (t))1/2 (E(t))1/2 .
E (t) + 2ν (x, t)
‚x
0

Owing to the Poincar´ inequality (12.16) with (a, b) = (0, 1) we obtain
e
ν
E(t) ¤ 2(F (t))1/2 (E(t))1/2 .
E (t) + 2 2
(CP )

By Young™s inequality (12.40) we have

1
2(F (t))1/2 (E(t))1/2 ¤ γE(t) + F (t),
γ

having set γ = ν/CP . Therefore, E (t) + γE(t) ¤ γ F (t), or, equivalently,
1
2

(eγt E(t)) ¤ γ eγt F (t). Then, integrating from 0 to t we get
1


t
1
E(t) ¤ e’γt E(0) + eγ(s’t) F (s)ds. (13.7)
γ
0

In particular, when f ≡ 0, (13.7) shows that the energy E(t) decays expo-
nentially fast in time.
584 13. Parabolic and Hyperbolic Initial Boundary Value Problems

13.2 Finite Di¬erence Approximation of the Heat
Equation
To solve the heat equation numerically we have to discretize both the x
and t variables. We can start by dealing with the x-variable, following the
same approach as in Section 12.2. We denote by ui (t) an approximation of
u(xi , t), i = 0, . . . , n and approximate the Dirichlet problem (13.1)-(13.3)
by the scheme
. (t) ’ ν (u i’1 (t) ’ 2ui (t) + ui+1 (t)) = fi (t), i = 1, . . . , n ’ 1, ∀t > 0,
u i
h2
∀t > 0,
u0 (t) = un (t) = 0,

ui (0) = u0 (xi ), i = 0, . . . , n,

where the upper dot indicates derivation with respect to time, and fi (t) =
f (xi , t). This is actually a semi-discretization of problem (13.1)-(13.3), and
is a system of ordinary di¬erential equations of the following form

u(t) = ’νAfd u(t) + f (t), ∀t > 0,

(13.8)
u(0) = u0

where u(t) = (u1 (t), . . . , un’1 (t))T is the vector of unknowns, f (t) =
(f1 (t), . . . , fn’1 (t))T , u0 = (u0 (x1 ), . . . , u0 (xn’1 ))T and Afd is the tridi-
agonal matrix introduced in (12.8). Note that for the derivation of (13.8)
we have assumed that u0 (x0 ) = u0 (xn ) = 0, which is coherent with the
boundary condition (13.2).
A popular scheme for the integration of (13.8) with respect to time is the
so-called θ’method. To construct the scheme, we denote by v k the value
of the variable v at time tk = k∆t, for ∆t > 0; then, the θ-method for the
time-integration of (13.8) is
± k+1
’ uk
u
 = ’νAfd (θuk+1 + (1 ’ θ)uk ) + θf k+1 + (1 ’ θ)f k ,

∆t
(13.9)

 k = 0, 1, . . .
0
u = u0

or, equivalently,

(I + νθ∆tAfd ) uk+1 = (I ’ ν∆t(1 ’ θ)Afd ) uk + gk+1 , (13.10)

where gk+1 = ∆t(θf k+1 + (1 ’ θ)f k ) and I is the identity matrix of order
n ’ 1.
For suitable values of the parameter θ, from (13.10) we can recover some
familiar methods that have been introduced in Chapter 11. For example, if
θ = 0 the method (13.10) coincides with the forward Euler scheme and we
13.2 Finite Di¬erence Approximation of the Heat Equation 585

can get uk+1 explicitly; otherwise, a linear system (with constant matrix
I + νθ∆tAfd ) needs be solved at each time-step.
Regarding stability, assume that f ≡ 0 (henceforth gk = 0 ∀k > 0),
so that from (13.5) the exact solution u(x, t) tends to zero for every x
as t ’ ∞. Then we would expect the discrete solution to have the same
behaviour, in which case we would call our scheme (13.10) asymptotically
stable, this being coherent with what we did in Chapter 11, Section 11.1
for ordinary di¬erential equations.
If θ = 0, from (13.10) it follows that

uk = (I ’ ν∆tAfd )k u0 , k = 1, 2, . . .

From the analysis of convergent matrices (see Section 1.11.2) we deduce
that uk ’ 0 as k ’ ∞ i¬

ρ(I ’ ν∆tAfd ) < 1. (13.11)

On the other hand, the eigenvalues of Afd are given by (see Exercise 3)

4
i = 1, . . . , n ’ 1.
sin2 (iπh/2),
µi = 2
h

Then (13.11) is satis¬ed i¬

12
∆t < h.


As expected, the forward Euler method is conditionally stable, and the
time-step ∆t should decay as the square of the grid spacing h.
In the case of the backward Euler method (θ = 1), we would have from
(13.10)

k
uk = (I + ν∆tAfd )’1 u0 , k = 1, 2, . . .

Since all the eigenvalues of the matrix (I + ν∆tAfd )’1 are real, positive and
strictly less than 1 for every value of ∆t, this scheme is unconditionally
stable. More generally, the θ-scheme is unconditionally stable for all the
values 1/2 ¤ θ ¤ 1, and conditionally stable if 0 ¤ θ < 1/2 (see Section
13.3.1).
As far as the accuracy of the θ-method is concerned, its local truncation
error is of the order of ∆t + h2 if θ = 1 while it is of the order of ∆t2 + h2 if
2
θ = 1 . The method corresponding to θ = 1/2 is frequently called the Crank-
2
Nicolson scheme and is therefore unconditionally stable and second-order
accurate with respect to both ∆t and h.
586 13. Parabolic and Hyperbolic Initial Boundary Value Problems

13.3 Finite Element Approximation of the Heat
Equation
The space discretization of (13.1) can also be accomplished using the Galerkin
¬nite element method by proceeding as in Chapter 12 in the elliptic case.
First, for all t > 0 we multiply (13.1) by a test function v = v(x) and
integrate over (0, 1). Then, we let V = H1 (0, 1) and ∀t > 0 we look for a
0
function t ’ u(x, t) ∈ V (brie¬‚y, u(t) ∈ V ) such that

1
‚u(t)
∀v ∈ V, (13.12)
vdx + a(u(t), v) = F (v)
‚t
0

1
with u(0) = u0 . Here, a(u(t), v) = 0 ν(‚u(t)/‚x) (‚v/‚x) dx and F (v) =
1
f (t)vdx are the bilinear form and the linear functional respectively as-
0
sociated with the elliptic operator L and the right hand side f . Notice that
a(·, ·) is a special case of (12.44) and that the dependence of u and f on
the space variable x will be understood henceforth.
Let Vh be a suitable ¬nite dimensional subspace of V . We consider the
following Galerkin formulation: ∀t > 0, ¬nd uh (t) ∈ Vh such that

1
‚uh (t)
vh dx + a(uh (t), vh ) = F (vh ) ∀vh ∈ Vh (13.13)
‚t
0


where uh (0) = u0 h and u0 h ∈ Vh is a convenient approximation of u0 .
Problem (13.13) is referred to as a semi-discretization of (13.12) since it is
only a space discretization of the heat equation.
Proceeding in a manner similar to that used to obtain the energy estimate
(13.7), we get the following a priori estimate for the discrete solution uh (t)
of (13.13)
t
1
Eh (t) ¤ e’γt Eh (0) + eγ(s’t) F (s)ds,
γ
0

1
where Eh (t) = 0 u2 (x, t) dx.
h
As for the ¬nite element discretization of (13.13), we introduce the ¬nite
element space Vh de¬ned in (12.57) and consequently a basis {•j } for Vh
as already done in Section 12.4.5. Then, the solution uh of (13.13) can be
sought under the form
Nh
uh (t) = uj (t)•j ,
j=1
13.3 Finite Element Approximation of the Heat Equation 587

where {uj (t)} are the unknown coe¬cients and Nh is the dimension of Vh .
Then, from (13.13) we obtain
« 
. (t)• • dx + a 
1 Nh Nh
uj (t)•j , •i  = F (•i ),
u i = 1, . . . , Nh
j j i
0 j=1 j=1


that is,

. (t) 1
Nh Nh
u •j •i dx + uj (t)a(•j , •i ) = F (•i ), i = 1, . . . , Nh .
j
j=1 j=1
0

Using the same notation as in (13.8) we obtain
.
Mu(t) + Afe u(t) = ffe (t) (13.14)
1
where Afe = (a(•j , •i )), ffe (t) = (F (•i )) and M = (mij ) = ( 0 •j •i dx) for
i, j = 1, . . . , Nh . M is called the mass matrix. Since it is nonsingular, the
system of ODEs (13.14) can be written in normal form as
.
u(t) = ’M’1 Afe u(t) + M’1 ffe (t). (13.15)

To solve (13.15) approximately we can still apply the θ-method and obtain

uk+1 ’ uk
+ Afe θuk+1 + (1 ’ θ)uk = θffe + (1 ’ θ)ffe .
k+1 k
M (13.16)
∆t
As usual, the upper index k means that the quantity at hand is computed at
time tk . As in the ¬nite di¬erence case, for θ = 0, 1 and 1/2 we respectively
obtain the forward Euler, backward Euler and Crank-Nicolson methods,
where the Crank-Nicolson method is the only one which is second-order
accurate with respect to ∆t.
For each k, (13.16) is a linear system whose matrix is
1
K= M + θAfe .
∆t
Since M and Afe are symmetric and positive de¬nite, the matrix K is also
symmetric and positive de¬nite. Thus, its Cholesky decomposition K =
HT H where H is upper triangular (see Section 3.4.2) can be carried out at
t = 0. Consequently, at each time step the following two linear triangular
systems, each of size equal to Nh , must be solved, with a computational
2
cost of Nh /2 ¬‚ops
±
 HT y = 1 M ’ (1 ’ θ)A uk + θf k+1 + (1 ’ θ)f k ,
fe fe
fe
∆t

Huk+1 = y.
588 13. Parabolic and Hyperbolic Initial Boundary Value Problems

When θ = 0, a suitable diagonalization of M would allow to decouple
the system equations (13.16). The procedure is carried out by the so-called
mass-lumping in which we approximate M by a nonsingular diagonal matrix
M. In the case of piecewise linear ¬nite elements M can be obtained using
the composite trapezoidal formula over the nodes {xi } to evaluate the
1
integrals 0 •j •i dx, obtaining mij = hδij , i, j = 1, . . . , Nh (see Exercise
˜
2).


Stability Analysis of the θ-Method
13.3.1
Applying the θ-method to the Galerkin problem (13.13) yields

<<

. 23
( 26)



>>