<< . .

. 82
( : 95)



. . >>


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

<< . .

. 82
( : 95)



. . >>