Superconvergence and asymptotic expansions for bilinear finite volume element approximation on non-uniform grids

Journal of Computational and Applied Mathematics 321 (2017) 323–335

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Superconvergence and asymptotic expansions for bilinear finite volume element approximation on non-uniform grids✩ Cunyun Nie a, *, Shi Shu b,c , Haiyuan Yu b,c , Wenhua Xia a a b c

The School of Science, Hunan Institute of Engineering, 411104, China Hunan Key Laboratory for Computation & Simulation in Science and Engineering, Xiangtan University, Hunan 411105, China Key Laboratory of Intelligent Computing & Information Processing of Ministry of Education, Xiangtan University, Hunan 411105, China

article

info

Article history: Received 9 March 2016 Received in revised form 1 November 2016 Keywords: Isoparametric bilinear finite volume element scheme Non-uniform grids Asymptotic expansion Superconvergence

a b s t r a c t We initially derive an asymptotic expansion and a high accuracy combination formula of the derivatives in the sense of pointwise for the isoparametric bilinear finite volume element scheme by employing the energy-embedded method on some non-uniform grids. Furthermore, we prove that the approximate derivatives are convergent of order two. Numerical examples confirm theoretical results. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Finite volume (FV) method has been one of the most popular numerical methods for solving partial differential equations because of its some attractive properties, such as preserving local conservation of mass, momentum and energy. The finite volume element (FVE) method is one important branch of FV method. The isoparametric bilinear finite volume element scheme is firstly introduced by Li and Zhu in 1982, and its error estimate in H 1 norm on quadrilateral grids is given in [1]. The box scheme is constructed by Schmidt and Kiel in 1993, and the saturated convergent order in H 1 norm and the superconvergent result on parallelogram grids are obtained in [2]. Later, ‘‘Covolume Method’’ is proposed and applied in computational fluid dynamic problems by Porsching and Chou in [3,4]. Simultaneously, some symmetric FVE schemes [5,6], high order FVE schemes [7–9] and new FVE schemes for three dimensional problems [10,11] were presented by some researchers. Recently, Lin and Liu gave some overview on FVE method in [12], and some new immersed finite volume element schemes appeared and became popular [13,14]. Recently, Li and his companion proved the optimal L2 error results for the isoparametric bilinear finite volume element scheme in [15–18]. Although the superconvergence for finite element methods is abundantly studied in [19–21], there are only some researches on superconvergence about the isoparametric bilinear finite volume element in [1,18,22,23]. Furthermore, they are almost in the sense of average instead of pointwise. Nie and Yu presented some superconvergence results in the sense of pointwise on uniform grids in [24,25]. To our knowledge, we have not found some discussions for this scheme about superconvergence in the pointwise on non-uniform grids. In present paper, we derive an asymptotic expansion for the isoparametric bilinear finite volume element solution on some non-uniform grids. In the derivation, we obtain the integral formula for the bilinear functional A(u − uI , v ), and ✩ This work was partially supported by NSFC Project (Grant No. 11031006, 11171281), and the Project of Scientific Research Fund of Hunan Provincial Education Department (Grant No. 14B044, 14A034) and Program for Changjiang Scholars and Innovative Research Team in University of China (No. IRT1179) Corresponding author. E-mail address: [email protected] (C.Y. Nie).

*

http://dx.doi.org/10.1016/j.cam.2016.12.024 0377-0427/© 2016 Elsevier B.V. All rights reserved.

324

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335

Fig. 1. (a) Quadrilateral partition Ωh . (b) Dual element bPi . (c) Transformation ψ .

introduce a proper auxiliary variational problem, the discrete Green function, and apply the energy-embedded method. Furthermore, we derive a high accuracy combination formula of the derivatives in the sense of pointwise, and prove that the approximate derivatives are convergent of order two. Numerical examples confirm theoretical results. The remainder of this paper is organized as follows. In Section 3, we introduce the isoparametric bilinear finite volume element scheme and some convergent results. In Section 3, we derive the asymptotic expansion for the finite volume element solution. In Section 4, we present a high accuracy combination formula of the approximate derivatives in the sense of pointwise on non-uniform grids and the corresponding superconvergence. Finally, we display numerical experiments to support our conclusions. 2. The isoparametric bilinear finite volume element scheme We mainly focus on the following model problem

⎧ ⎪ ⎨−∇ · (a(x)∇ u) = f (x), x = (x, y) ∈ Ω , u = α (x), x ∈ ∂ Ω1 , (1) ∂ u ⎪ ⎩ = β (x), x ∈ ∂ Ω2 , ∂n where Ω = Ω1 ∪ Ω1 ⊂ R2 is a convex polygonal domain with boundary ∂ Ω = ∂ Ω1 ∪ ∂ Ω1 , f (x) ∈ L2 (Ω ), and a(x) ∈ C 1 (Ω ) satisfies a(x) ≥ a0 , a0 is a positive constant, and α (x), β (x) are some given functions. Let Ωh = {Ki , 1 ≤ i ≤ M } be the quadrilateral partition of Ω (see Fig. 1(a)), and D = {Pi = (x1i , x2i ), 1 ≤ i ≤ N } be the set of partition nodes in Ωh , where M and N are, respectively, the numbers of elements and nodes. Denote Ωh∗ = {KP∗ , 1 ≤ i ≤ N } i as the dual partition of Ωh , where KP∗ is the dual element (also called control volume) about node Pi (see Fig. 1(b)). i In this paper, we always assume that the step size and partition number along xi , i = 1, 2 axis direction are hi , Ni , respectively, which satisfy N = (N1 + 1)(N2 + 1), N1 h1 = N2 h2 = 1, and quasi-uniform property c1 h1 ≤ h2 ≤ c2 h1 , where c1 , c2 are two positive constants. For any partition node Pˆ ij (xˆ i , yˆ j ) ∈ Ω h , 0 ≤ i ≤ N1 , 0 ≤ j ≤ N2 , where xˆ i = i h1 , yˆ j = j h2 . We will consider the following disturbed grids for any partition node Pij (xi , yj ) 0 ≤ i ≤ N1 , 0 ≤ j ≤ N2 satisfy some disturbed condition

|Pij − Pˆ ij | ≤ Ch3 ,

(2)

where C is some positive constant independent of h. In the following, we denote still the above disturbed partition as Ωh . We introduce an invertible bilinear transformation FK which maps referred element Kˆ onto convex element K (see Fig. 1(c)), where Kˆ = (0, 1) × (0, 1) is a unit square on the ξ − η plane, and x = x1 + a1 ξ + a2 η + a3 ξ η, y = y1 + b1 ξ + b2 η + b3 ξ η,

{

(3)

where (xi , yi ), i = 1, 2, 3, 4 are four nodal coordinates of element K , and a1 = x 2 − x 1 ,

a2 = x4 − x1 ,

a3 = x1 − x2 + x3 − x4 ,

b1 = y2 − y1 ,

b2 = y4 − y1 ,

b3 = y1 − y2 + y3 − y4 .

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335

325

One can easily see that

( ) xξ yη

( ) ξ = J(ξ , η) , η

(4)

where J(ξ , η) =

(

a1 + a3 η b1 + b3 η

a2 + a3 ξ b2 + b3 ξ

) (5)

is the Jacobi matrix. Lemma 2.1 (See [18]). Assume that for any element K ∈ Ωh whose four nodes are Pl (xl , yl ), l = 1, 2, 3, 4, h = max{h1 , h2 }, and denote the areas of K , △P1 P4 P3 and △P1 P4 P2 as S , S143 and S142 , respectively, and assume that JK (ξ , η)−1 ≜ det(J −1 (ξ , η)) is the determinant of the inverse of J(ξ , η). Then we have S = h1 h2 (1 + O(h2 )), S143 − S142 = O(h4 ),

|P1 P3 |2 − |P2 P4 |2 = O(h4 ), |P2 P3 |2 − |P1 P4 |2 = O(h4 ), and Jk−1 (ξ , η) = (h1 h2 )−1 (1 + O(h2 )). Furthermore, for integers k, k1 , k2 and smooth function u(x, y),

∂ ku ∂ ku k = h (1 + O(h2 )), ∂ξ k1 ∂ηk2 ∂ xk1 ∂ yk2 where k = k1 + k2 . Remark 2.1. To be convenient, we only discuss the case of Problem (1) with homogeneous Dirichlet boundary conditions in the following. In this case, ∂ Ω1 = ∂ Ω , ∂ Ω2 = Φ , and α (x) = 0. For nonhomogeneous Dirichlet boundary conditions, it can be transferred to the homogeneous case by introducing some transformation and extension function. For Neumann boundary conditions, although the discussion becomes complex because of the added line integral, the superconvergent results similarly hold. Numerical results in Section 5 confirm this conclusion. Let VB and Vh , respectively, be the trial and test function spaces VB = {v ∈ L2 (Ω ) : v|K ∗ = constant , KP∗i ∈ Ωh∗ }, Pi

and

¯ ) : uh |E = PKˆ ◦ FK−1 , PKˆ ∈ Pˆ 1,1 , K ∈ Ωh , uh |∂ Ω = 0}, Vh = {uh ∈ C (Ω where Pˆ 1,1 is the set of bilinear functions on Kˆ , and PKˆ ◦ FK−1 is the composite function of PKˆ and FK−1 . We introduce two interpolation operators in the following. The first one is Π ∗ : Vh ↦ → VB ,

Π ∗ v (x) =

{

v (Pi ), x ∈ KP∗i , ∀ v (x) ∈ Vh , 0, otherwise,

(6)

and the second one is Π : H01 (Ω ) ∩ H 2 (Ω ) → Vh ,

Πu =

N ∑

ui φi , u ∈ H01 (Ω ) ∩ H 2 (Ω ),

(7)

i=1

where ui = u(Pi ), and φi is the Lagrange interpolation basis function about node Pi . In the following, we also denote Π u as uI for the sake of convenience. Lemma 2.2 (See [18]). If u ∈ H01 (Ω ) ∩ H 2 (Ω ), then

|u − Π u|0 ≲ Ch2−m |u|2 , m = 0, 1,

|u − Π ∗ u|0 ≲ Ch|u|1 .

The isoparametric bilinear finite volume element solution uh ∈ Vh of problem (1) satisfies A(uh , v ) = (f , Π ∗ v ), ∀v ∈ Vh ,

(8)

326

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335

where A(uh , v ) = −

∑ ∫

a ∂ KP∗

KP∗ ∈Ωh∗ i

K

(f , Π ∗ v ) =

i=1

K ∗ ∈Ω ∗ Pi

h

a ∂ KP∗ ∩K i

∑ ∫ KP∗

(9)

K ∈Ωh

i

4 ∫ ⏐ ∑ ⏐ A(uh , v )⏐ = −

⏐ ∑ ∂ uh ∗ ⏐ Π v ds = A(uh , v )⏐ , K ∂n ∂ uh ∗ Π v ds, ∂n

f Π ∗ v dx.

i

We define the discrete H 1 semi-norm on space Vh as follows

) 12

( |w|1,h =

∑

|w|

, w ∈ Vh ,

2 1,h,K

(10)

K

where |w|21,h,E = (w2 − w1 )2 + (w3 − w2 )2 + (w4 − w3 )2 + (w1 − w4 )2 , wi = w (Pi ). Lemma 2.3 (See [18]). The semi-norm |w|1,h is equivalent with |w|1 , i.e., there exist two positive constants β1 , β2 independent of the mesh size h, such that

β1 |w|1,h ≤ |w|1 ≤ β2 |w|1,h ,

∀w ∈ Vh .

The following results (see [1,16,18]) hold. Lemma 2.4. The bilinear functional A(·, ·) defined by (9) satisfies A(uh , uh ) ≳ |uh |21 , A(uh , vh ) ≲ |uh |1 |vh |1 , ∀ uh , vh ∈ Vh .

(11)

Lemma 2.5. Let u ∈ H01 (Ω ) ∩ H 3 (Ω ) and uh ∈ Vh be the solutions of problems (1) and (8), respectively. Then we have

∥u − uh ∥1 ≲ h|u|2 , ∥u − uh ∥0 ≲ h2 ∥u∥3 .

(12)

3. Asymptotic expansion To obtain the asymptotic expansion, we firstly present some lemmas.

ˆ and uˆ I = Lemma 3.1 (See [26]). Assume that uˆ I is the bilinear interpolation function of uˆ on E, ˆ interpolation basis function about node Pl (see Fig. 2(b)). We have (1) If uˆ ∈ P2 (Kˆ ), then uˆ − uˆ I =

1 2

(

∑4

ˆ φˆ l , where φˆ l is the Lagrange

l=1 ul

) ∂ 2 uˆ ∂ 2 uˆ . ξ ( ξ − 1) + η ( η − 1) ∂ξ 2 ∂η2

(2) If uˆ ∈ P3 (Kˆ ), then uˆ − uˆ I = ξ (ξ − 1)

(

1 ∂ 3 uˆ 6 ∂ξ 3

where C1 , C2 are constants: C1 =

ξ+

1 ∂ 3 uˆ

)

2 ∂ξ 2 ∂η

1 ∂ 2 uˆ 1 2 ∂ξ 2 ξ = 3 ,η=0

|

η + C1 + η(η − 1)

, C2 =

(

1 ∂ 3 uˆ 6 ∂η3

η+

1 ∂ 3 uˆ 2 ∂ξ ∂η2

)

ξ + C2 ,

1 ∂ 2 uˆ 1. 2 ∂η2 ξ =0,η= 3

|

Lemma 3.2 (See [26]). Assume that a(x) ∈ C 1 (Ω ), u ∈ H 3 (Ω ) ∩ H01 (Ω ), and uI is the interpolation function defined by (7). Then, for any (x, y) ∈ K , we have

⏐∫ ⏐ ⏐ ∂ (u − uI ) ⏐⏐ 3 ⏐ ⏐ −→ (a(x, y) − a(xO , yO )) ∂ x dy⏐ ≲ h |u|3,K , l = 1, 3, OMl ⏐∫ ⏐ ⏐ ∂ (u − uI ) ⏐⏐ 3 ⏐ (a(x , y) − a(x , y )) dx O O ⏐ −→ ⏐ ≲ h |u|3,K , l = 2, 4, ∂ y OMl where (xO , yO ) is the center of K (see Fig. 2(a)).

(13) (14)

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335

327

ˆ l , 1 ≤ l ≤ 4. Fig. 2. (a) K and Dl , 1 ≤ l ≤ 4. (b) Kˆ and D

We will consider the estimate for A(u − uI , v ) (v ∈ Vh ) in the following. For ∀v ∈ Vh , we firstly derive the asymptotic expansion of A(u − uI , v )|K , where uI ≜ Π u. From (9), we have

∫ ∂ (u − uI ) ∂ (u − uI ) a ds + (v3 − v4 ) a ds A(u, v )|K = (v2 − v1 ) ∂ n ∂ n3 1 OM1 OM3 ∫ ∫ ∂ (u − uI ) ∂ (u − uI ) + (v 4 − v 1 ) a a ds + (v3 − v2 ) ds ∂ n4 ∂ n2 OM4 OM2 ∫

=

4 ∑

Ii ,

(15)

i=1

where n1 , n3 and n2 , n4 are, respectively, the normal directions of OM 1 , OM 3 , and OM 2 , OM 4 . Without loss of generality, we only discuss the first term I1 in (15). Let Eu ≜ u − uI and n1 =

1 l1

(β1 , −α1 ),

x +x −x −x

y +y −y −y

where α1 = 3 4 2 1 2 , β1 = 3 4 2 1 2 , l1 = By some simple calculations, we have

Jk

α12 + β12 .

|P1 P4 | + |P2 P3 |

= h2 (1 + O(h2 )), ds = l1 dη, ) 2 ( ) 1 1 , η = S + 2(S143 − S142 ) η − .

l1 =

(

√

2

(16) (17)

2

Furthermore, we also can obtain

β1

α1

)⏐ ⏐ l1 ⏐ = 1 , l1 l1 ⏐ξ = 1 Jk ( 2 , η ) 2 ( )⏐ [ ( )] β1 α1 ⏐⏐ −1 1 2 2 2 2 ηx − ηy = |P1 P3 | − |P2 P4 | + (|P2 P3 | − |P1 P4 | ) η − . l1 l1 ⏐ 1 2 l J ( 1 , η)

(

ξx

− ξy

ξ= 2

I1 = ( v 2 − v 1 )

∫ a OM 1

= (v 2 − v 1 )

∫

∂ Eu ds ∂ n1 (

a (Eu)x OM 1

= (v 2 − v 1 )

1 k 2

∫ 0

1 2

β1 l1

− (Eu)y

α1 l1

) ds

( ) ( )) ⏐ ⏐ β1 α1 β1 α1 ⏐ a (Eu)ξ ξx − ξy + (Eu)η ηx − ηy ⏐ (

= (v2 − v1 )(I11 + I12 ),

l1

l1

l1

l1

ξ = 21

dη

(18)

(19)

328

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335

where 1 2

∫

( )⏐ β1 α1 ⏐⏐ a(Eu)ξ ξx − ξy ⏐

I11 = 0 1 2

∫

(

I12 =

a(Eu)η 0

l1

l1

β1

α1

dη,

ξ = 12

)⏐ ⏐ ⏐ ηx − ηy dη. l1 l1 ⏐ξ = 1 2

We will estimate I11 , I12 , respectively. From Lemma 2.1, (16), (17) and (18), one can see that I11 = l21

1 2

∫

a(Eu)ξ Jk−1 (ξ , η)|ξ = 1 dη 2

0

= =

h2 h1 h2 h1

∫

1 2

a(Eu)ξ |ξ = 1 dη + O(h ) 2

2

0

ˆ a(O)

1 2

∫

1 2

∫

a(Eu)ξ |ξ = 1 dη 2

0 1 2

∫

h1

2

0

+ O(h2 )

h2

(Eu)ξ |ξ = 1 dη +

1 2

∫

ˆ (a(ξ , η) − a(O))(Eu) ξ | ξ = 1 dη 2

0

a(Eu)ξ |ξ = 1 dη, 2

0

ˆ = ( 1 , 1 ). where O 2 2 For the second term above, according to Lemma 3.2, one can get h2

1 2

∫

h1

0

3 ˆ (a(ξ , η) − a(O))(Eu) ξ |ξ = 1 dη ≲ h |u|3,K . 2

For the third term above, according to Lemma 3.1 and Bramble–Hilbert lemma, one can obtain

⏐∫ ⏐ ⏐ ⏐

1 2

0

⏐ ⏐ a(Eu)ξ |ξ = 1 dη⏐⏐ ≤ |u|2,Kˆ ≤ Ch|u|2,K . 2

Hence I11

ˆ h2 = a(O)

1 2

∫

h1

0

(Eu)ξ |ξ = 1 dη + O(h3 )∥u∥3,K .

(20)

2

From Lemma 2.1, (16), (17), (19) and Bramble–Hilbert lemma, one can get the estimate on I12 as follows.

|I12 | ≤ Ch2 |u|2,Kˆ ≤ Ch3 |u|2,K .

(21)

Combining (20) with (21), we can obtain

ˆ I1 = (v2 − v1 )a(O)

h2

1 2

∫

h1

(Eu)ξ |ξ = 1 dη + O(h3 )∥u∥3,K |v|1,K . 2

0

By similar induction,

ˆ h2 I2 = (v3 − v4 )a(O)

1

∫

h1

1 2

(Eu)ξ |ξ = 1 dη + O(h3 )∥u∥3,K |v|1,K . 2

Therefore,

ˆ I1 + I2 = (v2 − v1 + v3 − v4 )a(O)

∑4

Due to v ∈ P1,1 (Kˆ ) and v = simple calculations, we can obtain

∫ Kˆ

h2 h1

1 2

∫ 0

(Eu)ξ |ξ = 1 dη + O(h3 )∥u∥3,K |v|1,K .

(22)

2

v φˆ l , where φˆ 1 = (1 − ξ )(1 − η), φˆ 2 = ξ (1 − η), φˆ 3 = ξ η, φˆ 4 = (1 − ξ )η, by some

l=1 l

∂v 1 dξ dη = [(v1 − v2 ) + (v3 − v4 )]. ∂ξ 2

(23)

When u ∈ Π3 , form Lemma 3.1, we have 1 2

∫ 0

(Eu)ξ |ξ = 1 dη = − 2

1 ∂ 3u 48 ∂ ξ 3

−

∂ 3u . 24 ∂ξ ∂ 2 η 1

(24)

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335

329

By (23) and (24), the first term of the right side in (22) leads to

ˆ (v2 − v1 + v3 − v4 )a(O)

∫

h2 h1

1 2

(Eu)ξ |ξ = 1 dη = − 2

0

1 24

ˆ a(O)

h2 h1

∫ gˆ1 Kˆ

∂v dξ dη, ∂ξ

where

∂ 3u ∂ 3u + 2 . ∂ 3ξ ∂ξ ∂ 2 η

gˆ1 =

Furthermore, when u ∈ H 4 (Ω ) ∩ H01 (Ω ), by Bramble–Hilbert lemma,

∫ 1 2 h2 ˆ (v2 − v1 + v3 − v4 )a(O) (Eu)ξ |ξ = 1 dη 2 h1 0 ∫ 1 h2 ∂v ˆ = − a(O) gˆ1 dξ dη + O(h3 )∥u∥4,Kˆ |v|1,Kˆ 24 h1 Kˆ ∂ξ ∫ ∫ 1 h2 ∂v 1 h2 ˆ − a)gˆ1 ∂v dξ dη + O(h3 )∥u∥ ˆ |v| ˆ =− agˆ1 dξ dη − (a(O) 4,K 1 ,K 24 h1 Kˆ ∂ξ 24 h1 Kˆ ∂ξ ∫ 1 h2 ∂v =− agˆ1 dξ dη + O(h3 )∥u∥4,Kˆ |v|1,Kˆ . 24 h1 Kˆ ∂ξ Hence, (22) results in I1 + I2 = −

∫

1 h2

agˆ1

24 h1

Kˆ

∂v dξ dη + O(h3 )∥u∥4,Kˆ |v|1,Kˆ . ∂ξ

(25)

Analogous derivation leads to I3 + I4 = −

∫

1 h1

agˆ2

24 h2

Kˆ

∂v dξ dη + O(h3 )∥u∥4,kˆ |v|1,Kˆ , ∂η

(26)

where

∂ 3u ∂ 3u + 2 . ∂ 3η ∂η∂ 2 ξ

gˆ2 =

Substituting (25) and (26) into (15), we have A(u − uI , v )|K = −

1

[

24

h1

∫ agˆ2

h2

Kˆ

∂v h2 dξ dη + ∂η h1

∫ agˆ1 Kˆ

] ∂v dξ dη + O(h3 )∥u∥4,Kˆ |v|1,Kˆ . ∂ξ

From Lemma 2.1, we have

−

] ∫ ∂v ∂v h2 agˆ1 dξ dη + dξ dη 24 h2 Kˆ ∂η h1 Kˆ ∂ξ ) ( ) )] ∫ [( ( 3 h1 ∂ u ∂ 3u ∂v h2 ∂ 3 u ∂ 3u ∂v 1 +2 + +2 dξ dη =− a 24 Kˆ h2 ∂ 3 η ∂η∂ 2 ξ ∂η h1 ∂ 3 ξ ∂ξ ∂ 2 η ∂ξ ) ( ) )] ∫ [(( 3 3 ∂v ∂ 3u ∂v ∂ 3u 1 1∂ u −1 ∂ u = − h1 h2 a λ− + 2 λ + λ + 2 λ dxdy + O(h3 )∥u∥4,Kˆ |v|1,Kˆ , h h 3 h h 24 ∂ 3y ∂ y∂ 2 x ∂ y ∂ x ∂ x∂ 2 y ∂ x K 1

[

where λh =

h1

∫

agˆ2

h1 . h2

Hence, (15) leads to A(u − uI , v )|K

=−

1 24

[((

∫ h1 h2

a

−1 ∂

∂ 3u + 2 λ h ∂ 3y ∂ y∂ 2 x

λh

K

3

u

)

( ) )] 3 ∂v ∂ 3u ∂v −1 ∂ u + λh 3 + 2λh dxdy + O(h3 )∥u∥4,K |v|1,K . ∂y ∂ x ∂ x∂ 2 y ∂ x

Lemma 3.3. Assume that u is the solution of (1) and u ∈ H 4 (Ω ) ∩ H01 (Ω ), then A(u − uI , v ) =

1 24

∫ h1 h2

Ω

g Π ∗ v dxdy + O(h3 )∥u∥4 |v|1 ,

where A(u, v ) is defined by (9), and

( ( )) ( ( )) 3 3 ∂ ∂ 3u ∂ ∂ 3u −1 ∂ u −1 ∂ u a λh 3 + 2λh + a λh 3 + 2λh . g = ∂y ∂ y ∂ y∂ 2 x ∂x ∂ x ∂ x∂ 2 y

(27)

330

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335

Proof. Taking the sum for K , from Cauchy inequality and Green formula, we have A(u − uI , v ) =

M ∑

A(u − uI , v )|K

k=1

=−

1 24

h1 h2

M [∫ ∑

M ∑

a

−1 ∂

λh

K

k=1

+ O(h3 )

[((

∂ 3u + 2 λ h ∂ 3y ∂ y∂ 2 x 3

u

)

( ) )] ] 3 ∂v ∂ 3u ∂v −1 ∂ u + λh 3 + 2λh dxdy ∂y ∂ x ∂ x∂ 2 y ∂ x

∥u∥4,K |v|1,K

k=1

=−

1 24

[((

∫ h1 h2

a Ω

1 λ− h

∂ 3u ∂ 3u + 2λh 3 ∂ y ∂ y∂ 2 x

)

( ) )] 3 ∂v ∂ 3u ∂v 1 ∂ u + λh 3 + 2λ− dxdy h ∂y ∂ x ∂ x∂ 2 y ∂ x

+ O(h3 )∥u∥4 |v|1 )) ( ( ))] ∫ [ ( ( 3 3 ∂ ∂ 3u ∂ 3u 1 ∂ 1∂ u −1 ∂ u = h1 h2 a λ− + 2 λ + a λ + 2 λ v dxdy h h 3 h h 24 ∂ 3y ∂ y∂ 2 x ∂x ∂ x ∂ x∂ 2 y Ω ∂y ∫ 1 g v dxdy + O(h3 )∥u∥4 |v|1 , + O(h3 )∥u∥4 |v|1 = h1 h2 24

Ω

where g is defined as (27). Furthermore, we have A(u − uI , v ) =

=

1 24 1 24

∫ h1 h2

∫Ω

h1 h2

Ω

g Π ∗ v dxdy +

1 24

∫ h1 h2

Ω

g(v − Π ∗ v )dxdy + O(h3 )∥u∥4 |v|1

g Π v dxdy + O(h )∥u∥4 |v|1 . ∗

3

This completes the proof of this lemma. □ Now, we shall come to derive the asymptotic expansion. The discrete Green function ghz ∈ Vh satisfies A(vh , ghz ) = vh (z), vh ∈ Vh .

(28)

Lemma 3.4 (See [19,20]). The discrete Green function ghz ∈ Vh in (28) satisfies that 1

|ghz |1 ≲ |lnh| 2 , ghz ∈ Vh .

(29)

Lemma 3.5. If w ∈ H 2 (Ω ), then for any ϵ ∈ (0, 1), we have

∥w − wI ∥0,∞ ≲ h1−ϵ ∥w∥2 , where wI is defined by (7). If u ∈ H 4 (Ω ) ∩ H01 (Ω ), then we can introduce an auxiliary problem

{ −∇ · (a∇w) = g , x ∈ Ω , w|∂ Ω = 0,

(30)

where g is defined by (27). From the regularity of solution for problem (30), we have

∥w∥2 ≲ ∥g ∥0 ≲ |u|4 .

(31)

The weak solution of problem (30) w ∈ H 2 (Ω ) ∩ H01 (Ω ) satisfies A(w, v ) = (g , Π ∗ v ), ∀ v ∈ H01 (Ω ).

(32)

Assume that wh ∈ Vh is the isoparametric bilinear finite volume element solution of problem (30). Then we have A(wh , v h ) = (g , Π ∗ v h ), ∀v h ∈ Vh .

(33)

Setting v = v ∈ Vh in (32), and subtracting (33) from (32), we have h

A(w − wh , v h ) = 0, ∀v h ∈ Vh .

(34)

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335

331

Lemma 3.6. If u ∈ H 4 (Ω ) ∩ H01 (Ω ), then for any ϵ ∈ (0, 1), we have

∥w − wh ∥0,∞ ≲ O(h1−ϵ )|u|4 . From Lemma 3.3 and (33), the following result holds. Lemma 3.7. Assume that u ∈ H 4 (Ω ) ∩ H01 (Ω ) is the solution of (1), uh , wh ∈ Vh , respectively, are solutions of (8) and (33), and uI ∈ Vh is the interpolation function of u. Then we have

( A uh − uI −

h1 h2 24

wh , v

)

= O(h3 )∥u∥4 |v|1 , ∀ v ∈ Vh .

(35)

Theorem 3.1. Under the hypotheses of Lemma 3.7, we have h1 h2

w(z) + O(h3−ϵ )∥u∥4 , 24 where z is an inner node and ϵ ∈ (0, 1) is an arbitrary constant independent of the mesh size h. (uh − uI )(z) =

(36)

Proof. According to Lemma 3.7, we have

( A uh − uI −

h1 h2 24

wh , v

)

= O(h3 )∥u∥4 |v|1 , ∀ v ∈ Vh .

(37)

By the definition of the discrete Green function (28) and setting vh = ghz ∈ Vh in (37), and from Lemma 3.4, we have

( uh − uI −

h1 h2 24

) ( ) h1 h2 wh (z) = A uh − uI − wh , ghz 24

3

z 4 gh 1

= O(h )∥u∥ | | 1

= O(h3 |lnh| 2 )∥u∥4 . From Lemma 3.5, and noticing the fact that ϵ ∈ (0, 1) is an arbitrary constant independent of the mesh size h, which leads to 1

|lnh| 2 = O(h−ϵ ), we have (uh − uI )(z) =

h1 h2 24 h1 h2

1

(wh − w + w )(z) + O(h3 |lnh| 2 )∥u∥4

w(z) + O(h3−ϵ )∥u∥4 . 24 This completes the proof of Theorem 3.1. □ =

4. Superconvergence Let ulh = uh (Pl ) (l = 1, 2, 3, 4) and define h

∂ x u(P) :=

u2h − u1h

h

, ∂ y u(P) :=

u4h − u3h

, (38) 2h 2h to be the average partial derivative values at point P along x1 and x2 directions, respectively. Then we define the discrete average gradient operator

(( ( h ) )2 ( h )2 ) 21 h h ∇ h u(P) = ∂ x u, ∂ y u (P) and |∇ h u(P)| = ∂ x u(P) + ∂ y u(P) .

(39)

Similar to the proof of Lemma 3.5 or that in [27], we can prove that the following result holds. Lemma 4.1. Assume that u ∈ H 4 (Ω ) ∩ H01 (Ω ) is the solution of (1), then for any inner node P(xi , yi ), we have

|∇ h (u − uI )(P)| = O(h2−ϵ )∥u∥4 .

(40)

Lemma 4.2. Assume that u ∈ H 4 (Ω ) ∩ H01 (Ω ) is the solution of (1) and uh ∈ Vh is the isoparametric bilinear finite volume element solution. Then we have

|∇ h (uh − uI )(P)| = O(h2−ϵ )∥u∥4 , where ∇ h (uh − uI )(P) = (∂

h x (uh

− uI )(P), ∂

(41) h y (uh

− uI )(P)).

332

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335

From Lemmas 4.1 and 4.2, and noticing that

|(∇ h u − ∇ u)(P)| ≲ h2 |u|4 , one can obtain the following superconvergence result. Theorem 4.1. Assume that u ∈ H 4 (Ω ) ∩ H01 (Ω ) is the solution of problem (1) and uh ∈ Vh is solution of (8). Then, for any inner node P and any constant ϵ ∈ (0, 1), we have

|(∇ h uh − ∇ u)(P)| = O(h2−ϵ )∥u∥4 .

(42)

Remark 4.1. The result in Theorem 4.1 holds for the following parabolic problem (similar to that in the literature [25])

⎧ ut − ∇ · (a(x)∇ u) = f (x, t), x ∈ Ω , 0 < t ≤ T , ⎪ ⎪ ⎪u = α (x, t), ⎨ x ∈ ∂ Ω1 , 0 < t ≤ T , ∂u = β (x, t), x ∈ ∂ Ω2 , 0 < t ≤ T , ⎪ ⎪ ⎪ ⎩ ∂n u = u0 (x), x ∈ Ω , t = 0,

(43)

where u0 (x) is a given function, and other functions are similar to Problem (1). This fact can be confirmed by Example 5.3 in the following section. 5. Numerical experiment In this section, we will display three examples to verify theoretical results. The first two are for elliptic problems with homogeneous Dirichlet boundary condition and Dirichlet–Neumann boundary condition, and with the diffusion coefficients a(x, y) = 1 + x + y and 1 + x3 + y3 , respectively. The third one is for a linear parabolic problem with Dirichlet–Neumann boundary condition and a(x, y) = 1 + x3 + y3 . In these examples, Ω = (0, 1)2 . Let Ω h be a disturbed uniform quadrilateral partition, and N1 and N2 be, the partition numbers along x and y axis directions, respectively. In the following, we choose N1 = N2 = N. P(xi , yj ) = (i h1 , j h2 ), ˜ ˜ , y˜ ) = is disturbed by (∆x, ∆y) = (mod(i, 2)h21 , mod(j, 2)h22 ), where hi = N1 , i = 1, 2. Hence, the disturbed point P(x P P (xi + ∆x, yj + ∆y). It is obvious that the above disturbance is stronger than that in (2). The scheme (8) is employed. To verify the theoretical results, we carry on some numerical tests on two hands. On one hand, we choose and force two typical inner points for above partition to be static to obtain the exact (static) derivatives, and display some results about them P˜ 1 (0.75, 0.25), P˜ 2 (0.625, 0.875) in the second and fourth columns of all the following tables, where eh1 = |(∂x u − h h ∂ x uh )(P˜ i )|, eh2 = |(∂y u − ∂ y uh )(P˜ i )|, (1 ≤ i ≤ 2) are, respectively, the errors of partial derivatives about variables x and y at point P˜ i . On the other hand, we also carry on some numerical experiments on non-static points on the above disturbed parti-

tion (there is not any forced point), and choose the maximum error over all inner points to confirm the results in Theorem 4.1. h The results are shown in the sixth column of all the tables, where eh1,max = maxk∈DI |(∂x u−∂ x uh )(P˜ k )|, eh2,max = maxk∈DI |(∂y u− h

h ∂ y uh )(P˜ k )|, respectively, and DI is the set of all inner nodes on the disturbed partition. γ is the ratio of ehi and ei2 .

Example 5.1. We consider problem (1) and take a(x) = 1 + x + y, u(x, y) = sin(π x) sin(π y), ∂ Ω1 = ∂ Ω , ∂ Ω2 = Φ , and α (x) = 0. In this example, one can see that Problem (1) satisfies the homogeneous Dirichlet boundary condition from the solution function u(x, y). Numerical results are shown in Tables 1 and 2. From above tables, for the homogeneous Dirichlet boundary case, one can see that the approximations are convergent of order two about the partial derivatives for these two typical points. The maximum errors for the derivative along x and y directions over all inner points are also of order two. They confirm the results in Theorem 4.1. Example 5.2. We consider problem (1) and take a(x) = 1 + x3 + y3 , u(x, y) = cos(π x) sin(π y) + 2.0,

∂ Ω1 = {(x, y)|(0, 1) × {y = 0, y = 1}}, ∂ Ω2 = {(x, y)|{x = 0, x = 1} × (0, 1)}, α (x) = cos(π x) sin(π y) + 2.0, β (x) = sin(π x) sin(π y). In this example, one can see that Problem (1) satisfies the Dirichlet and Neumann boundary conditions from the solution function u(x, y). Numerical results are shown in Tables 3 and 4. From above tables, for the Dirichlet–Neumann boundary case, one can see that the approximations are all convergent of order two and similar to Example 5.1. It validates that the results in Theorem 4.1 hold for other type of boundary conditions.

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335

333

Table 1 The results for Dirichlet boundary condition case along x direction. N

16 32 64 128

P˜ 1 (0.75,0.25)

P˜ 2 (0.625,0.875)

eh1

γ

eh1

γ

1.35e−3 3.36e−4 8.38e−5 2.09e−5

4.07 3.99 4.01

2.48e−4 6.04e−5 1.50e−5 3.74e−6

4.01 4.00 4.00

eh1,max

γ

4.56e−2 1.13e−2 2.81e−3 7.03e−4

4.04 4.02 4.00

Table 2 The results for Dirichlet boundary condition case along y direction. N

16 32 64 128

P˜ 1 (0.75,0.25) eh2 6.89e−4 1.63e−4 4.03e−5 1.01e−5

P˜ 2 (0.625,0.875)

γ

eh2

γ

4.04 4.01 4.00

2.14e−3 5.42e−4 1.35e−4 3.37e−5

3.91 3.98 4.01

eh2,max

γ

4.55e−2 1.13e−2 2.82e−3 7.03e−4

4.02 4.01 4.01

Table 3 The results for Dirichlet–Neumann boundary condition case along x direction. P˜ 1 (0.75,0.25)

N

16 32 64 128

eh1 8.97e−3 2.23e−3 5.58e−4 1.39e−4

P˜ 2 (0.625,0.875)

γ

eh1

γ

4.02 4.00 4.01

2.87e−3 7.13e−4 1.78e−4 4.45e−5

4.02 4.00 4.00

eh1,max

γ

6.36e−2 1.63e−2 4.11e−3 1.03e−3

3.90 3.96 3.99

Table 4 The results for Dirichlet–Neumann boundary condition case along y direction. P˜ 1 (0.75,0.25)

N

16 32 64 128

eh2 6.89e−4 1.63e−4 4.03e−5 1.01e−5

P˜ 2 (0.625,0.875)

γ

eh2

γ

4.82 4.04 4.01 4.00

2.14e−3 5.42e−4 1.35e−4 3.37e−5

3.83 3.91 3.98 4.01

eh2,max

γ

5.15e−2 1.34e−2 3.41e−3 8.61e−4

3.84 3.93 3.96

Table 5 The results for x direction at t = 0.5. N

P˜ 1 (0.75,0.25) eh1

16 32 64 128

2.15e−2 5.26e−3 1.27e−3 3.11e−4

P˜ 2 (0.625,0.875)

γ

eh1

γ

4.07 4.14 4.08

8.53e−3 2.00e−3 5.03e−4 1.27e−4

4.26 3.97 3.96

eh1,max

γ

3.89e−2 9.92e−3 2.49e−3 6.21e−4

3.92 3.98 4.00

Example 5.3. We consider problem (43) and take a(x) = 1 + x3 + y3 , u(x, t) = e−t cos(π x) sin(π y) + 2.0, 0 ≤ t ≤ 1,

∂ Ω1 = {(x, y)|(0, 1) × {y = 0, y = 1}}, ∂ Ω2 = {(x, y)|{x = 0, x = 1} × (0, 1)}, α (x, t) = e−t cos(π x) sin(π y) + 2.0, β (x, t) = e−t sin(π x) sin(π y). In this example, for Problem (43), we choose T = 1 and take the partition for the time direction as follows ti =

i Nt

, i = 1, 2, . . . , Nt ,

where Nt is the number of partitions along t direction. The backward Euler method is used, and ∆t are chosen as the square of the space step size h2 to be convenient to confirm the convergent order O(c1 ∆t + c2 h2 ) = O(h2 ) for this parabolic problem. Numerical results we choose two typical time t = 0.5 and t = 1 are shown in Tables from 5 to 8, respectively. From the results in Tables 5–8, one can see that the approximations are convergent of order two about the partial derivatives for any inner points. Hence the results in Theorem 4.1 can be extended to parabolic problems.

334

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335 Table 6 The results for y direction at t = 0.5. N

16 32 64 128

P˜ 1 (0.75,0.25)

P˜ 2 (0.625,0.875)

eh2

γ

eh2

γ

1.82e−2 4.71e−3 1.21e−3 3.11e−4

3.86 3.89 3.89

1.99e−2 4.89e−3 1.24e−3 3.17-e4

4.06 3.94 3.91

eh2,max

γ

3.22e−2 8.46e−3 2.27e−3 5.82e−4

3.81 3.73 3.90

eh1,max

γ

2.36e−2 6.03e−3 1.54e−3 3.78e−4

3.91 3.92 4.07

eh2,max

γ

1.94e−2 5.04e−3 1.29e−3 3.25e−4

3.85 3.91 3.97

Table 7 The results for x direction at t = 1.0. N

16 32 64 128

P˜ 1 (0.75,0.25) eh1 1.31e−2 3.22e−3 7.92e−4 1.96e−4

P˜ 2 (0.625,0.875)

γ

eh1

γ

4.07 4.06 4.04

5.17e−3 1.21e−3 2.98e−4 7.46e−5

4.26 4.06 3.99

Table 8 The results for y direction at t = 1.0. N

16 32 64 128

P˜ 1 (0.75,0.25) eh2 1.13e−2 2.76e−3 6.97e−4 1.75e−4

P˜ 2 (0.625,0.875)

γ

eh2

γ

4.08 3.96 3.98

1.19e−2 2.85e−3 7.04−e4 1.77−e4

4.17 4.05 3.98

Acknowledgments The authors would like to thank the referees for the helpful and valuable suggestions to improve this paper greatly. References [1] R.H. Li, P.Q. Zhu, Generalized difference methods for second order elliptic partial differential equations quadrilateral grids II, Numer. Math. J. Chinese Univ. 4 (1982) 360–375. [2] T. Schmidt, Box schemes on quadrilateral meshes, Computing 51 (1993) 271–292. [3] T.A. Porsching, Error estimates for MAC-like approximations to the linear Navier–Stokes equations, Numer. Math. 29 (1987) 291–364. [4] S.H. Chou, Analysis and convergence of a covolume element method for the generalized Stokes problem, Math. Comp. 66 (1997) 85–104. [5] H.X. Rui, Symmetric modified finite volume element methods for self-adjoint elliptic and parabolic problems, J. Comput. Appl. Math. 146 (2002) 373–386. [6] X. Ma, S. Shu, A. Zhou, Symmetric finite volume discretization for parabolic problems, Comput. Meth. Appl. Mech. Eng. 192 (2003) 4467–4485. [7] Z.Q. Cai, J. Douglas, M. Park, Development and analysis of higher order finite volume methods over rectangles for elliptic equations, Adv. Comput. Math. 19 (2003) 3–33. [8] M. Yang, Cubic finite volume methods for second order elliptic equations with variable coefficients, Northeast. Math. J. 21 (2005) 146–152. [9] L. Chen, A new class of high order finite volume methods for second order elliptic equations, SIAM. J. Numer. Anal. 47 (2010) 4021–4043. [10] T.K. Wang, Alternating direction finite volume element methods for three-dimensional parabolic equations, Numer. Math. Theory Methods Appl. 3 (2010) 499–522. [11] A.A Abedini, R.A. Ghiassi, Three-dimensional finite volume model for shallow water flow simulation, Aust. J. Basic An. 4 (2010) 3208–3215. [12] Y.P. Lin, J.G. Liu, M. Yang, Finite volume element methods :an overview on recent developments, Int. J. Numer. Anal. Model. Ser. B 4 (1) (2013) 14–34. [13] X.M. He, T. Lin, Y. Lin, A bilinear immersed finite volume element method for the diffusion equation with discontinuous coefficient, Commun. Comput. Phys. 6 (1) (2013) 185–202. [14] L. Zh, Zhi, Y. Zh, Zh.L. Li, An immersed finite volume element method for 2D PDEs with discontinuous coefficients and non-homogeneous jump conditions, Comput. Math. Anal. 70 (2015) 89–103. [15] Z.Y. Chen, R.H. Li, A.H. Zhou, A note on the optimal L2-estimate of the finite volume element method, Adv. Comput. Math. 16 (2002) 291–303. [16] Y.H. Li, R.H. Li, Generalized difference methods on arbitrary quadrilateral networks, J. Comput. Math. 6 (1999) 653–672. [17] J.L. Lv, Y.H. Li, em L2 error estimate of the finite volume element methods on quadrilateral meshes, Adv. Comput. Math. 33 (2010) 129–148. [18] J.L. Lv, L2 Error Estimates and Superconvergence of the Finite Volume Element Methods on Quadrilateral Meshes, (Doctor thesis), JiLin University, 2009 (in Chinese). [19] C.M. Chen, Y.Q. Huang, High Accuracy of Finite Element Method, Hunan Science Press, China, 1995 (in chinese). [20] C.M. Chen, Structure Theory of Superconvergence of Finite Elements, Hunan Science Press, China, 2001 (in chinese). [21] Y.Q. Huang, H.F. Qin, D.S. Wang, Centroidal Voronoi Tessellation-based finite element superconvergence, Int. J. Numer. Methods Engrg. 76 (2008) 1819–1839. [22] E. Süli, Convergence of finite volume schemes for Poisson’s equation on nonuniform meshes, Int. J. Numer. Anal. Model. 3 (2006) 348–360. [23] S. Shu, H.Y. Yu, Y.Q. Huang, C.Y. Nie, A preserving-symmetry finite volume scheme and superconvergence on quadrangle grids, Int. J. Numer. Anal. Model. 3 (2006) 348–360.

C.Y. Nie et al. / Journal of Computational and Applied Mathematics 321 (2017) 323–335

335

[24] C.Y. Nie, S. Shu, H.Y. Yu, J. Wu, Superconvergence and asymptotic expansions for bilinear finite volume element approximations, Numer. Math. Theory Methods Appl. 6 (2013) 408–423. [25] C.Y. Nie, S. Shu, H.Y. Yu, Y.Y. Yang, Superconvergence and asymptotic expansion for semidiscrete bilinear finite volume element approximation of the parabolic problem, Comput. Math. Appl. 66 (2013) 91–104. [26] C.Y. Nie, Several Finite Volume Element Schemes and Some Applications in Radiation Heat Conduction Problems, (Doctor thesis), Xiangtan university, 2010 (in Chinese). [27] L. Zhang, On convergence of isoparametric bilinear finite elements, Comm. Numer. Meth. Eng. 12 (1996) 849–862.