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