Table of Contents Author Guidelines Submit a Manuscript
International Journal of Mathematics and Mathematical Sciences
VolumeΒ 2012, Article IDΒ 870402, 21 pages
http://dx.doi.org/10.1155/2012/870402
Research Article

A Second Order Characteristic Method for Approximating Incompressible Miscible Displacement in Porous Media

School of Mathematics, Shandong University, Jinan 250100, China

Received 25 March 2012; Revised 5 July 2012; Accepted 28 July 2012

Academic Editor: NistorΒ Victor

Copyright Β© 2012 Tongjun Sun and Keying Ma. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

An approximation scheme is defined for incompressible miscible displacement in porous media. This scheme is constructed by two methods. Under the regularity assumption for the pressure, cubic Hermite finite element method is used for the pressure equation, which ensures the approximation of the velocity smooth enough. A second order characteristic finite element method is presented to handle the material derivative term of the concentration equation. It is of second order accuracy in time increment, symmetric, and unconditionally stable. The optimal 𝐿 2 -norm error estimates are derived for the scalar concentration.

1. Introduction

In this paper, we will consider the following incompressible miscible displacement in porous media, which is governed by a coupled system of partial differential equations with initial and boundary values. The pressure is governed by an elliptic equation and the concentration is governed by a convection-diffusion equation [16] as follows: ξ‚΅ π‘˜ βˆ’ βˆ‡ β‹… πœ‡ ξ‚Ά πœ™ βˆ‡ 𝑝 ≑ βˆ‡ β‹… 𝑒 = π‘ž , π‘₯ ∈ Ξ© , 𝑑 ∈ 𝐽 , ( 1 . 1 a ) πœ• 𝑐 πœ• 𝑑 + 𝑒 β‹… βˆ‡ 𝑐 βˆ’ βˆ‡ β‹… ( 𝐷 βˆ‡ 𝑐 ) = ( Μƒ 𝑐 βˆ’ 𝑐 ) Μƒ π‘ž , π‘₯ ∈ Ξ© , 𝑑 ∈ 𝐽 , ( 1 . 1 b ) 𝑒 β‹… 𝜈 = ( 𝐷 ( π‘₯ ) βˆ‡ 𝑐 ) β‹… 𝜈 = 0 , π‘₯ ∈ πœ• Ξ© , 𝑑 ∈ 𝐽 , ( 1 . 1 c ) 𝑐 ( π‘₯ , 0 ) = 𝑐 0 ( π‘₯ ) , π‘₯ ∈ Ξ© , ( 1 . 1 d ) where Ξ© is a bounded domain in β„› 2 , 𝐽 = ( 0 , 𝑇 ] , and Μƒ π‘ž = m a x { π‘ž , 0 } is nonzero at injection wells only. The variables in (1.1a)–(1.1d) are the pressure 𝑝 ( π‘₯ , 𝑑 ) in the fluid mixture, the Darcy velocity 𝑒 = ( 𝑒 1 , 𝑒 2 ) ξ…ž , and the relative concentration 𝑐 ( π‘₯ , 𝑑 ) of the injected fluid. The 𝜈 is the unit outward normal vector on boundary πœ• Ξ© .

The coefficients and data in (1.1a)–(1.1d) are π‘˜ ( π‘₯ ) , the permeability of the porous media; πœ‡ ( π‘₯ ) , the viscosity of the fluid mixture; π‘ž ( π‘₯ , 𝑑 ) , representing flow rates at wells; πœ™ ( π‘₯ ) , the porosity of the rock; Μƒ 𝑐 ( π‘₯ , 𝑑 ) , the injected concentration at injection wells ( π‘ž > 0 ) and the resident concentration at production wells ( π‘ž < 0 ); Μƒ π‘ž = m a x ( π‘ž , 0 ) . Here, for the diffusion coefficient, we consider a dispersion-free case [5, 6] as follows: 𝐷 = πœ™ 𝑑 π‘š 𝐼 , ( 1 . 2 ) where 𝑑 π‘š is the molecular diffusivity and 𝐼 is a 2 Γ— 2 identity matrix. Furthermore, a compatibility condition ∫ Ξ© π‘ž ( π‘₯ , 𝑑 ) 𝑑 π‘₯ = 0 must be imposed to determine the pressure.

The pressure equation is elliptic and easily handled, but the concentration equation is parabolic and normally convection dominated. It is well known that standard Galerkin scheme applied to the convection-dominated problems does not work well and produces excessive numerical diffusion or nonphysical oscillation. The characteristic method has been introduced to obtain better approximations for (1.1a)–(1.1d), such as characteristic finite element method [37], characteristic finite difference method [8], the modified of characteristic finite element method (MMOC-Galerkin) [9], and the Eulerian-Lagrangian localized adjoint method (ELLAM) [10].

We had considered a combined numerical approximation for (1.1) in [11]. Standard mixed finite element was used for Darcy velocity equation and a characteristics-mixed finite element method was presented for approximating the concentration equation. Characteristic approximation was applied to handle the convection term, and a lowest order mixed finite element spatial approximation was adopted to deal with the diffusion term. Thus, the scalar unknown concentration and the diffusive flux can be approximated simultaneously. This approximation conserves mass globally. The optimal 𝐿 2 -norm error estimates were derived.

It should be pointed out that the works mentioned above all gave one order accuracy in time increment Ξ” 𝑑 . That is to say, the first order characteristic method in time was analyzed. As for higher order characteristic method in time, Rui and Tabata [12] used the second order Runge-Kutta method to approximate the material derivative term for convection-diffusion problems. The scheme presented was of second order accuracy in time increment Ξ” 𝑑 , symmetric, and unconditionally stable. Optimal error estimates were proved in the framework of 𝐿 2 -theory. Numerical analysis of convection-diffusion-reaction problems with higher order characteristic/finite elements were analyzed in [13, 14], which extended the work [12]. The 𝑙 ∞ ( 𝐿 2 ) error estimates of second order in time increment Ξ” 𝑑 were obtained.

The goal of this paper is to present a second order characteristic finite element method in time increment to handle the material derivative term of the concentration equation of (1.1a)–(1.1d). It is organized as follows. In Section 2, we formulate the second order characteristic finite element method for the concentration and cubic Hermite finite element method for the pressure, respectively. Then, we present a combined approximation scheme. In Section 3, we analyze the stability of the approximation scheme for the concentration equation. In Section 4, we derive the optimal-order 𝐿 2 -norm error estimates for the scalar concentration. They are of second order accuracy in time increment, symmetric, and unconditionally stable. We conclude our results in Section 5.

2. Formulation of the Method

2.1. Statements and Assumptions

In this paper, we adopt notations and norms of usual Sobolev spaces. For periodic functions, we use the notation as in [4] as follows. Let: π‘Š π‘˜ 𝑝 = ξ‚‹ π‘Š π‘˜ 𝑝 ξ‚» πœ• ( Ξ© ) = πœ“ ∢ 𝛼 πœ“ πœ• π‘₯ 𝛼 ∈ 𝐿 𝑝 ξ‚Ό ( Ξ© ) f o r | 𝛼 | ≀ π‘˜ , p e r i o d i c ( 2 . 1 ) be the periodic Sobolev space on Ξ© with the usual norm. If 𝑝 = 2 , we write 𝐻 π‘˜ =  𝐻 π‘˜ ξ‚‹ π‘Š ( Ξ© ) = π‘˜ 2 ( Ξ© ) ( 2 . 2 ) with norm β€– πœ“ β€– π‘˜  𝐻 = β€– πœ“ β€– π‘˜ ( Ξ© ) , β€– πœ“ β€– = β€– πœ“ β€– 0 = β€– πœ“ β€– 𝐿 2 ( Ξ© ) . ( 2 . 3 )

Moreover, we adopt some notations for the functional spaces involved, which were used in [1214]. For a Banach space 𝑋 and a positive integer π‘š , spaces 𝐢 π‘š ( [ 0 , 𝑇 ] , 𝑋 ) and 𝐻 π‘š ( ( 0 , 𝑇 ) , 𝑋 ) are abbreviated as 𝐢 π‘š ( 𝑋 ) and 𝐻 π‘š ( 𝑋 ) , respectively, and endowed with norms β€– πœ‘ β€– 𝐢 π‘š ( 𝑋 ) ∢ = m a x [ ] 𝑑 ∈ 0 , 𝑇 ξ‚» m a x 𝑗 = 0 , … , π‘š β€– β€– πœ‘ ( 𝑗 ) β€– β€– ( 𝑑 ) 𝑋 ξ‚Ό , β€– πœ‘ β€– 𝐻 π‘š ( 𝑋 )  ξ€œ ∢ = 𝑇 0 π‘š  𝑗 = 0 β€– β€– πœ‘ ( 𝑗 ) β€– β€– ( 𝑑 ) 2 𝑋 ξƒͺ 𝑑 𝑑 1 / 2 , ( 2 . 4 ) where πœ‘ ( 𝑗 ) denotes the 𝑗 th derivative of πœ‘ with respect to time. The Banach space 𝑍 π‘š is defined by 𝑍 π‘š = ξ€½ 𝑓 ∈ 𝐢 𝑗 ξ€· 𝐻 π‘š βˆ’ 𝑗 ξ€Έ ξ€Ύ ( Ξ© ) ; 𝑗 = 0 , … , π‘š , ( 2 . 5 ) equipped with the norm β€– πœ‘ β€– 𝑍 π‘š ξ€½ ∢ = m a x β€– πœ‘ β€– 𝐢 𝑗 ( 𝐻 π‘š βˆ’ 𝑗 ) ξ€Ύ . ; 0 ≀ 𝑗 ≀ π‘š ( 2 . 6 )

We also require the following assumptions on the coefficients in (1.1) [3]. Let π‘Ž βˆ— , π‘Ž βˆ— , πœ™ βˆ— , πœ™ βˆ— , and 𝐾 βˆ— be positive constants such that 0 < π‘Ž βˆ— ≀ π‘˜ ( π‘₯ ) πœ‡ ( π‘₯ ) ≀ π‘Ž βˆ— , 0 < πœ™ βˆ— ≀ πœ™ ( π‘₯ ) ≀ πœ™ βˆ— , | | | | + | | | Μƒ π‘ž ( π‘₯ , 𝑑 ) πœ• Μƒ π‘ž | | | πœ• 𝑑 ( π‘₯ , 𝑑 ) ≀ 𝐾 βˆ— . ( C ) Other assumptions will be made in individual theorems as necessary.

2.2. A Cubic Hermite Element Method for the Pressure

The variational form for the pressure equation (1.1a) is equal to the following: ξ‚΅ π‘˜ πœ‡ ξ‚Ά βˆ‡ 𝑝 , βˆ‡ πœ’ = ( π‘ž ( 𝑑 ) , πœ’ ) , πœ’ ∈ 𝐻 1 ( Ξ© ) , 𝑑 ∈ 𝐽 . ( 2 . 7 )

For β„Ž 𝑝 > 0 , we discretize (2.7) in space on a quasiuniform triangle mesh 𝒯 β„Ž 𝑝 of Ξ© with diameter of element ≀ β„Ž 𝑝 . Let π‘Š β„Ž 𝑝 βŠ‚ 𝐿 2 ( Ξ© ) be a cubic Hermite finite element space for this mesh. The finite element method for the pressure, given at a time 𝑑 ∈ 𝐽 , consists of 𝑃 ∈ π‘Š β„Ž 𝑝 such that ξ‚΅ π‘˜ πœ‡ ξ‚Ά βˆ‡ 𝑃 , βˆ‡ πœ’ = ( π‘ž ( 𝑑 ) , πœ’ ) , βˆ€ πœ’ ∈ π‘Š β„Ž 𝑝 . ( 2 . 8 )

Since the left-hand side of (2.7) ((2.8), resp.) is a continuous and coercive bilinear form and the right-hand side of (2.7) ((2.8), resp.) is a continuous linear functional, existence and uniqueness of 𝑝 ( 𝑃 , resp.) are ensured obviously [15].

By the theory of the cubic Hermit element, the finite element space π‘Š β„Ž 𝑝 possess the following approximation property [15, 16]: i n f πœ’ ∈ π‘Š β„Ž 𝑝 β€– 𝑔 βˆ’ πœ’ β€– 𝑠 ≀ 𝐾 β„Ž 𝑝 4 βˆ’ 𝑠 β€– 𝑔 β€– 4 , 𝑠 = 0 , 1 , βˆ€ 𝑔 ∈ 𝐻 4 ( Ξ© ) , ( 2 . 9 ) where 𝐾 is a positive constant independent of β„Ž 𝑝 .

Define Π: 𝐢 2 ( Ξ© ) β†’ π‘Š β„Ž 𝑝 is the interplant operator. Then, we have the following theorem.

Theorem 2.1. Let 𝑝 and 𝑃 be the solutions of (1.1a)–(1.1d) and (2.8), respectively, and assume 𝑝 ∈ 𝐻 4 ( Ξ© ) . Then, there exists a positive constant 𝐾 independent of β„Ž 𝑝 such that β€– βˆ‡ ( 𝑝 βˆ’ 𝑃 ) β€– ≀ 𝐾 β„Ž 3 𝑝 β€– 𝑝 β€– 4 . ( 2 . 1 0 )

Proof. By Céa lemma, we have π‘Ž β€– β€– βˆ‡ ( 𝑝 βˆ’ 𝑃 ) 2 ≀ ξ‚΅ π‘˜ πœ‡ ξ‚Ά = ξ‚΅ π‘˜ βˆ‡ ( 𝑝 βˆ’ 𝑃 ) , βˆ‡ ( 𝑝 βˆ’ 𝑃 ) πœ‡ ξ‚Ά + ξ‚΅ π‘˜ βˆ‡ ( 𝑝 βˆ’ 𝑃 ) , βˆ‡ ( 𝑝 βˆ’ Ξ  𝑝 ) πœ‡ ξ‚Ά = ξ‚΅ π‘˜ βˆ‡ ( 𝑝 βˆ’ 𝑃 ) , βˆ‡ ( Ξ  𝑝 βˆ’ 𝑃 ) πœ‡ ξ‚Ά βˆ‡ ( 𝑝 βˆ’ 𝑃 ) , βˆ‡ ( 𝑝 βˆ’ Ξ  𝑝 ) ≀ π‘Ž βˆ— β€– βˆ‡ ( 𝑝 βˆ’ 𝑃 ) β€– β€– βˆ‡ ( 𝑝 βˆ’ Ξ  𝑃 ) β€– . ( 2 . 1 1 ) Then by (2.9), we can derive β€– βˆ‡ ( 𝑝 βˆ’ 𝑃 ) β€– ≀ 𝐾 β„Ž 3 𝑝 β€– 𝑝 β€– 4 . ( 2 . 1 2 )

2.3. A Second Order Characteristic Method for the Concentration

Now, we define the characteristic lines associated with vector field 𝑒 and recall some classical properties satisfied by them. Thus, for given ( π‘₯ , 𝑑 ) ∈ Ξ© Γ— [ 0 , 𝑇 ] , the characteristic line through ( π‘₯ , 𝑑 ) is the vector function 𝑋 𝑒 ( π‘₯ , 𝑑 ; β‹… ) solving the following initial value problem: πœ• 𝑋 𝑒 ξ€· 𝑋 πœ• 𝜏 ( π‘₯ , 𝑑 ; 𝜏 ) = 𝑣 𝑒 ξ€Έ ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 , 𝑋 𝑒 ( π‘₯ , 𝑑 ; 𝜏 ) = π‘₯ , ( 2 . 1 3 ) where 𝑣 ( 𝑋 𝑒 ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 ) ∢ = 𝑒 ( 𝑋 𝑒 ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 ) / πœ™ .

Next, assuming they exist, we denote by 𝐹 𝑒 (by 𝐿 , resp.) the gradient of 𝑋 𝑒 (of 𝑣 , resp.) with respect to the space variable π‘₯ , that is, ξ€· 𝐹 𝑒 ξ€Έ π‘Ÿ 𝑠 ( πœ• ξ€· 𝑋 π‘₯ , 𝑑 ; 𝜏 ) ∢ = 𝑒 ξ€Έ π‘Ÿ πœ• π‘₯ 𝑠 ( π‘₯ , 𝑑 ; 𝜏 ) , 𝐿 π‘Ÿ 𝑠 ( π‘₯ , 𝑑 ) ∢ = πœ• 𝑣 π‘Ÿ πœ• π‘₯ 𝑠 ( π‘₯ , 𝑑 ) . ( 2 . 1 4 ) We adopted some propositions and lemmas from [13].

Proposition 2.2. If 𝑣 ∈ 𝐢 0 ( 𝐢 𝑛 ( Ξ© ) ) for 𝑛 β‰₯ 1 an integer, then 𝑋 𝑒 ∈ 𝐢 0 ( Ξ© Γ— [ 0 , 𝑇 ] Γ— [ 0 , 𝑇 ] ) and it is 𝐢 𝑛 with respect to the π‘₯ variable.

In order to compute second order approximations of matrices 𝐹 𝑒 and 𝐹 𝑒 βˆ’ 1 , we need the following equations: πœ• 𝐹 𝑒 ξ€· 𝑋 πœ• 𝜏 ( π‘₯ , 𝑑 ; 𝜏 ) = 𝐿 𝑒 ξ€Έ 𝐹 ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 𝑒 πœ• ( π‘₯ , 𝑑 ; 𝜏 ) , 2 𝐹 𝑒 πœ• 𝜏 2 ξ‚€ ( π‘₯ , 𝑑 ; 𝜏 ) = βˆ‡ πœ• 𝑣  ξ€· 𝑋 πœ• 𝑑 + 𝐿 𝑣 𝑒 ξ€Έ 𝐹 ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 𝑒 ( π‘₯ , 𝑑 ; 𝜏 ) . ( 2 . 1 5 )

Proposition 2.3. If 𝑣 ∈ 𝐢 0 ( 𝐢 1 ( Ξ© ) ) , then β€– β€– 𝐹 𝑒 ( β€– β€– π‘₯ , 𝑑 ; 𝜏 ) ≀ 𝑒 β€– 𝑣 β€– 𝐢 0 1 ( ( 𝐢 Ξ© ) ) | 𝜏 βˆ’ 𝑑 | [ ] . , βˆ€ π‘₯ ∈ Ξ© , 𝑑 , 𝜏 ∈ 0 , 𝑇 ( 2 . 1 6 )

Proposition 2.4. If 𝑣 ∈ 𝐢 0 ( 𝐢 2 ( Ξ© ) ) ∩ 𝐢 1 ( 𝐢 1 ( Ξ© ) ) , then 𝐹 𝑒 satisfies the Taylor expansions as 𝐹 𝑒 + ξ€œ ( π‘₯ , 𝑑 ; 𝑠 ) = 𝐼 + ( 𝑠 βˆ’ 𝑑 ) 𝐿 ( π‘₯ , 𝑑 ) 𝑑 𝑠 ( ξ‚€ 𝜏 βˆ’ 𝑠 ) βˆ‡ πœ• 𝑣  ξ€· 𝑋 πœ• 𝑑 + 𝐿 𝑣 𝑒 ( ξ€Έ 𝐹 π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 𝑒 ( π‘₯ , 𝑑 ; 𝜏 ) 𝑑 𝜏 , ( 2 . 1 7 ) and its inverse, 𝐹 𝑒 βˆ’ 1 , satisfies the Liouville theorem as 𝐹 𝑒 βˆ’ 1 ξ€· 𝑋 ( π‘₯ , 𝑑 ; 𝑠 ) = 𝐼 + ( 𝑠 βˆ’ 𝑑 ) 𝐿 𝑒 ξ€Έ βˆ’ ξ€œ ( π‘₯ , 𝑑 ; 𝑠 ) , 𝑠 𝑑 𝑠 ξ‚€ ( 𝜏 βˆ’ 𝑑 ) βˆ‡ πœ• 𝑣  ξ€· 𝑋 πœ• 𝑑 + 𝐿 𝑣 𝑒 ξ€Έ 𝐹 ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 𝑒 ξ€· 𝑋 𝑒 ξ€Έ ( π‘₯ , 𝑑 ; 𝑠 ) , 𝑠 ; 𝜏 𝑑 𝜏 . ( 2 . 1 8 )

By using Liouville’s theorem and the chain rule, we obtain πœ• πœ• 𝜏 d e t 𝐹 𝑒 βˆ’ 1 ( π‘₯ , 𝑑 ; 𝜏 ) = βˆ’ d e t 𝐹 𝑒 βˆ’ 1 ξ€· 𝑋 ( π‘₯ , 𝑑 ; 𝜏 ) d i v 𝑣 𝑒 ξ€Έ , πœ• ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 2 πœ• 𝜏 2 d e t 𝐹 𝑒 βˆ’ 1 ( π‘₯ , 𝑑 ; 𝜏 ) = βˆ’ d e t 𝐹 𝑒 βˆ’ 1 ξ€· ( π‘₯ , 𝑑 ; 𝜏 ) ( d i v 𝑣 ) 2 ξ€· 𝑋 𝑒 ξ€Έ ξ‚€ ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 Γ— d i v πœ• 𝑣  ξ€· 𝑋 πœ• 𝑑 + 𝐿 𝑣 𝑒 ξ€Έ βˆ’ ξ€· ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 𝐿 β‹… 𝐿 𝑇 𝑋 ξ€Έ ξ€· 𝑒 . ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 ξ€Έ ξ€Έ ( 2 . 1 9 )

Proposition 2.5. If 𝑣 ∈ 𝐢 0 ( 𝐢 1 ( Ξ© ) ) , then d e t 𝐹 𝑒 ( π‘₯ , 𝑑 ; 𝜏 ) ≀ 𝑒 β€– 𝑣 β€– 𝐢 0 1 ( ( 𝐢 Ξ© ) ) | 𝜏 βˆ’ 𝑑 | [ ] . , βˆ€ π‘₯ ∈ Ξ© , 𝑑 , 𝜏 ∈ 0 , 𝑇 ( 2 . 2 0 )

Proposition 2.6. If 𝑣 ∈ 𝐢 0 ( 𝐢 2 ( Ξ© ) ) ∩ 𝐢 1 ( 𝐢 1 ( Ξ© ) ) , then d e t 𝐹 𝑒 βˆ’ 1 satisfies d e t 𝐹 𝑒 βˆ’ 1 ( ξ€œ π‘₯ , 𝑑 ; 𝑠 ) = 1 βˆ’ ( 𝑠 βˆ’ 𝑑 ) d i v 𝑣 ( π‘₯ , 𝑑 ) + 𝑑 𝑠 ( πœ• 𝜏 βˆ’ 𝑠 ) 2 πœ• 𝜏 2 ξ€· d e t 𝐹 𝑒 βˆ’ 1 ξ€Έ ( π‘₯ , 𝑑 ; 𝜏 ) 𝑑 𝜏 . ( 2 . 2 1 )

Variational Formulation
From the definition of the characteristic curves and by using the chain rule, it follows that πœ™ 𝑑 𝑐 ξ€· 𝑋 𝑑 𝜏 𝑒 ξ€Έ ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 = πœ™ πœ• 𝑐 ξ€· 𝑋 πœ• 𝑑 𝑒 ξ€Έ ξ€· 𝑋 ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 + 𝑒 𝑒 ξ€Έ ξ€· 𝑋 ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 β‹… βˆ‡ 𝑐 𝑒 ξ€Έ . ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 ( 2 . 2 2 ) By writing (1.1b) at point 𝑋 𝑒 ( π‘₯ , 𝑑 ; 𝜏 ) and time 𝜏 and using (2.22), we have πœ™ 𝑑 𝑐 ξ€· 𝑋 𝑑 𝜏 𝑒 ξ€Έ ξ€· 𝑋 ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 βˆ’ βˆ‡ β‹… ( 𝐷 βˆ‡ 𝑐 ) 𝑒 ξ€Έ ξ€· 𝑋 ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 = ( ( Μƒ 𝑐 βˆ’ 𝑐 ) Μƒ π‘ž ) 𝑒 ξ€Έ . ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 ( 2 . 2 3 )
Before giving a week formulation of (2.23), we adopted a lemma from [13], which can be considered as Green’s formula.

Lemma 2.7. Let 𝑋 ∢ Ξ© β†’ 𝑋 ( Ξ© ) , 𝑋 ∈ 𝐢 2 ( Ξ© ) be an invertible vector valued function. Let 𝐹 = βˆ‡ 𝑋 and assume that 𝐹 βˆ’ 1 ∈ 𝐢 1 ( Ξ© ) . Then ξ€œ Ξ© ξ€œ d i v 𝐰 ( 𝑋 ( π‘₯ ) ) πœ“ ( π‘₯ ) 𝑑 π‘₯ = Ξ“ 𝐹 βˆ’ 𝑇 𝑛 ( π‘₯ ) β‹… 𝐰 ( 𝑋 ( π‘₯ ) ) πœ“ ( π‘₯ ) 𝑑 𝐀 π‘₯ βˆ’ ξ€œ Ξ© 𝐹 βˆ’ 1 ξ€œ 𝐰 ( 𝑋 ( π‘₯ ) ) β‹… βˆ‡ πœ“ ( π‘₯ ) 𝑑 π‘₯ βˆ’ Ξ© d i v 𝐹 βˆ’ 𝑇 β‹… 𝐰 ( 𝑋 ( π‘₯ ) ) πœ“ ( π‘₯ ) 𝑑 π‘₯ , ( 2 . 2 4 ) with 𝐰 ∈ 𝐻 1 ( 𝑋 ( Ξ© ) ) a vector valued function and πœ“ ∈ 𝐻 1 ( Ξ© ) a scalar function.

Now, we can multiply (2.23) by a test function πœ“ ∈ 𝐻 1 ( Ξ© ) , integrate in Ξ© , and apply the usual Green’s formula and (2.24) with 𝑋 ( π‘₯ ) = 𝑋 𝑒 ( π‘₯ , 𝑑 ; 𝜏 ) , obtaining ξ€œ Ξ© πœ™ 𝑑 𝑐 ξ€· 𝑋 𝑑 𝜏 𝑒 ( ξ€Έ ξ€œ π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 πœ“ ( π‘₯ ) 𝑑 π‘₯ + Ξ© 𝐹 𝑒 βˆ’ 1 ( ξ€· 𝑋 π‘₯ , 𝑑 ; 𝜏 ) 𝐷 βˆ‡ 𝑐 𝑒 ( ξ€Έ + ξ€œ π‘₯ , 𝑑 ; 𝜏 ) β‹… βˆ‡ πœ“ ( π‘₯ ) 𝑑 π‘₯ Ξ© d i v 𝐹 𝑒 βˆ’ 𝑇 ξ€· 𝑋 ( π‘₯ , 𝑑 ; 𝜏 ) β‹… 𝐷 βˆ‡ 𝑐 𝑒 ξ€Έ = ξ€œ ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 πœ“ ( π‘₯ ) 𝑑 π‘₯ Ξ© ξ€· 𝑋 ( ( Μƒ 𝑐 βˆ’ 𝑐 ) Μƒ π‘ž ) 𝑒 ξ€Έ ( π‘₯ , 𝑑 ; 𝜏 ) , 𝜏 πœ“ ( π‘₯ ) 𝑑 π‘₯ . ( 2 . 2 5 )

2.4. The Combined Approximation Scheme

We now present our sequential time-stepping procedure that combines (2.8) and (2.25). Part 𝐽 into 0 = 𝑑 0 < 𝑑 1 < β‹― < 𝑑 𝑁 = 𝑇 , with Ξ” 𝑑 𝑛 𝑐 = 𝑑 𝑛 βˆ’ 𝑑 𝑛 βˆ’ 1 . The analysis is valid for variable time steps, but we drop the superscript from Ξ” 𝑑 𝑐 for convenience. For functions 𝑓 on Ξ© Γ— 𝐽 , we write 𝑓 𝑛 ( π‘₯ ) for 𝑓 ( π‘₯ , 𝑑 𝑛 ) . As in [3], let us part 𝐽 into pressure time steps 0 = 𝑑 0 < 𝑑 1 < β‹― < 𝑑 𝑀 = 𝑇 , with Ξ” 𝑑 π‘š 𝑝 = 𝑑 π‘š βˆ’ 𝑑 π‘š βˆ’ 1 . Each pressure step is also a concentration step, that is, for each π‘š there exists 𝑛 such that 𝑑 π‘š = 𝑑 𝑛 , in general, Ξ” 𝑑 𝑝 > Ξ” 𝑑 𝑐 . We may vary Ξ” 𝑑 𝑝 , but except for Ξ” 𝑑 1 𝑝 we drop the superscript. For functions 𝑓 on Ξ© Γ— 𝐽 , we write 𝑓 π‘š ( π‘₯ ) for 𝑓 ( π‘₯ , 𝑑 π‘š ) ; thus, subscripts refer to pressure steps and superscripts to concentration steps.

If concentration step 𝑑 𝑛 relates to pressure steps by 𝑑 π‘š βˆ’ 1 < 𝑑 𝑛 ≀ 𝑑 π‘š , we require a velocity approximation for (2.25) based on π‘ˆ π‘š βˆ’ 1 and earlier values. If π‘š β‰₯ 2 , take the linear extrapolation of π‘ˆ π‘š βˆ’ 1 and π‘ˆ π‘š βˆ’ 2 defined by 𝐸 π‘ˆ 𝑛 = ξ‚΅ 𝑑 1 + 𝑛 βˆ’ 𝑑 π‘š βˆ’ 1 𝑑 π‘š βˆ’ 1 βˆ’ 𝑑 π‘š βˆ’ 2 ξ‚Ά π‘ˆ π‘š βˆ’ 1 βˆ’ 𝑑 𝑛 βˆ’ 𝑑 π‘š βˆ’ 1 𝑑 π‘š βˆ’ 1 βˆ’ 𝑑 π‘š βˆ’ 2 π‘ˆ π‘š βˆ’ 2 ; ( 2 . 2 6 ) if π‘š = 1 , set 𝐸 π‘ˆ 𝑛 = π‘ˆ 0 . ( 2 . 2 7 ) We retain the superscript on Ξ” 𝑑 1 𝑝 because 𝐸 π‘ˆ 𝑛 is first-order correct in time during the first pressure step and second-order during later steps.

Let 𝑋 ∢ ( 0 , 𝑇 ) β†’ 𝑅 2 be a solution of the ordinary differential equation 𝑑 𝑋 = Μ‚ Μ‚ 𝑑 𝑑 𝑣 ( 𝑋 , 𝑑 ) , w h e r e 𝑣 ( 𝑋 , 𝑑 ) ∢ = 𝐸 π‘ˆ ( 𝑋 , 𝑑 ) πœ™ . ( 2 . 2 8 ) Then, we can write πœ™ 𝑑 𝑑 𝑑 𝑐 ( 𝑋 , 𝑑 ) = πœ™ πœ• 𝑐 πœ• 𝑑 + 𝐸 π‘ˆ ( 𝑋 , 𝑑 ) β‹… βˆ‡ 𝑐 . ( 2 . 2 9 ) Subject to an initial condition 𝑋 ( 𝑑 𝑛 + 1 ) = π‘₯ , we get approximate values of 𝑋 at 𝑑 𝑛 by the Euler method and the second order Runga-Kutta method, respectively, 𝑋 𝑛 𝐸 ( π‘₯ ) = π‘₯ βˆ’ Ξ” 𝑑 𝑐 Μ‚ 𝑣 𝑛 + 1 ( 𝑋 π‘₯ ) , 𝑛 𝑅 𝐾 ( π‘₯ ) = π‘₯ βˆ’ Ξ” 𝑑 𝑐 Μ‚ 𝑣 𝑛 + ( 1 / 2 ) ξ‚΅ π‘₯ βˆ’ Ξ” 𝑑 𝑐 2 Μ‚ 𝑣 𝑛 + 1 ξ‚Ά . ( π‘₯ ) ( 2 . 3 0 )

Next, assuming they exist, we denote by 𝐹 𝑛 𝐸 (resp., by 𝐿 𝐸 ) the gradient of 𝑋 𝑛 𝐸 (resp., of Μ‚ 𝑣 ( π‘₯ ) ) with respect to the space variable π‘₯ , that is, [13, 14] 𝐹 𝑛 𝐸 ( π‘₯ ) ∢ = βˆ‡ 𝑋 𝑛 𝐸 ( π‘₯ ) = 𝐼 ( π‘₯ ) βˆ’ Ξ” 𝑑 𝑐 𝐿 𝐸 𝑛 + 1 ( π‘₯ ) . ( 2 . 3 1 )

Hypothesis 2.8. For convenience, we assume that (1.1a)–(1.1d) is Ξ© -periodic ([3]), that is, all functions will be assumed to be spatially Ξ© -periodic throughout the rest of this paper.

This is physically reasonable, because no-flow conditions (1.1c) are generally treated by reflection, and, in general, interior flow patterns are much more important than boundary effects in reservoir simulation. Thus, the boundary conditions (1.1c) can be dropped.

Lemma 2.9. Under Hypothesis 2.8, if β€– Μ‚ 𝑣 β€– 𝐢 0 ( π‘Š 1 , ∞ ( Ξ© ) ) Ξ” 𝑑 𝑐 < 1 / 2 , we can see that 𝑋 𝑛 𝐸 ξ‚€ Ξ©  = 𝑋 𝑛 𝑅 𝐾 ξ‚€ Ξ©  = Ξ© . ( 2 . 3 2 )

Lemma 2.10. Under Hypothesis 2.8, if β€– Μ‚ 𝑣 β€– 𝐢 0 ( π‘Š 1 , ∞ ( Ξ© ) ) Ξ” 𝑑 𝑐 < 1 / 2 , we have ξ€· 𝐹 𝑛 𝐸 ξ€Έ βˆ’ 1 ( π‘₯ ) = 𝐼 + Ξ” 𝑑 𝑐 𝐿 𝐸 𝑛 + 1 ξ€· ( π‘₯ ) + Ξ” 𝑑 𝑐 ξ€Έ 2 ξ€· 𝐿 𝐸 𝑛 + 1 ξ€Έ ( π‘₯ ) 2 ξ‚€ ξ€· + 𝑂 Ξ” 𝑑 𝑐 ξ€Έ 3  . ( 2 . 3 3 )

Corollary 2.11. Under the assumptions of Lemma 2.10, f o r a l l π‘₯ ∈ Ξ© , we have ξ€· 𝐹 d e t 𝑛 𝐸 ξ€Έ βˆ’ 1 ( π‘₯ ) = 1 + Ξ” 𝑑 𝑐 Μ‚ 𝑣 d i v 𝑛 + 1 ξ‚€ ξ€· ( π‘₯ ) + 𝑂 Ξ” 𝑑 𝑐 ξ€Έ 2  , | | | ξ€· 𝐹 d e t 𝑛 𝐸 ξ€Έ βˆ’ 1 | | | ( π‘₯ ) ≀ 1 + Ξ” 𝑑 𝑐 β€– β€– Μ‚ 𝑣 𝑛 + 1 β€– β€– ( π‘₯ ) 𝐢 0 ( π‘Š 1 , ∞ ( Ξ© ) ) ξ‚€ ξ€· + 𝑂 Ξ” 𝑑 𝑐 ξ€Έ 2  . ( 2 . 3 4 )

The time difference for (2.25) will be combined with a standard finite element in the space variables. For β„Ž 𝑐 > 0 , we discrete (2.25) in space on a quasi-uniform mesh 𝒯 β„Ž 𝑐 of Ξ© with diameter of element ≀ β„Ž 𝑐 . Let 𝑀 β„Ž 𝑐 βŠ‚ 𝐻 1 ( Ξ© ) be a finite element space.

A characteristic discretization of the weak form (2.25) is given by { 𝐢 𝑛 } 𝑁 𝑛 = 1 ∈ 𝑀 β„Ž 𝑐 such that  πœ™ 𝐢 𝑛 + 1 βˆ’ 𝐢 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 Ξ” 𝑑 𝑐 ξƒͺ +  𝐷 , πœ‘ βˆ‡ 𝐢 𝑛 + 1 + ( βˆ‡ 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 2 ξƒͺ + , βˆ‡ πœ‘ Ξ” 𝑑 𝑐 2 ξ€· 𝐷 ( 𝐿 𝑛 βˆ‡ 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 ξ€Έ + , βˆ‡ πœ‘ Ξ” 𝑑 𝑐 2 Μ‚ 𝑣 ξ€· ξ€· βˆ‡ d i v 𝑛 β‹… 𝐷 βˆ‡ 𝐢 𝑛 ξ€Έ ∘ 𝑋 𝑛 𝐸 ξ€Έ + 1 , πœ‘ 2 ξ€· Μƒ π‘ž 𝑛 + 1 𝐢 𝑛 + 1 + ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 ξ€Έ = 1 , πœ‘ 2 ξ€· Μƒ π‘ž 𝑛 + 1 Μƒ 𝑐 𝑛 + 1 + ( Μƒ π‘ž 𝑛 Μƒ 𝑐 𝑛 ) ∘ 𝑋 𝑛 𝐸 ξ€Έ , πœ‘ , βˆ€ πœ‘ ∈ 𝑀 β„Ž 𝑐 , 𝐢 ( 2 . 3 5 a ) 0 =  𝐢 0 , βˆ€ π‘₯ ∈ Ξ© , ( 2 . 3 5 b ) where 𝑔 ∘ 𝑋 𝑛 𝐸 and 𝑔 ∘ 𝑋 𝑛 𝑅 𝐾 are compositions ξ€· 𝑔 ∘ 𝑋 𝑛 𝐸 ξ€Έ ξ€· 𝑋 ( π‘₯ ) = 𝑔 𝑛 𝐸 ξ€Έ , ξ€· ( π‘₯ ) 𝑔 ∘ 𝑋 𝑛 𝑅 𝐾 ξ€Έ ξ€· 𝑋 ( π‘₯ ) = 𝑔 𝑛 𝑅 𝐾 ξ€Έ ( π‘₯ ) , ( 2 . 3 6 )  𝐢 0 is an initial approximation of exact solution 𝑐 0 ( π‘₯ ) into 𝑀 β„Ž 𝑐 , which will be defined in Section 5.

At each pressure time step 𝑑 π‘š , we define π‘ˆ π‘š = π‘˜ πœ‡ βˆ‡ 𝑃 π‘š , 𝐢 ( 2 . 3 7 ) βˆ— = m i n ( m a x ( 𝐢 , 0 ) , 1 ) ( 2 . 3 8 ) is the truncation of 𝐢 to [ 0 , 1 ] . Then at 𝑑 π‘š , (2.8) is the following: ξ‚΅ π‘˜ πœ‡ βˆ‡ 𝑃 π‘š ξ‚Ά = ξ€· π‘ž , βˆ‡ πœ’ π‘š ξ€Έ , πœ’ , βˆ€ πœ’ ∈ π‘Š β„Ž 𝑝 ξ€· 𝑃 , ( 2 . 3 9 a ) π‘š ξ€Έ , 1 = 0 . ( 2 . 3 9 b ) The steps of calculation are as follows.

Step 1. 𝐢 0 known β†’ solve ( π‘ˆ 0 , 𝑃 0 ) by (2.37), (2.39a) and (2.39b);

Step 2. by (2.35a) and (2.35b) to solve 𝐢 1 β†’ then by (2.35a) and (2.35b) to solve 𝐢 2 ;

Step 3. analogously, { 𝐢 𝑗 βˆ’ 1 } 𝑛 1 𝑗 = 1 known β†’ { 𝐢 𝑗 } 𝑛 1 𝑗 = 1 such that 𝑑 𝑛 1 = 𝑑 1 ;

Step 4. then by (2.37), (2.39a) and (2.39b) for ( π‘ˆ 1 , 𝑃 1 ) ;

Step 5. calculate the approximations in turn analogously to get the pressure, velocity, and concentration at other time step, respectively.

Throughout the analysis, 𝐾 will denote a generic positive constant, independent of β„Ž 𝑐 , β„Ž 𝑝 , Ξ” 𝑑 𝑐 , and Ξ” 𝑑 𝑝 , but possibly depending on constants in ( C ) . Similarly, πœ€ will denote a generic small positive constant.

3. Stability for the Concentration Equation

In this section, we derive the stability for the concentration equation. For a given series of functions { πœ‘ } 𝑁 𝑛 = 0 , we define the following norms and seminorm: β€– πœ‘ β€– 𝑙 ∞ ( 𝐿 2 ) = m a x { β€– πœ‘ 𝑛 } , β€– ; 0 ≀ 𝑛 ≀ 𝑁 β€– πœ‘ 𝑛 β€– 𝑙 2 ( 𝐿 2 ) = ξƒ― Ξ” 𝑑 𝑐 𝑁  𝑛 = 0 β€– πœ‘ 𝑛 β€– 2 ξƒ° 1 / 2 , β€– πœ‘ β€– ξ…ž 𝑙 2 ξ€· 𝐻 1 ξ€Έ = ξƒ― Ξ” 𝑑 𝑐 𝑁 βˆ’ 1  𝑛 = 0 β€– β€– βˆ‡ πœ‘ 𝑛 + 1 + ( βˆ‡ πœ‘ 𝑛 ) ∘ 𝑋 𝑛 𝐸 β€– β€– 2 𝐷 ξƒ° 1 / 2 , β€– πœ‘ β€– 2 πœ™ = ( πœ™ πœ‘ , πœ‘ ) , β€– πœ‘ β€– 2 𝐷 = ( 𝐷 πœ‘ , πœ‘ ) , β€– πœ‘ β€– 2 𝐷 = ( π‘ž ( π‘₯ , 𝑑 ) πœ‘ , πœ‘ ) . ( 3 . 1 ) In the following sections, we use positive constants as 𝑐 0 = 𝑐 0 ξ€· β€– 𝐸 π‘ˆ β€– 𝐢 0 ( 𝐿 ∞ ) ξ€Έ , 𝑐 1 = 𝑐 1 ξ€· β€– 𝐸 π‘ˆ β€– 𝑙 ∞ ( π‘Š ( 1 , ∞ ) ) ξ€Έ , 𝑐 2 = 𝑐 2 ξ€· β€– 𝐸 π‘ˆ β€– 𝐢 0 ( π‘Š ( 2 , ∞ ) ) ∩ 𝐢 1 ( 𝐿 ∞ ) ξ€Έ . ( 3 . 2 ) In our analysis, we need some lemmas.

Lemma 3.1. Under the definitions (2.26) and (2.30), for 𝑛 = 0 , … , 𝑁 and πœ‘ ∈ 𝐿 2 ( Ξ© ) , it holds that β€– β€– πœ‘ ∘ 𝑋 𝑛 𝑖 β€– β€– 2 ≀ ξ€· 1 + 𝑐 1 Ξ” 𝑑 𝑐 ξ€Έ β€– πœ‘ β€– 2 , 𝑖 = 1 , 2 . ( 3 . 3 )

Proof. We only need to show a proof in the case 𝑖 = 1 . Let 𝐽 1 be the Jacobian matrix of the transformation 𝑦 = 𝑋 𝑛 𝐸 ( π‘₯ ) = π‘₯ βˆ’ 𝑣 𝑛 ( π‘₯ ) Ξ” 𝑑 𝑐 as 𝐽 1 = βŽ› ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 βˆ’ πœ• 𝑣 𝑛 1 ( π‘₯ ) πœ• π‘₯ 1 Ξ” 𝑑 𝑐 βˆ’ πœ• 𝑣 𝑛 1 ( π‘₯ ) πœ• π‘₯ 2 Ξ” 𝑑 𝑐 βˆ’ πœ• 𝑣 𝑛 2 ( π‘₯ ) πœ• π‘₯ 1 Ξ” 𝑑 𝑐 1 βˆ’ πœ• 𝑣 𝑛 2 ( π‘₯ ) πœ• π‘₯ 2 Ξ” 𝑑 𝑐 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . ( 3 . 4 ) According to the proof of (3.23) in [4], we have | | | | πœ• 𝑣 𝑛 𝑗 ( π‘₯ ) πœ• π‘₯ 𝑖 Ξ” 𝑑 𝑐 | | | | ≀ 𝐾 β„Ž βˆ’ 1 Ξ” 𝑑 𝑐 . ( 3 . 5 ) Since Ξ” 𝑑 𝑐 = π‘œ ( β„Ž 𝑐 ) , for β„Ž 𝑐 sufficiently small, we see that | | d e t 𝐽 1 | | βˆ’ 1 ≀ 𝑐 1 Ξ” 𝑑 𝑐 . ( 3 . 6 )
Following (3.6), we have β€– β€– πœ‘ ∘ 𝑋 𝑛 𝐸 β€– β€– 2 = ξ€œ Ξ© πœ‘ ξ€· 𝑋 𝑛 𝐸 ( ξ€Έ π‘₯ ) 2 ξ€œ 𝑑 π‘₯ = Ξ© πœ‘ ( 𝑦 ) 2 ξ€· d e t 𝐽 1 ξ€Έ βˆ’ 1 𝑑 𝑦 . ( 3 . 7 ) We complete the proof.

Suppose that { 𝐢 𝑛 } 𝑁 𝑛 = 0 βŠ‚ 𝐻 1 ( Ξ© ) and Μƒ 𝑐 ∈ 𝐢 0 ( 𝐿 2 ) be given. We define linear forms π’œ β„Ž 𝑛 + 1 / 2 and β„± β„Ž 𝑛 + 1 / 2 on 𝑀 β„Ž 𝑐 for 𝑛 = 0 , … , 𝑁 βˆ’ 1 by  π’œ β„Ž 𝑛 + 1 / 2  ≑  πœ™ 𝐢 𝐢 , πœ‘ 𝑛 + 1 βˆ’ 𝐢 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 Ξ” 𝑑 𝑐 ξƒͺ +  𝐷 , πœ‘ βˆ‡ 𝐢 𝑛 + 1 + ( βˆ‡ 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 2 ξƒͺ + , βˆ‡ πœ‘ Ξ” 𝑑 𝑐 2 ξ€· 𝐷 ( 𝐿 𝑛 βˆ‡ 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 ξ€Έ + , βˆ‡ πœ‘ Ξ” 𝑑 𝑐 2 Μ‚ 𝑣 ξ€· ξ€· βˆ‡ d i v 𝑛 β‹… 𝐷 βˆ‡ 𝐢 𝑛 ξ€Έ ∘ 𝑋 𝑛 𝐸 ξ€Έ + 1 , πœ‘ 2 ξ€· Μƒ π‘ž 𝑛 + 1 𝐢 𝑛 + 1 + ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 ξ€Έ ,  β„± , πœ‘ β„Ž 𝑛 + 1 / 2  ≑ 1 , πœ‘ 2 ξ€· Μƒ π‘ž 𝑛 + 1 Μƒ 𝑐 𝑛 + 1 + ( Μƒ π‘ž 𝑛 Μƒ 𝑐 𝑛 ) ∘ 𝑋 𝑛 𝐸 ξ€Έ , , πœ‘ ( 3 . 8 ) where πœ‘ ∈ 𝑀 β„Ž 𝑐 .

Lemma 3.2. Let { 𝐢 𝑛 } 𝑁 𝑛 = 0 be a solution of (2.35a) and (2.35b). Then it holds that  π’œ β„Ž 𝑛 + 1 / 2 𝐢 , 𝐢 𝑛 + 1  β‰₯ 𝐷 Ξ” 𝑑 𝑐 ξ‚΅ 1 2 β€– 𝐢 𝑛 β€– 2 πœ™ + Ξ” 𝑑 𝑐 4 β€– βˆ‡ 𝐢 𝑛 β€– 2 𝐷 + Ξ” 𝑑 𝑐 π‘ž βˆ— 4 β€– 𝐢 𝑛 β€– 2 ξ‚Ά + 1 2 Ξ” 𝑑 𝑐 β€– β€– 𝐢 𝑛 + 1 βˆ’ 𝐢 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 β€– β€– 2 πœ™ + 1 4 β€– β€– βˆ‡ 𝐢 𝑛 + 1 + βˆ‡ 𝐢 𝑛 ∘ 𝑋 𝑛 𝐸 β€– β€– 2 𝐷 + 1 4 π‘ž βˆ— β€– β€– Μƒ π‘ž 𝑛 + 1 𝐢 𝑛 + 1 + ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 β€– β€– 2 βˆ’ ξ‚» 𝑐 1 2 β€– 𝐢 𝑛 β€– 2 πœ™ + 𝑐 1 Ξ” 𝑑 𝑐 4 β€– βˆ‡ 𝐢 𝑛 β€– 2 𝐷 βˆ’ π‘ž βˆ— β€– 𝐢 𝑛 β€– 2 + 1 + 𝑐 1 Ξ” 𝑑 𝑐 π‘ž βˆ— β€– Μƒ π‘ž 𝑛 𝐢 𝑛 β€– 2 + 𝑐 1 Ξ” 𝑑 𝑐 2 ξ‚€ β€– β€– βˆ‡ 𝐢 𝑛 + 1 β€– β€– 2 𝐷 + β€– βˆ‡ 𝐢 𝑛 β€– 2 𝐷 + β€– β€– 𝐢 𝑛 + 1 β€– β€– 2 𝐷  ξ‚Ό , ( 3 . 9 ) where 𝐷 Ξ” 𝑑 𝑐 is the forward difference operator defined by 𝐷 Ξ” 𝑑 𝑐 πœ‘ 𝑛 = ξ€· πœ‘ 𝑛 + 1 βˆ’ πœ‘ 𝑛 ξ€Έ Ξ” 𝑑 𝑐 . ( 3 . 1 0 )

Proof. Substituting πœ‘ = 𝐢 𝑛 + 1 into (3.8), we have  π’œ β„Ž 𝑛 + 1 / 2 𝐢 , 𝐢 𝑛 + 1  ≑  πœ™ 𝐢 𝑛 + 1 βˆ’ 𝐢 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 Ξ” 𝑑 𝑐 , 𝐢 𝑛 + 1 ξƒͺ +  𝐷 βˆ‡ 𝐢 𝑛 + 1 + ( βˆ‡ 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 2 , βˆ‡ 𝐢 𝑛 + 1 ξƒͺ + Ξ” 𝑑 𝑐 2 ξ€· 𝐷 ( 𝐿 𝑛 βˆ‡ 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 , βˆ‡ 𝐢 𝑛 + 1 ξ€Έ + Ξ” 𝑑 𝑐 2 Μ‚ 𝑣 ξ€· ξ€· βˆ‡ d i v 𝑛 β‹… 𝐷 βˆ‡ 𝐢 𝑛 ξ€Έ ∘ 𝑋 𝑛 𝐸 , 𝐢 𝑛 + 1 ξ€Έ + 1 2 ξ€· Μƒ π‘ž 𝑛 + 1 𝐢 𝑛 + 1 + ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 , 𝐢 𝑛 + 1 ξ€Έ = 𝐼 1 + 𝐼 2 + 𝐼 3 + 𝐼 4 + 𝐼 5 . ( 3 . 1 1 ) Lemma 3.1 implies that 𝐼 1 =  πœ™ 𝐢 𝑛 + 1 βˆ’ 𝐢 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 Ξ” 𝑑 𝑐 , 𝐢 𝑛 + 1 + 𝐢 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 2 + 𝐢 𝑛 + 1 βˆ’ 𝐢 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 2 ξƒͺ β‰₯ 𝐷 Ξ” 𝑑 𝑐 ξ‚€ 1 2 β€– β€– 𝐢 𝑛 + 1 β€– β€– 2 πœ™  βˆ’ 𝑐 1 2 β€– 𝐢 𝑛 β€– 2 πœ™ + 1 2 Ξ” 𝑑 𝑐 β€– β€– 𝐢 𝑛 + 1 βˆ’ 𝐢 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 β€– β€– 2 πœ™ , 𝐼 ( 3 . 1 2 ) 2 =  𝐷 βˆ‡ 𝐢 𝑛 + 1 + ( βˆ‡ 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 2 , βˆ‡ 𝐢 𝑛 + 1 + ( βˆ‡ 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 2 + βˆ‡ 𝐢 𝑛 + 1 βˆ’ ( βˆ‡ 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 2 ξƒͺ β‰₯ 𝐷 Ξ” 𝑑 𝑐 ξ‚΅ Ξ” 𝑑 𝑐 4 β€– βˆ‡ 𝐢 𝑛 β€– 2 𝐷 ξ‚Ά βˆ’ 𝑐 1 Ξ” 𝑑 𝑐 4 β€– βˆ‡ 𝐢 𝑛 β€– 2 𝐷 + 1 4 β€– β€– βˆ‡ 𝐢 𝑛 + 1 + ( βˆ‡ 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 β€– β€– 2 𝐷 . ( 3 . 1 3 ) Next, by using 𝑐 1 Ξ” 𝑑 𝑐 < 1 , we obtain 𝐼 3 = Ξ” 𝑑 𝑐 2 ξ€· 𝐷 ( 𝐿 𝑛 βˆ‡ 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 , βˆ‡ 𝐢 𝑛 + 1 ξ€Έ ≀ Ξ” 𝑑 𝑐 2 √ 1 + 𝑐 1 Ξ” 𝑑 𝑐 β€– 𝐿 𝑛 βˆ‡ 𝐢 𝑛 β€– 𝐷 β€– β€– βˆ‡ 𝐢 𝑛 + 1 β€– β€– 𝐷 ≀ 𝑐 1 Ξ” 𝑑 𝑐 2  β€– βˆ‡ 𝐢 𝑛 β€– 2 𝐷 + β€– β€– βˆ‡ 𝐢 𝑛 + 1 β€– β€– 2 𝐷  . ( 3 . 1 4 ) Then when 𝐼 3 β‰₯ 0 and 𝐼 3 < 0 , we have 𝐼 3 𝑐 β‰₯ βˆ’ 1 Ξ” 𝑑 𝑐 2  β€– βˆ‡ 𝐢 𝑛 β€– 2 𝐷 + β€– β€– βˆ‡ 𝐢 𝑛 + 1 β€– β€– 2 𝐷  . ( 3 . 1 5 ) Similarly, for 𝐼 4 , we obtain the estimate 𝐼 4 𝑐 β‰₯ βˆ’ 1 Ξ” 𝑑 𝑐 2  β€– βˆ‡ 𝐢 𝑛 β€– 2 𝐷 + β€– β€– 𝐢 𝑛 + 1 β€– β€– 2 𝐷  . ( 3 . 1 6 ) Analogous computations to term 𝐼 2 give 𝐼 5 =  Μƒ π‘ž 𝑛 + 1 𝐢 𝑛 + 1 + ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 1 2 , Μƒ π‘ž 𝑛 + 1 𝐢 𝑛 + 1 + ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 1 2 Μƒ π‘ž 𝑛 + 1 + Μƒ π‘ž 𝑛 + 1 𝐢 𝑛 + 1 βˆ’ ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 1 2 Μƒ π‘ž 𝑛 + 1 ξƒͺ β‰₯ 1 4 ξ‚΅ π‘ž βˆ— β€– β€– 𝐢 𝑛 + 1 β€– β€– 2 βˆ’ 1 π‘ž βˆ— β€– β€– ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 1 β€– β€– 2 ξ‚Ά + 1 4 π‘ž βˆ— β€– β€– Μƒ π‘ž 𝑛 + 1 𝐢 𝑛 + 1 + ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 1 β€– β€– 2 β‰₯ 𝐷 Ξ” 𝑑 𝑐 ξ‚΅ Ξ” 𝑑 𝑐 π‘ž βˆ— 4 β€– 𝐢 𝑛 β€– 2 ξ‚Ά + π‘ž βˆ— 4 β€– 𝐢 𝑛 β€– 2 βˆ’ 1 + 𝑐 1 Ξ” 𝑑 𝑐 π‘ž βˆ— β€– Μƒ π‘ž 𝑛 𝐢 𝑛 β€– 2 + 1 4 π‘ž βˆ— β€– β€– Μƒ π‘ž 𝑛 + 1 𝐢 𝑛 + 1 + ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 β€– β€– 2 β‰₯ 𝐷 Ξ” 𝑑 𝑐 ξ‚΅ Ξ” 𝑑 𝑐 π‘ž βˆ— 4 β€– 𝐢 𝑛 β€– 2 ξ‚Ά + π‘ž βˆ— 4 β€– 𝐢 𝑛 β€– 2 βˆ’ 2 π‘ž βˆ— β€– Μƒ π‘ž 𝑛 𝐢 𝑛 β€– 2 + 1 4 π‘ž βˆ— β€– β€– Μƒ π‘ž 𝑛 + 1 𝐢 𝑛 + 1 + ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 β€– β€– 2 , ( 3 . 1 7 ) which completes the proof.

From Lemma 3.1, we have  β„± β„Ž 𝑛 + 1 / 2 , 𝐢 𝑛 + 1  = 1 2 ξ€· Μƒ π‘ž 𝑛 + 1 Μƒ 𝑐 𝑛 + 1 + ( Μƒ π‘ž 𝑛 Μƒ 𝑐 𝑛 ) ∘ 𝑋 𝑛 𝐸 , 𝐢 𝑛 + 1 ξ€Έ ≀ 1 β€– β€– 𝐢 2 πœ™ 𝑛 + 1 β€– β€– 2 πœ™ + 1 2  β€– β€– Μƒ π‘ž 𝑛 + 1 Μƒ 𝑐 𝑛 + 1 β€– β€– 2 + ξ€· 1 + 𝑐 1 Ξ” 𝑑 𝑐 ξ€Έ β€– Μƒ π‘ž 𝑛 Μƒ 𝑐 𝑛 β€– 2  . ( 3 . 1 8 ) Combining (3.18) with (3.8), we get 𝐷 Ξ” 𝑑 𝑐 ξ‚΅ 1 2 β€– 𝐢 𝑛 β€– 2 πœ™ + Ξ” 𝑑 𝑐 4 β€– βˆ‡ 𝐢 𝑛 β€– 2 𝐷 + Ξ” 𝑑 𝑐 π‘ž βˆ— 4 β€– 𝐢 𝑛 β€– 2 ξ‚Ά + 1 2 Ξ” 𝑑 𝑐 β€– β€– 𝐢 𝑛 + 1 βˆ’ 𝐢 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 β€– β€– 2 πœ™ + 1 4 β€– β€– βˆ‡ 𝐢 𝑛 + 1 + βˆ‡ 𝐢 𝑛 ∘ 𝑋 𝑛 𝐸 β€– β€– 2 𝐷 + 1 4 π‘ž βˆ— β€– β€– Μƒ π‘ž 𝑛 + 1 𝐢 𝑛 + 1 + ( Μƒ π‘ž 𝑛 𝐢 𝑛 ) ∘ 𝑋 𝑛 𝐸 β€– β€– 2 + π‘ž βˆ— β€– 𝐢 𝑛 β€– 2 ≀ 𝑐 1 2 β€– 𝐢 𝑛 β€– 2 πœ™ + 𝑐 1 Ξ” 𝑑 𝑐 4 β€– βˆ‡ 𝐢 𝑛 β€– 2 𝐷 + 2 π‘ž βˆ— β€– Μƒ π‘ž 𝑛 𝐢 𝑛 β€– 2 + 𝑐 1 Ξ” 𝑑 𝑐 2 ξ‚€ β€– β€– βˆ‡ 𝐢 𝑛 + 1 β€– β€– 2 𝐷 + β€– βˆ‡ 𝐢 𝑛 β€– 2 𝐷 + β€– β€– 𝐢 𝑛 + 1 β€– β€– 2 𝐷  + 1 β€– β€– 𝐢 2 πœ™ 𝑛 + 1 β€– β€– 2 πœ™ + 1 2  β€– β€– Μƒ π‘ž 𝑛 + 1 Μƒ 𝑐 𝑛 + 1 β€– β€– 2 + ξ€· 1 + 𝑐 1 Ξ” 𝑑 𝑐 ξ€Έ β€– Μƒ π‘ž 𝑛 Μƒ 𝑐 𝑛 β€– 2  , ( 3 . 1 9 ) which completes the proof of the stability by virtue of Gronwall’s inequality.

Theorem 3.3 3.3 (stability). Let { 𝐢 𝑛 } 𝑁 𝑛 = 0 be the solution of (2.35a) and (2.35b) subject to the initial value 𝐢 0 . Then there exists a positive constant 𝑐 1 = 𝑐 1 ( β€– 𝐸 π‘ˆ β€– 𝑙 ∞ ( π‘Š 1 , ∞ ) ) independent of β„Ž 𝑐 and Ξ” 𝑑 𝑐 such that β€– 𝐢 β€– 𝑙 ∞ ( 𝐿 2 ) + √ Ξ” 𝑑 𝑐 β€– βˆ‡ 𝐢 β€– 𝑙 ∞ ( 𝐿 2 ) + √ Ξ” 𝑑 𝑐 β€– 𝐢 β€– 𝑙 ∞ ( 𝐿 2 ) + β€– 𝐢 β€– ξ…ž 𝑙 2 ξ€· 𝐻 1 ξ€Έ ≀ 𝑐 1  β€– β€– 𝐢 0 β€– β€– + √ Ξ” 𝑑 𝑐 β€– β€– βˆ‡ 𝐢 0 β€– β€– 2 + β€– Μƒ π‘ž Μƒ 𝑐 β€– 𝑙 2 ( 𝐿 2 )  . ( 3 . 2 0 )

4. Error Estimate Theorem

Now, we turn to derive an optimal priori error estimate in 𝐿 2 -norm for the concentration of approximation (2.35a) and (2.35b). In order to state error estimates, we need the following Lagrange interpolation operator [15] Ξ  β„Ž ∢ 𝐢 0 ( Ξ© ) β†’ 𝑀 β„Ž 𝑐 .

Lemma 4.1. There exists a positive integer π‘˜ such that β€– β€– Ξ  β„Ž β€– β€– 𝑐 βˆ’ 𝑐 𝑠 ≀ 𝐾 β„Ž π‘˜ + 1 βˆ’ 𝑠 β€– 𝑐 β€– π‘˜ + 1 , 𝑠 = 0 , 1 , βˆ€ 𝑐 ∈ 𝐻 π‘˜ + 1 ( Ξ© ) ∩ 𝐢 0 ξ‚€ Ξ©  . ( 4 . 1 ) Let 𝑒 𝑛 = 𝐢 𝑛 βˆ’ Ξ  β„Ž 𝑐 𝑛 and πœ‚ 𝑛 = 𝑐 𝑛 βˆ’ Ξ  β„Ž 𝑐 𝑛 . By Lemma 4.1, it holds that β€– πœ‚ 𝑛 β€– 1 ≀ 𝐾 β„Ž π‘˜ β€– 𝑐 𝑛 β€– π‘˜ + 1 , β€– β€– 𝐷 Ξ” 𝑑 𝑐 πœ‚ 𝑛 β€– β€– ≀ 𝐾 β„Ž π‘˜ √ Ξ” 𝑑 𝑐 β€– β€– β€– πœ• 𝑐 β€– β€– β€– πœ• 𝑑 𝐿 2 ( ( 𝑑 𝑛 , 𝑑 𝑛 + 1 ) ; 𝐻 π‘˜ ) . ( 4 . 2 )

Let 𝑐 ∈ 𝐢 1 ( 𝐿 2 ) ∩ 𝐢 0 ( 𝐻 2 ) , 𝑒 ∈ 𝐢 0 ( 𝐿 ∞ ) , and Μƒ π‘ž ∈ 𝐢 0 ( 𝐿 2 ) be given. Corresponding to π’œ β„Ž 𝑛 + 1 / 2 and β„± β„Ž 𝑛 + 1 / 2 , we introduce linear forms π’œ 𝑛 + 1 / 2 and β„± 𝑛 + 1 / 2 on ( 𝐻 1 ( Ξ© ) ) ξ…ž for 𝑛 = 0 , … , 𝑁 βˆ’ 1 by (2.25) as follows:  π’œ 𝑛 + 1 / 2  ≑ ξ‚΅ πœ™ ξ‚€ 𝑐 , πœ‘ 𝑑 𝑐  𝑑 𝑑 𝑛 + ( 1 / 2 ) ∘ 𝑋 𝑒 𝑛 + ( 1 / 2 ) ξ‚Ά + ξ‚΅ ξ‚€ 𝐹 , πœ‘ 𝑒 𝑛 + ( 1 / 2 )  βˆ’ 1 ξ€· 𝐷 βˆ‡ 𝑐 𝑛 + ( 1 / 2 ) ξ€Έ ∘ 𝑋 𝑒 𝑛 + ( 1 / 2 ) ξ‚Ά + ξ‚΅ ξ‚€ 𝐹 , βˆ‡ πœ‘ d i v 𝑒 𝑛 + ( 1 / 2 )  βˆ’ 𝑇 ξ€· 𝐷 βˆ‡ 𝑐 𝑛 + ( 1 / 2 ) ξ€Έ ∘ 𝑋 𝑒 𝑛 + ( 1 / 2 ) ξ‚Ά + ξ‚€ ξ€· , πœ‘ Μƒ π‘ž 𝑛 + ( 1 / 2 ) 𝑐 𝑛 + ( 1 / 2 ) ξ€Έ ∘ 𝑋 𝑒 𝑛 + ( 1 / 2 )  ,  β„± , πœ‘ 𝑛 + 1 / 2  ≑ ξ‚€ ξ€· , πœ‘ Μƒ π‘ž 𝑛 + ( 1 / 2 ) Μƒ 𝑐 𝑛 + ( 1 / 2 ) ξ€Έ ∘ 𝑋 𝑒 𝑛 + ( 1 / 2 )  , πœ‘ , βˆ€ πœ‘ ∈ 𝐻 1 ( Ξ© ) . ( 4 . 3 ) If 𝑐 is the solution of (1.1a)–(1.1d), we have for 𝑛 = 0 , … , 𝑁 βˆ’ 1 ξ€· π’œ 𝑛 + ( 1 / 2 ) ξ€Έ = ξ€· β„± 𝑐 , πœ‘ 𝑛 + ( 1 / 2 ) ξ€Έ , πœ‘ , βˆ€ πœ‘ ∈ 𝐻 1 ( Ξ© ) . ( 4 . 4 )

We decompose 𝑒 as π’œ β„Ž 𝑛 + ( 1 / 2 ) ξ‚€ β„± 𝑒 = β„Ž 𝑛 + ( 1 / 2 ) βˆ’ β„± 𝑛 + ( 1 / 2 )  + ξ‚€ π’œ 𝑛 + ( 1 / 2 ) βˆ’ π’œ β„Ž 𝑛 + ( 1 / 2 )  𝑐 + π’œ β„Ž 𝑛 + ( 1 / 2 ) πœ‚ . ( 4 . 5 ) In order to estimate the terms of the right-hand side of (4.5), we need the following lemmas.

Lemma 4.2 (see [13]). Assume the above Hypotheses hold, and that the coefficients of the problem (1.1a)–(1.1d) satisfy 𝑣 ∈ 𝐢 0 ( π‘Š 3 , ∞ ( Ξ© ) ) ∩ 𝐢 1 ( π‘Š 2 , ∞ ( Ξ© ) ) , π‘ž ∈ π‘Š 2 , ∞ ( Ξ© ) , and β€– Μ‚ 𝑣 β€– 𝐢 0 ( π‘Š 1 , ∞ ( Ξ© ) ) Ξ” 𝑑 𝑐 < 1 / 2 . Let the solution of (4.4) satisfy 𝑐 ∈ 𝑍 3 , βˆ‡ 𝑐 ∈ 𝑍 3 . Then, for each 𝑛 = 0 , 1 , … , 𝑁 βˆ’ 1 , there exists function πœ‰ 𝐴 𝑛 + ( 1 / 2 ) ∢ Ξ© β†’ 𝑅 , such that π’œ  ξ‚€ 𝑛 + ( 1 / 2 ) βˆ’ π’œ β„Ž 𝑛 + ( 1 / 2 )  ξ‚­ = ξ‚€ πœ‰ 𝑐 , πœ‘ 𝐴 𝑛 + ( 1 / 2 )  , πœ‘ , πœ‘ ∈ 𝐻 1 ( Ξ© ) . ( 4 . 6 ) Moreover, πœ‰ 𝐴 𝑛 + ( 1 / 2 ) ∈ 𝐿 2 ( Ξ© ) and the following estimate holds: β€– β€– πœ‰ 𝐴 𝑛 + ( 1 / 2 ) β€– β€– 0 ≀ Μƒ 𝑐 1 ξ€· Ξ” 𝑑 𝑐 ξ€Έ 2 ξ€· β€– 𝑐 β€– 𝑍 3 + β€– βˆ‡ 𝑐 β€– 𝑍 3 ξ€Έ , ( 4 . 7 ) where Μƒ 𝑐 denotes a constant independent of Ξ” 𝑑 𝑐 .

Lemma 4.3 (see [13]). Assume the above Hypotheses hold, and that the coefficients of the problem (1.1a)–(1.1d) satisfy 𝑣 ∈ 𝐢 0 ( π‘Š 2 , ∞ ( Ξ© ) ) ∩ 𝐢 1 ( π‘Š 1 , ∞ ( Ξ© ) ) , π‘ž ∈ π‘Š 2 , ∞ ( Ξ© ) , and β€– Μ‚ 𝑣 β€– 𝐢 0 ( π‘Š 1 , ∞ ( Ξ© ) ) Ξ” 𝑑 𝑐 < 1 / 2 . Let the solution of (4.4) satisfy 𝑐 ∈ 𝑍 3 , βˆ‡ 𝑐 ∈ 𝑍 3 . Then, for each 𝑛 = 0 , 1 , … , 𝑁 βˆ’ 1 , there exists function πœ‰ 𝑓 𝑛 + ( 1 / 2 ) ∢ Ξ© β†’ 𝑅 , such that 𝐹  ξ‚€ 𝑛 + ( 1 / 2 ) βˆ’ β„± β„Ž 𝑛 + ( 1 / 2 ) ξ‚­ = ξ‚€ πœ‰ , πœ‘ 𝑓 𝑛 + ( 1 / 2 )  , πœ‘ , πœ‘ ∈ 𝐻 1 ( Ξ© ) . ( 4 . 8 ) Moreover, πœ‰ 𝑓 𝑛 + ( 1 / 2 ) ∈ 𝐿 2 ( Ξ© ) and the following estimate holds: β€– β€– πœ‰ 𝑓 𝑛 + ( 1 / 2 ) β€– β€– 0 ≀ Μƒ 𝑐 1 ξ€· Ξ” 𝑑 𝑐 ξ€Έ 2 β€– Μƒ π‘ž Μƒ 𝑐 β€– 𝑍 2 , ( 4 . 9 ) where Μƒ 𝑐 1 denotes a constant independent of Ξ” 𝑑 𝑐 .

Lemma 4.4. Let πœ‚ = 𝑐 βˆ’ Ξ  β„Ž 𝑐 . There exists  π’œ β„Ž 𝑛 + 1 / 2 πœ‚ , 𝑒 𝑛 + 1  ≀ 1 8 β€– β€– βˆ‡ 𝑒 𝑛 + 1 + βˆ‡ 𝑒 𝑛 ∘ 𝑋 𝑛 𝐸 β€– β€– 2 + 1 2 ξ€½ ξ€· 𝐷 βˆ‡ πœ‚ 𝑛 + 1 , βˆ‡ 𝑒 𝑛 + 1 ξ€Έ βˆ’ ( 𝐷 βˆ‡ πœ‚ 𝑛 , βˆ‡ 𝑒 𝑛 ) ξ€Ύ + 𝑐 1 ξ‚» β„Ž 2 π‘˜ ξ‚΅ 1 Ξ” 𝑑 𝑐 β€– β€– β€– πœ• 𝑐 β€– β€– β€– πœ• 𝑑 2 𝐿 2 ξ€· 𝑑 𝑛 , 𝑑 𝑛 + 1 ; 𝐻 π‘˜ ξ€Έ + β€– 𝑐 β€– 2 𝐢 0 ξ€· 𝐻 π‘˜ + 1 ξ€Έ ξ‚Ά + β„Ž 2 π‘˜ + 2 β€– 𝑐 β€– 2 𝐢 0 ξ€· 𝐻 π‘˜ + 1 ξ€Έ + Ξ” 𝑑 𝑐 ξ‚€ β€– βˆ‡ 𝑒 𝑛 β€– 2 + β€– β€– βˆ‡ 𝑒 𝑛 + 1 β€– β€– 2  ξ‚Ό . ( 4 . 1 0 )

Proof. From the definition of π’œ β„Ž 𝑛 + 1 / 2 , we have  π’œ β„Ž 𝑛 + 1 / 2 πœ‚ , 𝑒 𝑛 + 1  ≑  πœ™ πœ‚ 𝑛 + 1 βˆ’ πœ‚ 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 Ξ” 𝑑 𝑐 , 𝑒 𝑛 + 1 ξƒͺ +  𝐷 βˆ‡ πœ‚ 𝑛 + 1 + ( βˆ‡ πœ‚ 𝑛 ) ∘ 𝑋 𝑛 𝐸 2 , βˆ‡ 𝑒 𝑛 + 1 ξƒͺ + Ξ” 𝑑 𝑐 2 ξ€· 𝐷 ( 𝐿 𝑛 βˆ‡ πœ‚ 𝑛 ) ∘ 𝑋 𝑛 𝐸 , βˆ‡ 𝑒 𝑛 + 1 ξ€Έ + Ξ” 𝑑 𝑐 2 ξ€· ( βˆ‡ d i v 𝑣 𝑛 β‹… 𝐷 βˆ‡ πœ‚ 𝑛 ) ∘ 𝑋 𝑛 𝐸 , 𝑒 𝑛 + 1 ξ€Έ + 1 2 ξ€· Μƒ π‘ž 𝑛 + 1 πœ‚ 𝑛 + 1 + ( Μƒ π‘ž 𝑛 πœ‚ 𝑛 ) ∘ 𝑋 𝑛 𝐸 , 𝑒 𝑛 + 1 ξ€Έ = 𝐼 1 + 𝐼 2 + 𝐼 3 + 𝐼 4 + 𝐼 5 . ( 4 . 1 1 ) Since it holds that β€– β€– β€– πœ™ πœ‚ 𝑛 + 1 βˆ’ πœ‚ 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 Ξ” 𝑑 𝑐 β€– β€– β€– ≀ 𝑐 1 √ Ξ” 𝑑 𝑐 ξ‚» β€– β€– β€– πœ• πœ‚ β€– β€– β€– πœ• 𝑑 𝐿 2 ( 𝑑 𝑛 , 𝑑 𝑛 + 1 ; 𝐿 2 ) + β€– πœ‚ β€– 𝐿 2 ( 𝑑 𝑛 , 𝑑 𝑛 + 1 ; 𝐻 1 ) ξ‚Ό , ( 4 . 1 2 ) we have 𝐼 1 ≀ ξ€· 1 + 𝑐 1 Ξ” 𝑑 𝑐 ξ€Έ ξ€· 1 + 2 𝑐 1 ξ€Έ 2 𝐾 2 β„Ž 2 π‘˜ 2 Ξ” 𝑑 𝑐 ξ‚΅ β€– β€– β€– πœ• 𝑐 β€– β€– β€– πœ• 𝑑 2 𝐿 2 ξ€· 𝑑 𝑛 , 𝑑 𝑛 + 1 ; 𝐻 π‘˜ ξ€Έ + β€– 𝑐 β€– 2 𝐿 2 ξ€· 𝑑 𝑛 , 𝑑 𝑛 + 1 ; 𝐻 π‘˜ + 1 ξ€Έ ξ‚Ά + 1 2 β€– β€– 𝑒 𝑛 + 1 β€– β€– 2 . ( 4 . 1 3 ) To estimate 𝐼 2 , we divide it into three parts as 𝐼 2 = 1 2 ξ€Ί ξ€· 𝐷 βˆ‡ πœ‚ 𝑛 + 1 , βˆ‡ 𝑒 𝑛 + 1 ξ€Έ βˆ’ ( 𝐷 βˆ‡ πœ‚ 𝑛 , βˆ‡ 𝑒 𝑛 ) ξ€» + 𝐼 2 1 + 𝐼 2 2 , ( 4 . 1 4 ) where 𝐼 2 1 = 1 2 ξ€Ί ξ€· ( 𝐷 βˆ‡ πœ‚ 𝑛 ) ∘ 𝑋 𝑛 𝐸 , βˆ‡ 𝑒 𝑛 + 1 + βˆ‡ 𝑒 𝑛 ∘ 𝑋 𝑛 𝐸 , 𝐼 ξ€Έ ξ€» 2 2 = 1 2 ξ€Ί ( 𝐷 βˆ‡ πœ‚ 𝑛 , βˆ‡ 𝑒 𝑛 ξ€· ) βˆ’ ( 𝐷 βˆ‡ πœ‚ 𝑛 ) ∘ 𝑋 𝑛 𝐸 , βˆ‡ 𝑒 𝑛 ∘ 𝑋 𝑛 𝐸 . ξ€Έ ξ€» ( 4 . 1 5 ) The term 𝐼 2 1 is estimated as 𝐼 2 1 ≀ 1 8 β€– β€– βˆ‡ 𝑒 𝑛 + 1 + βˆ‡ 𝑒 𝑛 ∘ 𝑋 𝑛 𝐸 β€– β€– 2 β€– β€– + 2 ( 𝐷 βˆ‡ πœ‚ 𝑛 ) ∘ 𝑋 𝑛 𝐸 β€– β€– 2 ≀ 1 8 β€– β€– βˆ‡ 𝑒 𝑛 + 1 + βˆ‡ 𝑒 𝑛 ∘ 𝑋 𝑛 𝐸 β€– β€– 2 + ξ€· 1 + 𝑐 1 Ξ” 𝑑 𝑐 ξ€Έ 𝐾 ξ€· 𝑑 βˆ— ξ€Έ 2 β„Ž 2 π‘˜ β€– 𝑐 𝑛 β€– 2 𝐻 π‘˜ + 1 . ( 4 . 1 6 ) Using the transformation 𝑦 = 𝑋 𝑛 𝐸 ( π‘₯ ) , we have 𝐼 2 2 = 1 2 ξ€œ Ξ© 𝐷 βˆ‡ πœ‚ 𝑛 β‹… βˆ‡ 𝑒 𝑛 ξƒ― ξ‚΅ 1 βˆ’ d e t πœ• 𝑦 ξ‚Ά πœ• π‘₯ βˆ’ 1 ξƒ° ≀ 𝑐 𝑑 𝑦 𝑑 π‘₯ 1 ξ€· 𝑑 βˆ— ξ€Έ 2 Ξ” 𝑑 𝑐 4 ξ€· 𝐾 2 β„Ž 2 π‘˜ β€– 𝑐 𝑛 β€– 2 𝐻 π‘˜ + 1 + β€– βˆ‡ 𝑒 𝑛 β€– 2 ξ€Έ . ( 4 . 1 7 ) It is clear that 𝐼 3 ≀ ξ€· 1 + 𝑐 1 Ξ” 𝑑 𝑐 ξ€Έ 𝑐 2 1 ξ€· 𝑑 βˆ— ξ€Έ 2 Ξ” 𝑑 𝑐 𝐾 2 β„Ž 2 π‘˜ 4 β€– 𝑐 𝑛 β€– 2 𝐻 π‘˜ + 1 + 1 4 β€– β€– βˆ‡ 𝑒 𝑛 + 1 β€– β€– 2 , 𝐼 ( 4 . 1 8 ) 4 ≀ ξ€· 1 + 𝑐 1 Ξ” 𝑑 𝑐 ξ€Έ 𝑐 2 1 ξ€· 𝑑 βˆ— ξ€Έ 2 Ξ” 𝑑 𝑐 𝐾 2 β„Ž 2 π‘˜ 4 β€– 𝑐 𝑛 β€– 2 𝐻 π‘˜ + 1 + 1 4 β€– β€– βˆ‡ 𝑒 𝑛 + 1 β€– β€– 2 . ( 4 . 1 9 )
Term 𝐼 5 can be divided it into three parts like term 𝐼 2 as 𝐼 5 = 𝐼 5 1 + 𝐼 5 2 + 𝐼 5 3 , ( 4 . 2 0 ) where 𝐼 5 1 = 1 2  ξ‚Έ ξ‚΅ Μƒ π‘ž 𝑛 + 1 πœ‚ 𝑛 + 1 ,  Μƒ π‘ž 𝑛 + 1 𝑒 𝑛 + 1 ξ‚Ά βˆ’ ξ‚€ √ Μƒ π‘ž 𝑛 πœ‚ 𝑛 , √ Μƒ π‘ž 𝑛 𝑒 𝑛  ξ‚Ή , 𝐼 5 2 = 1 2 √  ξ‚€ ξ‚€ Μƒ π‘ž 𝑛 πœ‚ 𝑛  ∘ 𝑋 𝑛 𝐸 , √ Μƒ π‘ž 𝑛 ∘ 𝑋 𝑛 𝐸 ξ€· 𝑒 𝑛 + 1 + 𝑒 𝑛 ∘ 𝑋 𝑛 𝐸 ξ€Έ , 𝐼  ξ‚„ 5 3 = 1 2 √  ξ‚€ Μƒ π‘ž 𝑛 πœ‚ 𝑛 , √ Μƒ π‘ž 𝑛 𝑒 𝑛  βˆ’ √ ξ‚€ ξ‚€ Μƒ π‘ž 𝑛 πœ‚ 𝑛  ∘ 𝑋 𝑛 𝐸 , ξ‚€ √ Μƒ π‘ž 𝑛 𝑒 𝑛  ∘ 𝑋 𝑛 𝐸 .  ξ‚„ ( 4 . 2 1 ) Thus, the same kind of computations used for 𝐼 2 leads to the following estimate: 𝐼 5 ≀ 𝐼 5 1 + ξ‚΅ 𝑐 1 𝑐 3 Ξ” 𝑑 𝑐 4 + ξ€· 1 + 𝑐 1 Ξ” 𝑑 𝑐 ξ€Έ 𝑐 3 ξ‚Ά 𝐾 2 β„Ž 2 π‘˜ β€– 𝑐 𝑛 β€– 2 𝐻 π‘˜ + 1 + ξ€· 2 𝑐 1 + 𝑐 2 1 Ξ” 𝑑 𝑐 ξ€Έ 𝑐 3 Ξ” 𝑑 𝑐 8  β€– β€– √ Μƒ π‘ž 𝑛 πœ‚ 𝑛 β€– β€– 2 + β€– β€– β€–  Μƒ π‘ž 𝑛 + 1 πœ‚ 𝑛 + 1 β€– β€– β€– 2 ξƒͺ + 1 8 β€– β€– β€–  Μƒ π‘ž 𝑛 + 1 πœ‚ 𝑛 + 1 + ξ‚€ √ Μƒ π‘ž 𝑛 πœ‚ 𝑛  ∘ 𝑋 𝑛 𝐸 β€– β€– β€– 2 . ( 4 . 2 2 ) Gathering (4.11)–(4.19) together, we complete the proof.

We now turn to estimate 𝑒 𝑛 . From the definition of π’œ β„Ž 𝑛 + ( 1 / 2 ) and π’œ 𝑛 + ( 1 / 2 ) , we see  π’œ β„Ž 𝑛 + ( 1 / 2 ) 𝑒 , 𝑒 𝑛 + 1 ξ‚­ =  π’œ β„Ž 𝑛 + ( 1 / 2 ) πœ‚ , 𝑒 𝑛 + 1 ξ‚­ + π’œ  ξ‚€ 𝑛 + ( 1 / 2 ) βˆ’ π’œ β„Ž 𝑛 + ( 1 / 2 )  𝑐 , 𝑒 𝑛 + 1 ξ‚­ + β„±  ξ‚€ β„Ž 𝑛 + ( 1 / 2 ) βˆ’ β„± 𝑛 + ( 1 / 2 )  , 𝑒 𝑛 + 1 ξ‚­ . ( 4 . 2 3 ) From Lemmas 4.2, 4.3 and 4.4, we obtain 𝐷 Ξ” 𝑑 𝑐 ξ‚΅ 1 2 β€– β€– 𝑒 𝑛 + 1 β€– β€– 2 πœ™ + Ξ” 𝑑 𝑐 4 β€– βˆ‡ 𝑒 𝑛 β€– 2 𝐷 + Ξ” 𝑑 𝑐 4 β€– 𝑒 𝑛 β€– 2 π‘ž ξ‚Ά + 1 2 Ξ” 𝑑 𝑐 β€– β€– 𝑒 𝑛 + 1 βˆ’ 𝑒 𝑛 ∘ 𝑋 𝑛 𝑅 𝐾 β€– β€– 2 + 1 8 β€– β€– βˆ‡ 𝑒 𝑛 + 1 + βˆ‡ 𝑒 𝑛 ∘ 𝑋 𝑛 𝐸 β€– β€– 2 𝐷 + 1 4 β€– β€– Μƒ π‘ž 𝑛 + 1 𝑒 𝑛 + 1 + Μƒ π‘ž 𝑛 𝑒 𝑛 ∘ 𝑋 𝑛 𝐸 β€– β€– 2 ≀ 𝑐 1  β€– β€– 𝑒 𝑛 + 1 β€– β€– 2 + Ξ” 𝑑 𝑐 β€– β€– βˆ‡ 𝑒 𝑛 + 1 β€– β€– 2 𝐷 + β€– 𝑒 𝑛 β€– 2 + Ξ” 𝑑 𝑐 β€– βˆ‡ 𝑒 𝑛 β€– 2 𝐷 + Ξ” 𝑑 𝑐 β€– 𝑒 𝑛 β€– 2 π‘ž  + 1 2 ξ€½ ξ€· 𝐷 βˆ‡ πœ‚ 𝑛 + 1 , βˆ‡ 𝑒 𝑛 + 1 ξ€Έ βˆ’ ( 𝐷 βˆ‡ πœ‚ 𝑛 , βˆ‡ 𝑒 𝑛 ) ξ€Ύ + 𝑐 2 ξ‚» β„Ž 2 π‘˜ ξ‚΅ 1 Ξ” 𝑑 𝑐 β€– β€– β€– πœ• 𝑐 β€– β€– β€– πœ• 𝑑 2 𝐿 2 ξ€· 𝑑 𝑛 , 𝑑 𝑛 + 1 ; 𝐻 π‘˜ ξ€Έ + β€– 𝑐 β€– 2 𝐢 0 ξ€· 𝐻 π‘˜ + 1 ξ€Έ ξ‚Ά + β„Ž 2 π‘˜ + 2 β€– 𝑐 β€– 2 𝐢 0 ξ€· 𝐻 π‘˜ + 1 ξ€Έ + ξ€· Ξ” 𝑑 𝑐 ξ€Έ 4 ξ‚€ β€– β€– | | | | β€– β€– βˆ‡ 𝑐 2 3 + β€– | 𝑐 | β€– 2 3  ξ‚Ό . ( 4 . 2 4 ) Summing up the above equation about time from 𝑑 = 0 to 𝑑 𝑛 , we get π‘Ž 𝑛 + 1 + Ξ” 𝑑 𝑐 𝑛  𝑗 = 0 𝑏 𝑗 ≀ 𝑐 1 Ξ” 𝑑 𝑐 𝑛  𝑗 = 0 π‘Ž 𝑗 + Ξ” 𝑑 𝑐 2 ξ€½ ξ€· 𝐷 βˆ‡ πœ‚ 𝑛 + 1 , βˆ‡ 𝑒 𝑛 + 1 ξ€Έ βˆ’ ξ€· 𝐷 βˆ‡ πœ‚ 0 , βˆ‡ 𝑒 0 ξ€Έ ξ€Ύ + 𝑐 2 ξ‚» β„Ž 2 π‘˜ ξ‚΅ 1 Ξ” 𝑑 𝑐 β€– β€– β€– πœ• 𝑐 β€– β€– β€– πœ• 𝑑 2 𝐿 2 ξ€· 𝑑 𝑛 , 𝑑 𝑛 + 1 ; 𝐻 π‘˜ ξ€Έ + β€– 𝑐 β€– 2 𝐢 0 ξ€· 𝐻 π‘˜ + 1 ξ€Έ ξ‚Ά + β„Ž 2 π‘˜ + 2 β€– 𝑐 β€– 2 𝐢 0 ξ€· 𝐻 π‘˜ + 1 ξ€Έ + ξ€· Ξ” 𝑑 𝑐 ξ€Έ 4 ξ‚€ β€– β€– | | | | β€– β€– βˆ‡ 𝑐 2 3 + β€– | 𝑐 | β€– 2 3  ξ‚Ό , ( 4 . 2 5 ) where π‘Ž 𝑗 = 1 2 β€– β€– 𝑒 𝑗 β€– β€– 2 πœ™ + Ξ” 𝑑 𝑐 4 β€– β€– βˆ‡ 𝑒 𝑗 β€– β€– 2 𝐷 + Ξ” 𝑑 𝑐 4 β€– β€– 𝑒 𝑗 β€– β€– 2 π‘ž , 𝑏 𝑗 = 1 2 Ξ” 𝑑 𝑐 β€– β€– 𝑒 𝑗 + 1 βˆ’ 𝑒 𝑗 ∘ 𝑋 𝑗 𝑅 𝐾 β€– β€– 2 + 1 8 β€– β€– βˆ‡ 𝑒 𝑗 + 1 + βˆ‡ 𝑒 𝑗 ∘ 𝑋 𝑗 𝐸 β€– β€– 2 𝐷 + 1 4 β€– β€– Μƒ π‘ž 𝑗 + 1 𝑒 𝑗 + 1 + Μƒ π‘ž 𝑗 𝑒 𝑗 ∘ 𝑋 𝑗 𝐸 β€– β€– 2 . ( 4 . 2 6 ) Using the estimates Ξ” 𝑑 𝑐 2 ξ€· 𝐷 βˆ‡ πœ‚ 𝑗 , βˆ‡ 𝑒 𝑗 ξ€Έ ≀ Ξ” 𝑑 𝑐 8 β€– β€– βˆ‡ 𝑒 𝑗 β€– β€– 2 + 𝑐 1 Ξ” 𝑑 𝑐 β„Ž 2 π‘˜ β€– β€– 𝑐 𝑗 β€– β€– 2 π‘˜ + 1 ( 4 . 2 7 ) for 𝑗 = 0 and 𝑛 + 1 , and Gronwall’s inequality, we derive the following error estimate.

Theorem 4.5 (error estimate). Let { 𝐢 𝑛 } 𝑁 𝑛 = 0 be the solution of (2.35a) and (2.35b) subject to the initial value 𝐢 0 . Then there exists a positive constant 𝑐 1 = 𝑐 1 ( β€– 𝐸 π‘ˆ β€– 𝐢 0 ( π‘Š 1 , ∞ ) ) independent of β„Ž 𝑐 and Ξ” 𝑑 𝑐 such that β€– 𝑐 βˆ’ 𝐢 β€– 𝑙 ∞ ( 𝐿 2 ) + √ Ξ” 𝑑 𝑐 β€– β€– βˆ‡ ( 𝑐 βˆ’ 𝐢 ) 𝑙 ∞ ( 𝐿 2 ) + √ Ξ” 𝑑 𝑐 β€– 𝑐 βˆ’ 𝐢 β€– 𝑙 ∞ ( 𝐿 2 ) + β€– 𝑐 βˆ’ 𝐢 β€– ξ…ž 𝑙 2 ξ€· 𝐻 1 ξ€Έ ≀ 𝑐 2 ξ‚» β„Ž π‘˜ ξ‚΅ β€– β€– β€– πœ• 𝑐 β€– β€– β€– πœ• 𝑑 𝐿 2 ( β„Ž π‘˜