ńņš. 23 |

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 =

2Īµ

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

6Āµ

ļ£“

ļ£³ 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.

2Ī½

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 |