# A posteriori error estimates of characteristic mixed finite elements for convection-diffusion control problems

• Yuelong Tang and Yuchun Hua
From the journal Open Mathematics

## Abstract

In this article, we consider fully discrete characteristic mixed finite elements for convection-diffusion optimal control problems. We use the characteristic line method to treat the hyperbolic part of the state equation as a directional derivative. The state and the co-state are discretized by the lowest order Raviart-Thomas mixed finite element spaces and the control is approximated by piecewise constant functions. Using some proper duality problems, we derive a posteriori error estimates for the scalar functions. Such estimates are not available in the literature. A numerical example is presented to validate the theoretical results.

MSC 2010: 49J20; 65N30

## 1 Introduction

In the recent 30 years, optimal control problems (OCPs) have been extensively utilized in many aspects of modern life such as social, economic, scientific, and engineering numerical simulation [1]. Thus, they must be solved successfully with efficient numerical methods. Among these numerical methods, the finite element method (FEM) is a good choice for solving OCPs governed by partial differential equations (PDEs). There have been extensive studies on convergence or superconvergence of FEMs for OCPs, see [2,3,4, 5,6,7, 8,9]. A systematic introduction of FEM for PDE OCPs can be found in [10,11, 12,13].

Adaptive FEM has been investigated extensively. It has become one of the most popular methods in cientific computation and numerical modeling. It ensures a higher density of nodes in a certain area of the given domain, where the solution is more difficult to approximate, indicated by a posteriori error estimators. Hence, it is an important approach to boost the accuracy and efficiency of finite element discretization. There is a lot of work concentrating on the adaptivity of various OCPs. For example, [14,15, 16,17,18, 19,20].

In flow control problem or temperature control problem, their objective functional contains not only the state variable but also its gradient. At this moment, the accuracy of the gradient is important in the numerical discretization of the coupled state equations. The mixed finite element method (MFEM) is appropriate for the state equations in such cases since both the scalar variable and its flux variable can be approximated to the same accuracy, for example, [21,22,23, 24,25].

Recently, more and more attention has been paid to elliptic convection-diffusion OCPs [26,27,28]. Fu and Rui have studied a priori and a posteriori error estimates of characteristic finite element methods (CFEMs) for time-dependent convection-diffusion OCPs in [29] and [30], respectively. They also considered a priori error estimates of characteristic mixed finite element method (CMFEM) for time-dependent convection-diffusion OCPs in [31]. In [32], Sun and Ma investigated a nonoverlapping domain decomposition method combined with the characteristic method for OCPs governed by parabolic convection-diffusion equations.

In this work, we consider a posteriori error estimates of fully discrete CMFEM for OCPs governed by convection-diffusion equations. The characteristic approximation is applied to handle the convection term, and the lowest order RT mixed finite element spatial approximation is adopted to deal with the diffusion term. The OCPs that we are interested in are as follows:

(1) min u K 1 2 0 T ( p p d 2 + y y d 2 + u 2 ) d t ,

(2) y t + b y + div p + c y = f + u , x Ω , t J ,

(3) p = A y , x Ω , t J ,

(4) y = 0 , x Ω , t J ,

(5) y ( x , 0 ) = y 0 ( x ) , x Ω ,

where Ω R 2 is a bounded convex polygon with Lipschitz continuous boundary Ω , J = [ 0 , T ] . Let p d ( L 2 ( J ; H 1 ( Ω ) ) ) 2 , y d L 2 ( J ; H 1 ( Ω ) ) , f L 2 ( J ; L 2 ( Ω ) ) , y 0 ( x ) H 0 1 ( Ω ) , c W 1 , ( Ω ) , A = ( a i j ( x ) ) 2 × 2 C ( Ω ¯ ; R 2 × 2 ) such that X T A X c 1 X R 2 2 , X R 2 , where c 1 > 0 . Usually, b = ( b 1 ( x , t ) , b 2 ( x , t ) ) T denotes a velocity field in the flow control, we assume that b ( C ( J ; C 0 1 ( Ω ¯ ) ) ) 2 and the flow is incompressible, i.e.,

b = 0 , x Ω , t J .

Moreover, we assume that the constraint on the control is an obstacle such that

K = { u L 2 ( J ; L 2 ( Ω ) ) : u ( x , t ) 0 , a.e. in Ω × J } .

In this article, we adopt the standard notation W m , p ( Ω ) for Sobolev spaces on Ω with a norm m , p given by v m , p p = α m D α v L p ( Ω ) p , a semi-norm m , p given by v m , p p = α = m D α v L p ( Ω ) p . We set W 0 m , p ( Ω ) = { v W m , p ( Ω ) : v Ω = 0 } . For p = 2 , we denote H m ( Ω ) = W m , 2 ( Ω ) , H 0 m ( Ω ) = W 0 m , 2 ( Ω ) , and m = m , 2 , = 0 , 2 .

We denote by L s ( 0 , T ; W m , p ( Ω ) ) the Banach space of all L s integrable functions from J into W m , p ( Ω ) with norm v L s ( J ; W m , p ( Ω ) ) = 0 T v W m , p ( Ω ) s d t 1 s for s [ 1 , ) , and the standard modification for s = . Similarly, one can define the spaces H 1 ( J ; W m , p ( Ω ) ) and C k ( J ; W m , p ( Ω ) ) . In addition, C denotes a general positive constant independent of h and k , where h is the spatial mesh size for the control and state discretization and k is the time increment.

The plan of this article is as follows. In Section 2, we shall construct a CMFEM approximation for the OCPs (1)–(5). In Section 3, we derive a posteriori L 2 ( 0 , T ; L 2 ( Ω ) ) error estimates by using two duality problems. A numerical example is provided in Section 4.

## 2 CMFEM approximation of convection-diffusion OCPs

In this section, we shall construct a CMFEM and backward Euler discretization approximation of convection-diffusion OCPs (1)–(5). To fix the idea, we shall take the state spaces L = L 2 ( J ; V ) and Q = H 1 ( J ; W ) , where V and W are defined as follows:

V = H ( div ; Ω ) = { v ( L 2 ( Ω ) ) 2 , div v L 2 ( Ω ) } , W = L 2 ( Ω ) .

The Hilbert space V is equipped with the following norm:

v H ( div ; Ω ) = ( v 0 , Ω 2 + div v 0 , Ω 2 ) 1 / 2 .

Now, we take account of the hyperbolic part of (2), namely y t + b y , as a directional derivative. Let s denote the unit vector in the direction of ( b 1 , b 2 , 1 ) in Ω × J , and set

λ = ( b 2 + 1 ) 1 / 2 = ( b 1 2 + b 2 2 + 1 ) 1 / 2 .

Then equation (2) is equivalent to the following form:

λ y s + div p + c y = f + u .

We recast (1)–(5) as the following weak form: find ( p , y , u ) L × Q × K such that

(6) min u K 1 2 0 T ( p p d 2 + y y d 2 + u 2 ) d t

(7) ( λ y s , w ) + ( div p , w ) + ( c y , w ) = ( f + u , w ) , w W , t J ,

(8) ( A 1 p , v ) ( y , div v ) = 0 , v V , t J ,

(9) y ( x , 0 ) = y 0 ( x ) , x Ω .

It follows from [1] that the OCPs (6)–(9) have a unique solution ( p , y , u ) L × Q × K and that a triplet ( p , y , u ) is the solution of (6)–(9) if and only if there is a co-state ( q , z ) L × Q such that ( p , y , q , z , u ) satisfies the following optimality conditions:

(10) ( λ y s , w ) + ( div p , w ) + ( c y , w ) = ( f + u , w ) , w W , t J ,

(11) ( A 1 p , v ) ( y , div v ) = 0 , v V , t J ,

(12) y ( x , 0 ) = y 0 ( x ) , x Ω ,

(13) ( λ z s , w ) + ( div q , w ) + ( c z , w ) = ( y y d , w ) , w W , t J ,

(14) ( A 1 q , v ) ( z , div v ) = ( p p d , v ) , v V , t J ,

(15) z ( x , T ) = 0 , x Ω ,

(16) ( u + z , u ˜ u ) 0 , u ˜ K ,

where ( , ) is the inner product of L 2 ( Ω ) .

We now consider the time discretization. Let N Z + , k = T / N , and t i = i k , i = 0 , 1 , , N . Approximation y s i ( x ) = y s ( x , t i ) by a backward-Euler difference quotient in the s -direction, that is,

y s i ( x ) y i ( x ) y i 1 ( x b ( x , t i ) k ) k 1 + b ( x , t i ) 2 .

Let G ( x , t ; t ) be the approximate characteristic curve passing through point x at time t , which is defined by

G ( x , t ; t ) x b ( x , t ) ( t t ) .

We denote by x ¯ = G ( x , t i ; t i 1 ) the foot at time t i 1 of the characteristic curve with head x at time t i , and f ¯ ( x ) = f ( x ¯ ) , then

λ i y s i y i y ¯ i 1 k .

Next, we use the MFEM for spatial discretization. Let T h be regular triangulations of Ω . h τ is the diameter of element τ and h = max τ T h { h τ } . Let V h × W h V × W denote the lowest order Raviart-Thomas space [33] associated with the triangulations T h of Ω . P k denotes the space of polynomials of total degree at most k ( k 0 ). Let V ( τ ) = { v P 0 2 ( τ ) + x P 0 ( τ ) } , W ( τ ) = P 0 ( τ ) . We define

V h { v h V : v h τ V ( τ ) , τ T h } , W h { w h W : w h τ W ( τ ) , τ T h } , K h { w h W h : w h τ = max { 0 , w h } , τ T h } .

Then a fully discrete CMFEM approximation of (6)–(9) is to find ( p h i , y h i , u h i ) V h × W h × K h , i = 1 , 2 , , N , such that

(17) min u h i K h 1 2 i = 1 N k ( p h i p d i 2 + y h i y d i 2 + u h i 2 )

(18) y h i y ¯ h i 1 k , w h + ( div p h i , w h ) + ( c y h i , w h ) = ( f i + u h i , w h ) , w h W h ,

(19) ( A 1 p h i , v h ) ( y h i , div v h ) = 0 , v h V h ,

(20) y h 0 ( x ) = y 0 h ( x ) , x Ω ,

where y 0 h ( x ) = R h y 0 ( x ) and R h is the L 2 ( Ω ) -projection operator, which will be specific later.

It follows from [31] that the control problem (17)–(20) has a unique solution ( p h i , y h i , u h i ) , i = 1 , 2 , , N , and that a triplet ( p h i , y h i , u h i ) V h × W h × K h is the solution of (17)–(20) if and only if there is a co-state ( q h i 1 , z h i 1 ) V h × W h , i = N , N 1 , , 1 such that ( p h i , y h i , q h i 1 , z h i 1 , u h i ) ( V h × W h ) 2 × K h satisfies the following optimality conditions:

(21) y h i y ¯ h i 1 k , w h + ( div p h i , w h ) + ( c y h i , w h ) = ( f i + u h i , w h ) , w h W h ,

(22) ( A 1 p h i , v h ) ( y h i , div v h ) = 0 , v h V h ,

(23) y h 0 ( x ) = y 0 h ( x ) , x Ω ,

(24) z h i 1 z ¯ ¯ h i J i k , w h + ( div q h i 1 , w h ) + ( c z h i 1 , w h ) = ( y h i y d i , w h ) , w h W h ,

(25) ( A 1 q h i 1 , v h ) ( z h i 1 , div v h ) = ( p h i p d i , v h ) , v h V h ,

(26) z h N ( x ) = 0 , x Ω ,

(27) ( u h i + z h i 1 , u ˜ h u h i ) 0 , u ˜ h K h ,

where z ¯ ¯ h i = z h i ( x ¯ ¯ ) , and x ¯ ¯ represents the head of the characteristic curve with foot x at time t i 1 , namely,

x = G ( x ¯ ¯ , t i ; t i 1 ) .

We denote by J i det D G ( x , t i ; t i 1 ) 1 the determinant of the Jacobian transformation from G to x . Similar to the discussion in [29], we know that

det D G ( x ¯ ¯ , t i ; t i 1 ) = 1 ( b i ) k + O ( k 2 ) ,

since the flow is incompressible, i.e., b = 0 , x Ω , t ( 0 , T ) . Thus, J i = 1 for sufficiently small k .

For i = 0 and i = N , we let

(28) ( A 1 p h 0 , v h ) ( y h 0 , div v h ) = 0 , v h V h ,

(29) ( A 1 q h N , v h ) ( z h N , div v h ) = ( p h N p d N , v h ) , v h V h .

For i = 1 , 2 , , N , let

Y h ( t i 1 , t i ] = ( ( t i t ) y h i 1 + ( t t i 1 ) y h i ) / k , Z h ( t i 1 , t i ] = ( ( t i t ) z h i 1 + ( t t i 1 ) z h i ) / k , P h ( t i 1 , t i ] = ( ( t i t ) p h i 1 + ( t t i 1 ) p h i ) / k , Q h ( t i 1 , t i ] = ( ( t i t ) q h i 1 + ( t t i 1 ) q h i ) / k , U h ( t i 1 , t i ] = u h i .

For any function w C ( J ; L 2 ( Ω ) ) and i = 1 , 2 , , N , let

w ˆ ( x , t ) t ( t i 1 , t i ] = w ( x , t i ) , w ˜ ( x , t ) t ( t i 1 , t i ] = w ( x , t i 1 ) .

Moreover, we let

p ¯ d ( t i 1 , t i ] = ( ( t i t ) p d i + ( t t i 1 ) p d i + 1 ) / k , i = 1 , 2 , , N 1 , p ¯ d ( t N 1 , t N ] = p d N , P ¯ h ( t i 1 , t i ] = ( ( t i t ) p h i + ( t t i 1 ) p h i + 1 ) / k , i = 1 , 2 , , N 1 , P ¯ h ( t N 1 , t N ] = p h N .

Then the optimality conditions (21)–(27) satisfying

(30) ( Y h , t , w h ) + ( div P ˆ h , w h ) + ( c Y ˆ h , w h ) = f ˆ + U h Y ˜ h Y ˜ h ( x ¯ ) k , w h , w h W h ,

(31) ( A 1 P ˆ h , v h ) ( Y ˆ h , div v h ) = 0 , v h V h ,

(32) Y h ( x , 0 ) = y 0 h ( x ) , x Ω ,

(33) ( Z h , t , w h ) + ( div Q ˜ h , w h ) + ( c Z ˜ h , w h ) = Y ˆ h y ˆ d + Z ˆ h ( x ¯ ¯ ) J ˆ Z ˆ h k , w h , w h W h ,

(34) ( A 1 Q ˜ h , v h ) ( Z ˜ h , div v h ) = ( P ˆ h p ˆ d , v h ) , v h V h ,

(35) Z h ( x , T ) = 0 , x Ω ,

(36) ( U h + Z ˜ h , u ˜ h U h ) 0 , u ˜ h K h .

In the rest of the article, we shall use some intermediate variables. For any control function U h K h , we first define the state solution ( p ( U h ) , y ( U h ) , q ( U h ) , z ( U h ) ) satisfies

(37) ( y t ( U h ) + b y ( U h ) , w ) + ( div p ( U h ) , w ) + ( c y ( U h ) , w ) = ( f + U h , w ) , w W ,

(38) ( A 1 p ( U h ) , v ) ( y ( U h ) , div v ) = 0 , v V ,

(39) y ( U h ) ( x , 0 ) = y 0 ( x ) , x Ω ,

(40) ( z t ( U h ) + b z ( U h ) , w ) + ( div q ( U h ) , w ) + ( c z ( U h ) , w ) = ( y ( U h ) y d , w ) , w W ,

(41) ( A 1 q ( U h ) , v ) ( z ( U h ) , div v ) = ( p ( U h ) p d , v ) , v V ,

(42) z ( U h ) ( x , T ) = 0 , x Ω .

Let R h : W W h be the L 2 ( Ω ) -projection into W h [34], which satisfies for any w W

(43) ( R h w w , χ ) = 0 , w W , χ W h ,

(44) R h w w 0 , q C h w 1 , q , w W W 1 , q ( Ω ) .

Let Π h : V V h be the Raviart-Thomas projection operator [35], which satisfies: for any v V

(45) E w h ( v Π h v ) ν E d s = 0 , w h W h , E h ,

(46) τ ( v Π h v ) v h d x d y = 0 , v h V h , τ T h ,

where h denotes the set of element sides in T h .

We have the commuting diagram property

(47) div Π h = R h div : V W h and div ( I Π h ) V W h ,

where and after, I denotes the identity operator.

Furthermore, the interpolation operator Π h satisfies a local error estimate:

(48) v Π h v 0 , Ω C h v 1 , T h , v V H 1 ( T h ) .

## 3 A posteriori error estimates

In this section, we derive a posteriori error estimates of fully discrete CMFEM for parabolic convection-diffusion OCPs. In order to the following analysis, we divide the domain Ω into three parts:

Ω = { x Ω : Z ˜ h ( x ) 0 } , Ω 0 = { x Ω : Z ˜ h ( x ) > 0 , U h ( x ) = 0 } , Ω + = { x Ω : Z ˜ h ( x ) > 0 , U h ( x ) > 0 } .

It is easy to see that the partition of the above three subsets is dependent on t . For all t , the three subsets are not intersected each other, and

Ω ¯ = Ω ¯ Ω ¯ 0 Ω ¯ + .

First, let us derive the a posteriori error estimates for the control u .

## Theorem 3.1

Let ( y , p , z , q , u ) and ( Y h , P h , Z h , Q h , U h ) be the solutions of (10)–(16) and (30)–(36), respectively. Then we have

(49) u U h L 2 ( J ; L 2 ( Ω ) ) 2 C η 1 2 + Z ˜ h z ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 ,

where

η 1 2 = U h + Z ˜ h L 2 ( J ; L 2 ( Ω Ω + ) ) 2 .

## Proof

It follows from (19) that

(50) u U h L 2 ( J ; L 2 ( Ω ) ) 2 = 0 T ( u U h , u U h ) d t = 0 T ( u + z , u U h ) d t + 0 T ( U h + Z ˜ h , U h u ) d t + 0 T ( Z ˜ h z ( U h ) , u U h ) d t + 0 T ( z ( U h ) z , u U h ) d t 0 T ( U h + Z ˜ h , U h u ) d t + 0 T ( Z ˜ h z ( U h ) , u U h ) d t + 0 T ( z ( U h ) z , u U h ) d t I 1 + I 2 + I 3 .

We first estimate I 1 . Note that

(51) I 1 = 0 T ( U h + Z ˜ h , U h u ) d t = 0 T Ω Ω + ( U h + Z ˜ h ) ( U h u ) d x d t + 0 T Ω 0 ( U h + Z ˜ h ) ( U h u ) d x d t .

It follows from Hölder’s inequality and Young’s inequality that

(52) 0 T Ω Ω + ( U h + Z ˜ h ) ( U h u ) d x d t C ( δ ) U h + Z ˜ h L 2 ( J ; L 2 ( Ω Ω + ) ) 2 + δ u U h L 2 ( J ; L 2 ( Ω Ω + ) ) 2 = C ( δ ) η 1 2 + δ u U h L 2 ( J ; L 2 ( Ω ) ) 2 ,

Furthermore, we have that

U h + Z ˜ h Z ˜ h > 0 , U h u = 0 u 0 on Ω 0 .

It yields that

(53) 0 T Ω 0 ( U h + Z ˜ h ) ( U h u ) d x d t 0 .

Then (51)–(53) imply that

(54) I 1 C ( δ ) η 1 2 + δ u U h L 2 ( J ; L 2 ( Ω ) ) 2 .

Moreover, it is clear that

(55) I 2 = 0 T ( Z ˜ h z ( U h ) , u U h ) d t C ( δ ) Z ˜ h z ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 + δ u U h L 2 ( J ; L 2 ( Ω ) ) 2 .

Now we turn to I 3 . Note that

y ( x , 0 ) = y ( U h ) ( x , 0 ) = y 0 ( x ) and z ( x , T ) = z ( U h ) ( x , T ) = 0 .

Then from (10)–(15) and (37)–(42), we have

(56) I 3 = 0 T ( z ( U h ) z , u U h ) d t = 0 T ( ( y ( U h ) y , y y ( U h ) ) + ( p ( U h ) p , p p ( U h ) ) ) d t 0 .

It follows from (50) and (54)–(56), we obtain (49).□

In order to estimate the error Z ˜ h z ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 , we need the following well-known stability results (see [36] for the details) for the following dual equations:

(57) φ t + b φ div ( A φ ) + c φ = F , x Ω , t J , φ Ω = 0 , t J , φ ( x , 0 ) = 0 , x Ω ,

and

(58) ψ t b ψ div ( A ψ ) + c ψ = F , x Ω , t J , ψ Ω = 0 , t J , ψ ( x , T ) = 0 , x Ω .

## Lemma 3.1

[36] Let φ and ψ be the solutions of (57) and (58), respectively. Let Ω be a convex domain. Then, for ϕ = φ or ϕ = ψ

0 T Ω ϕ ( x , t ) 2 d x C F L 2 ( J ; L 2 ( Ω ) ) 2 , 0 T Ω ϕ 2 d x d t C F L 2 ( J ; L 2 ( Ω ) ) 2 , 0 T Ω D 2 ϕ 2 d x d t C F L 2 ( J ; L 2 ( Ω ) ) 2 , 0 T Ω ϕ t 2 d x d t C F L 2 ( J ; L 2 ( Ω ) ) 2 ,

where D 2 ϕ = max { 2 ϕ / x i x j , 1 i , j 2 } .

## Lemma 3.2

[37] Let f and g be piecewise continuous nonnegative functions defined on 0 t T , g being nondecreasing. If for each t J ,

(59) f ( t ) g ( t ) + 0 t f ( s ) d s ,

then f ( t ) e t g ( t ) .

In the following two theorems, we shall estimate the error Z ˜ h z ( U h ) L 2 ( J ; L 2 ( Ω ) ) .

## Theorem 3.2

Let ( Y h , P h , Z h , Q h , U h ) and ( y ( U h ) , p ( U h ) , z ( U h ) , q ( U h ) , U h ) be the solutions of (30)–(36) and (37)–(42), respectively. Then we have

(60) Y h y ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 C i = 2 8 η i 2 ,

where

η 2 2 = 0 T τ h τ 2 τ Y h , t + div P ˆ h + c Y ˆ h f ˆ U h + Y ˜ h Y ˜ h ( x ¯ ) k 2 d x d t ; η 3 2 = 0 T Ω Y h , t b A 1 P h Y ˆ h Y ˜ h ( x ¯ ) k 2 d x d t ; η 4 2 = 0 T τ h τ 2 τ ( A 1 P h ) 2 d x d t ; η 5 2 = P ˆ h P h L 2 ( J ; L 2 ( Ω ) ) 2 ; η 6 2 = f ˆ f L 2 ( J ; L 2 ( Ω ) ) 2 ; η 7 2 = Y ˆ h Y h L 2 ( J ; L 2 ( Ω ) ) 2 ; η 8 2 = y 0 h ( x ) y 0 ( x ) L 2 ( Ω ) 2 .

## Proof

From (31), we can obtain the following equality:

(61) ( A 1 P h , v h ) ( Y h , div v h ) = 0 , v h V h .

Let ψ be the solution of (58) with F = Y h y ( U h ) , using (30)–(32), (37)–(39), and (45)–(48), we infer that

(62) Y h y ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 = 0 T ( Y h y ( U h ) , F ) d t = 0 T ( Y h y ( U h ) , ψ t div ( A ψ ) b ψ + c ψ ) d t = 0 T ( ( ( Y h y ( U h ) ) t , ψ ) ( Y h , div ( Π h ( A ψ ) ) ) + ( p ( U h ) , ψ ) ) d t + 0 T ( ( b ( y ( U h ) Y h ) , ψ ) + ( c ( Y h y ( U h ) ) , ψ ) ) d t + ( ( Y h y ( U h ) ) ( x , 0 ) , ψ ( x , 0 ) ) = 0 T ( ( ( Y h y ( U h ) ) t , ψ ) ( A 1 P h , Π h ( A ψ ) ) ( b y ( U h ) , ψ ) ( div p ( U h ) , ψ ) ) d t + 0 T ( ( Y h , ( b ψ ) ) + ( c ( Y h y ( U h ) ) , ψ ) ) d t + ( y 0 h ( x ) y 0 ( x ) , ψ ( x , 0 ) ) = 0 T ( ( Y h , t , ψ ) + ( A 1 P h , A ψ Π h ( A ψ ) ) ( P ˆ h P h , ψ ) ( div P ˆ h , ψ ) ) d t + 0 T ( ( Y h , ( Π h ( b ψ ) ) + ( c Y h f U h , ψ ) ) d t + ( y 0 h ( x ) y 0 ( x ) , ψ ( x , 0 ) ) = 0 T Y h , t + div P ˆ h + c Y ˆ h f ˆ U h + Y ˜ h Y ˜ h ( x ¯ ) k , ψ d t + 0 T ( ( f ˆ f , ψ ) + ( c ( Y h Y ˆ h ) , ψ ) + ( P ˆ h P h , ψ ) ) d t + 0 T ( A 1 P h , A ψ Π h ( A ψ ) ) d t + ( y 0 h ( x ) y 0 ( x ) , ψ ( x , 0 ) ) 0 T ( Y h , ( Π h ( b ψ ) ) ) d t 0 T Y ˜ h Y ˜ h ( x ¯ ) k , ψ d t L 1 + L 2 + + L 6 .

Using (30), (43), Cauchy inequality, and Lemma 3.1, we have

(63) L 1 = 0 T Y h , t + div P ˆ h + c Y ˆ h f ˆ U h + Y ˜ h Y ˜ h ( x ¯ ) k , ψ R h ψ d t C ( δ ) η 2 2 + δ ψ L 2 ( J ; H 1 ( Ω ) ) 2 C η 2 2 + 1 6 Y h y ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 .

Similarly, using Cauchy inequality, and Lemma 3.1, we derive

(64) L 2 C ( η 5 2 + η 6 2 + η 7 2 ) + 1 6 Y h y ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 ,

(65) L 3 C η 4 2 + 1 6 Y h y ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 ,

(66) L 4 C η 8 2 + 1 6 Y h y ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 .

Finally, for L 5 and L 6 , using (61), Cauchy inequality, and Lemma 3.1, we find that

(67) L 5 + L 6 = 0 T ( A 1 P h , Π h ( b ψ ) ) d t 0 T Y ˜ h Y ˜ h ( x ¯ ) k , ψ d t = 0 T Y h , t Y ˆ h Y ˜ h ( x ¯ ) k b A 1 P h , ψ d t + 0 T ( A 1 P h , b ψ Π h ( b ψ ) ) d t C η 3 2 + C η 4 2 + 1 6 Y h y ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 .

Hence, using (61)–(67), we obtain (60).□

## Theorem 3.3

Let ( y , p , z , q , u ) and ( Y h , P h , Z h , Q h , U h ) be the solutions of (10)–(16) and (30)–(36), respectively. Let ( y ( U h ) , p ( U h ) , z ( U h ) , q ( U h ) , U h ) be defined as in (37)–(42). Then we have the following error estimate:

(68) Z ˜ h z ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 C η 4 2 + η 7 2 + i = 9 16 η i 2 + C Y h y ( U h ) L 2 ( J ; L 2 ( Ω ) ) 2 ,

where

η 9 2 = 0 T τ h τ 2 τ Z h , t + div Q ˜ h Z ˆ h ( x ¯ ¯ ) J ˆ Z ˆ h k + c Z ˜ h Y ˆ h + y ˆ d 2 d x d t ; η 10 2 = 0 T Ω Z h , t b A 1 Q h b P ¯ h + b p ¯ d Z ˆ h ( x ¯ ¯ ) J ˆ Z ˜ h k 2 d x d t ; η 11 2 = 0 T τ h τ 2 τ ( A 1 Q h + P ¯ h