Analysis of the iterative penalty method for the Stokes equations

Applied Mathematics Letters 19 (2006) 1024–1028 www.elsevier.com/locate/aml

Analysis of the iterative penalty method for the Stokes equations✩ Cheng Xiao-liang a,∗ , Abdul Wasim Shaikh b,c a Department of Mathematics, Zhejiang University, Hangzhou 310028, PR China b Department of Mathematics, Zhejiang University, PR China c University of Sindh, Jamshoro, Pakistan

Received 18 August 2004; received in revised form 8 June 2005; accepted 11 October 2005

Abstract In this work we propose an iterative penalty method for addressing the Stokes equations. We can use a “not very small” penalty parameter to avoid the unstable computation by iteration. The Numerical experiments show that the algorithm is very effective. c 2005 Elsevier Ltd. All rights reserved.  Keywords: Penalty method; Stokes equations; Stable computation; Iterative scheme; Error estimates

1. Introduction In this work, we consider the numerical solution of the Stokes problem for viscous incompressible fluid flows using the penalty finite element method. The model problem to be studied is finding the velocity u = (u 1 , u 2 ) ∈ [H01(Ω )]2 and the pressure p ∈ L 20 (Ω ) such that
(∇ u, ∇ v) − (div v, p) = ( f, v), ∀ v ∈ [H01(Ω )]2 , (1.1) (div u, q) = 0, ∀ q ∈ L 20 (Ω ). Here, Ω ⊂ R 2 is a bounded domain, f is the given body force. The symbol (·, ·) denotes the usual inner product in L 2 (Ω ) or [L 2 (Ω )]2 . The space L 20 (Ω ) consists of all the L 2 (Ω )-functions whose mean values in Ω are zero. The Sobolev space H01(Ω ) is the set of all the L 2 (Ω )-functions whose first order partial derivatives are also in L 2 (Ω ), and whose traces on the boundary ∂Ω vanish. We use | · |1 for its semi-norm, which is equivalent to the norm  · 1 in Sobolev space H01(Ω ). For f ∈ [H −1(Ω )]2 , the Stokes problem (1.1) has a unique solution, cf. [1,2]. It is easy to see that the system (1.1) is equivalent to the following constrained minimization formulation: Find u ∈ H01(Ω ) and div u = 0 such that min

v∈H01 (Ω ),div

1 (∇ u, ∇ v ) − ( f, v). v=0 2

✩ The project is supported by the National Natural Science Foundation of China (Grant No. 10471129). ∗ Corresponding author.

E-mail address: [email protected] (X.-l. Cheng). c 2005 Elsevier Ltd. All rights reserved. 0893-9659/$ - see front matter  doi:10.1016/j.aml.2005.10.021

X.-l. Cheng, A.W. Shaikh / Applied Mathematics Letters 19 (2006) 1024–1028

1025

We consider the penalty method related to reduced integration for the above problem: Find uε ,h ∈ Vh ⊂ [H01(Ω )]2 such that 1 ∀ vh ∈ Vh , (1.2) (∇ uε,h , ∇ vh ) + I (div uε,h , div vh ) = ( f, vh ), ε where I is a numerical integration operator, ε > 0 is a small parameter. Under certain assumptions, the reduced integration penalty method (1.2) is equivalent to the following mixed finite element method: find ( u ε,h , pε,h ) ∈ Vh × Q h such that
(∇ uε,h , ∇ vh ) − (div vh , pε,h ) = ( f, vh ), ∀ vh ∈ Vh , (1.3) ε( pε,h , qh ) + (div u ε,h , qh ) = 0, ∀ qh ∈ Q h , for some finite element space Q h ⊂ L 20 (Ω ). We suppose that the spaces Vh and Q h satisfy the Babuˇska–Brezzi condition [3,4] β0 qh 0 ≤

sup 0 =vh ∈Vh

|(qh , div vh )| , | vh |1

∀qh ∈ Q h ,

(1.4)

where β0 > 0 is a constant independent of h. It is well known [1,2] that the error estimate of the penalty method (1.2) or (1.3) for approximating (1.1) is u − vh 1 +infqh ∈Q h  p−qh 0 . Thus, ε must be sufficiently small to yield an accurate O(ε+ Rh ), where Rh = infvh ∈Vh  approximation of the solution of (1.1). On the other hand, the condition number of the numerical discretizations of (1.2) is O(ε−1 h −2 ) so the scheme becomes unstable due to round-off error, if ε is too small. In [5], a modified penalty method via the extrapolation umn,h = uεn ,h − εn

uεm ,h − uεn ,h , εm − εn

(1.5)

where uεm ,h , uεn ,h are the solutions of penalty method (1.2) with respect to εm and εn is proposed. It is proved that the error estimates are O(εm εn + Rh ) in H 1-norm. This allows us to use a “not very small” penalty parameter. In this work, we propose the iterative penalty method using (0)

(∇ uε,h , ∇ vh ) +

1 (0) I (div uε,h , div vh ) = ( f, vh ), ε

∀ vh ∈ Vh ,

(1.6)

and for k = 1, 2, . . . (k)

(∇ uε,h , ∇ vh ) +

1 (k) (k−1) I (div uε,h , div vh ) = (∇ uε,h , ∇ vh ), ε

∀ vh ∈ Vh .

(1.7)

We will prove that the error estimate is O(εk+1 + Rh ) for any positive integer k. Eqs. (1.6) and (1.7) are the same linear system with a different right hand side. For a not very small fix parameter ε, we can use the fast Poisson solver as the preconditioner of the linear system (1.6) and (1.7). It is an effective algorithm for the Stokes equations. Numerical experiments will be presented in the last section. 2. Error estimates To analyze the algorithm (1.6) and (1.7), we first rewrite it in the equivalent mixed penalty formulation: Find (0) (0) ( u ε,h , pε,h ) ∈ Vh × Q h such that  (0) h ) − (div vh , pε,h ) = ( f, vh ), ∀ vh ∈ Vh , (∇ u(0) ε,h , ∇ v (2.1) (0) (0) ε( pε,h , qh ) + (div u ε,h , qh ) = 0, ∀ qh ∈ Q h , (k)

(k)

and for k = 1, 2, . . ., we find ( u ε,h , pε,h ) ∈ Vh × Q h such that  (k) (∇ u(k) h ) − (div vh , pε,h ) = ( f, vh ), ∀ vh ∈ Vh , ε,h , ∇ v (k) (k) (k−1) ε( pε,h , qh ) + (div u ε,h , qh ) = ε( pε,h , qh ), ∀ qh ∈ Q h .

(2.2)

1026

X.-l. Cheng, A.W. Shaikh / Applied Mathematics Letters 19 (2006) 1024–1028

We state the result for the error estimates for the initial value of the algorithm. This is the standard penalty method for the Stokes equation (1.1). The proof can be found in [1,2,6–8]. (0)

(0)

(k)

(k)

Theorem 2.1. Let ( u , p) be the solution of equation (1.1) and ( u ε,h , pε,h ) be the solution of Eq. (2.1), and let the BB-condition (1.4) hold, then   (0) (0)  u − u ε,h 1 +  p − pε,h 0 ≤ C0 inf  (2.3) u − vh 1 + inf  p − qh 0 + ε p0 , qh ∈Q h

vh ∈Vh

where the constant C0 is independent of h and parameter ε. Our main results are the following error estimates. Theorem 2.2. Let ( u , p) be the solution of equation (1.1) and ( u ε,h , pε,h ) (k = 0, 1, 2, . . .) be the solution of Eq. (2.2), and let the BB-condition (1.4) hold, then for k ≥ 1   (k) (k) (k−1) u − vh 1 + inf  p − qh 0 + ε p − pε,h 0 ,  u − u ε,h 1 +  p − pε,h 0 ≤ C0 inf  (2.4) qh ∈Q h

vh ∈Vh

where the constant C0 is independent of h and parameter ε. h , the projections onto space Vh , Q h , Proof. Let ( u , p) be the solution of equation (1.1). We denote by πVh , π Q respectively, that satisfy the following relations:

u − vh 1 ,  u − πVh u 1 = inf 

(2.5)

h  p − πQ p0 = inf  p − qh 0 .

(2.6)

vh ∈Vh

qh ∈Q h

Let (k) ξ (k) = u − u ε,h ,

(k)

η(k) = p − pε,h .

(2.7)

From (1.1), (2.1) and (2.2), for k ≥ 1, we have the following relations:
∀ vh ∈ Vh , (∇ ξ (k) , ∇ vh ) − (div vh , η(k) ) = 0, ε(η(k) , qh ) + (div ξ (k) , qh ) = ε(η(k−1) , qh ), ∀ qh ∈ Q h .

(2.8)

Then we have from the first equation of (2.8) 2

h h |πVh u − u(k)  − u(k)  − u(k) ε,h |1 = (∇(π V u ε,h ), ∇(π V u ε,h )) h u − u(k)  − u (k) ≤ (∇(πVh u − u), ∇(πVh u − u(k) ε,h )) + (∇( ε,h ), ∇(π V u ε,h )) (k)

(k)

≤ CπVh u − u1 πVh u − uε,h 1 + (η(k) , div (πVh u − uε,h )).

(2.9)

Furthermore, from the second equation of (2.8), we obtain h h  − u(k) (η(k) , div (πVh u − u(k) ε,h )) = ( p − π Q p, div(π V u ε,h )) (k) (k) h h + (π Q p − pε,h , div (πVh u − u )) + (π Q p − pε,h , div ξ (k) ) (k) h h h ≤ C( p − π Q p0 πVh u − u (k)  − u1 ) ε,h 1 + π Q p − pε,h 0 π V u (k) (k) h h − ε(η(k) , π Q p − pε,h ) + ε(η(k−1) , π Q p − pε,h )

(2.10)

and (k) 2

(k)

(k)

h h h h −ε(η(k) , π Q p − pε,h ) = −επ Q p − pε,h 0 − ε( p − π Q p, π Q p − pε,h ).

(2.11)

From (2.9)–(2.11), we get (k) 2

(k) 2

h |πVh u − uε,h |1 + επ Q p − pε,h 0 (k)

h ≤ C(πVh u − u1 +  p − π Q p0 )πVh u − u ε,h 1 (k)

h h + (CπVh u − u1 + ε p − π Q p0 + εη(k−1) 0 )π Q p − pε,h 0 .

(2.12)

X.-l. Cheng, A.W. Shaikh / Applied Mathematics Letters 19 (2006) 1024–1028

1027

Using the BB-condition (1.4) and the first equation of (2.8), we get h p π Q

−

(k) pε,h 0

1 ≤ β0

h p − p (k) , div |(π Q vh )| ε,h

sup

 vh 1

0 =vh ∈Vh h |(π Q p

− p, div vh ) + (∇ ξ (k) , ∇ vh )|

≤

1 β0

≤

1 h ( p − π Q p0 + ξ (k) 1 ). β0

sup

 vh 1

0 =vh ∈Vh

(2.13)

Thus, the inequality (2.12) implies (k) 2

(k) 2

h p − pε,h 0 |πVh u − uε,h |1 + επ Q (k)

h ≤ C(πVh u − u1 +  p − π Q p0 )πVh u − u ε,h 1 h h + (CπVh u − u1 + ε p − π Q p0 + εη(k−1) 0 )( p − π Q p0 + ξ (k) 1 ).

By the identity ab ≤ ca 2 +

1 2 4c b

(2.14)

for any c > 0, for example, the last term in (2.14) becomes

1 2 2 εη(k−1) 0 ξ (k) 1 ≤ Cε2 η(k−1) 0 + |ξ (k) |1 , 8 where we have used the equivalence of the norm and semi-norm in space H01(Ω ), ξ (k) 1 ≤ C|ξ (k) |1 . By applying a similar argument to the other terms in (2.14), we can prove that (using triangle inequalities | u − u (k) ε,h |1 ≤ (k) (k) h h | u − πVh u |1 + |πVh u − u(k) ε,h |1 and  p − pε,h 0 ≤  p − π Q p0 + π Q p − pε,h 0 )   (k) 2 (k) 2 2 2 2 (k−1) 2  u − u ε,h 1 +  p − pε,h 0 ≤ C inf  u − vh 1 + inf  p − qh 0 + ε η 0 , qh ∈Q h

vh ∈Vh

where the constant C (and thus C0 = theorem.


(2.15)

√ C) is independent of h, parameter ε and k. This implies the results of the

Remark 2.3. From the proof in Theorem 2.2, we can also get the result of Theorem 2.1. If we choose the penalty parameter ε such that C0 ε < 1,

(2.16)

then we have the error estimates (k)

(k)

 u − u ε,h 1 +  p − pε,h 0 ≤

C0 1 − C0 ε



 u − vh 1 + inf  p − qh 0 + (C0 ε)k+1  p0 . inf 

vh ∈Vh

qh ∈Q h

(2.17)

Remark 2.4. For the Navier–Stokes equation, we can also propose the following iterative penalty method based on the formulation of [6,9,10] 1 (∇ u(0) h ) + (( u (0) u (0) h ) + (∇ · u (0) (0) h ) ε,h , ∇ v ε,h · ∇) ε,h , v ε,h , u ε,h · v 2 1 + I (div u (0) h ) = ( f, vh ), ∀ vh ∈ Vh , ε,h , div v ε and for k = 1, 2, . . .

(2.18)

1 1 h ) + (( u (k) u (k) h ) + (∇ · u (k) (k) h ) + I (div u(k) h ) (∇ u(k) ε,h , ∇ v ε,h · ∇) ε,h , v ε,h , u ε,h · v ε,h , div v 2 ε (k−1) (k−1) (k−1) u ε,h · ∇) u ε,h , vh ) = (∇ uε,h , ∇ vh ) + (( 1 (k−1) (k−1) + (∇ · u ε,h , uε,h · vh ), 2

∀ vh ∈ Vh .

(2.19)

1028

X.-l. Cheng, A.W. Shaikh / Applied Mathematics Letters 19 (2006) 1024–1028

Table 3.1 The errors of the iteration penalty method with ε = 0.1

E (0) E (1) E (2) E (3) E (4) E (5) E (6) E (7) E (8) E (9) E (10)

h = 1/4

h = 1/8

h = 1/16

0.01580 0.005975 0.004094 0.003726 0.003662 0.003668 0.003667 0.003664 0.003663 0.003663 0.003662

0.01459 0.003518 0.001391 0.00095551 0.0008412 0.0008054 0.0007909 0.0007835 0.0007791 0.0007762 0.0007742

0.01402 0.003022 0.0008487 0.0003412 0.0002312 0.0002021 0.0001949 0.0001969 0.0001974 0.0001974 0.0001973

For the nonlinear system (2.18) and (2.19), we can use some linearization method such as the Newton method. See [6] for the analysis and numerical results. 3. Numerical test In [7,6,5] there are many numerical results for the penalty method. Here we only show the influence of the penalty parameter. We use the bilinear-constant velocity–pressure element. In the penalty formulation, we let Vh be the bilinear element and have one-point Gauss integration in each element. This bilinear-constant pair satisfies the BB-condition (1.4), but the spurious pressure modes are filtered out automatically [11–14]. We consider the exact velocity   2 x (1 − x)2 (2y − 6y 2 + 4y 3 ) u = y 2 (1 − y)2 (−2x + 6x 2 − 4x 3 ) and the pressure p = x 2 − y2 (k)

on the domain Ω = (0, 1) × (0, 1). We present the maximum error E (k) = max | u − uε,h |, ε = 0.1 for h = 1/4, h = 1/8 and h = 1/16, see Table 3.1. We can see that after a few iterations the errors tend to the fixed errors which are caused by the finite element approximation of mesh size h. This is consistent with the error of the penalty method (1.2) as ε tends to 0. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]

F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, 1991. V. Girault, P.A. Raviart, Finite Element Methods for Navier–Stokes Equations, Theory and Algorithms, Springer-Verlag, 1986. I. Babuˇska, The finite element method with Lagrangian multipliers, Numer. Math. 20 (1973) 179–192. F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers, RAIRO Numer. Anal. 8 (1974) 129–151. S.Y. Lin, Y.S. Chin, T.M. Wu, A modified penalty method for Stokes equations and its applications to Navier–Stokes equations, SIAM J. Sci. Comput. 16 (1995) 1–19. G.F. Carey, R. Krishnan, Penalty approximation of Stokes flow, Comput. Methods Appl. Mech. Engrg. 35 (1982) 169–206. J.T. Oden, N. Kikuchi, Y.J. Song, Penalty finite element methods for the analysis of Stokesian flows, Comput. Methods Appl. Mech. Engrg. 31 (1982) 297–329. R. Stenberg, Error analysis of some finite element methods for the Stokes problem, Math. Comput. 54 (1990) 495–508. R. Temam, Navier–Stokes Equations, Theory and Numerical Analysis, North-Holland, Amsterdam, 1977. G.F. Carey, J.T. Oden, Finite Elements, Fluid Mechanics—volume VI, Prentice-Hall, Englewood Cliffs, NJ, 1986. J.M. Boland, R.A. Nicolaides, On the stability of bilinear-constant velocity-pressure finite elements, Numer. Math. 44 (1984) 219–222. C. Johnson, J. Pitk¨aranta, Analysis of some mixed finite element methods related to reduced integration, Math. Comput. 38 (1982) 375–400. X.L. Cheng, W. Han, H.C. Huang, Analysis of some mixed elements for the Stokes problem, J. Comput. Appl. Math. 85 (1997) 19–35. J.-s. Jiang, X.-l. Cheng, On the stability of biquadratic–bilinear velocity–pressure finite element, J. Comput. Math. 10 (1992) 167–170.