A mass-conservative characteristic FE scheme for optimal control problems governed by convection–diffusion equations

Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 82–92

Contents lists available at SciVerse ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

A mass-conservative characteristic FE scheme for optimal control problems governed by convection–diffusion equations Hongfei Fu a, Hongxing Rui b,⇑ a b

Department of Computational and Applied Mathematics, China University of Petroleum, Qingdao 266580, China School of Mathematics, Shandong University, Jinan 250100, China

a r t i c l e

i n f o

Article history: Received 17 October 2009 Received in revised form 9 February 2011 Accepted 26 May 2012 Available online 15 June 2012 Keywords: Mass-conservative characteristic finite element Optimal control problems Convection–diffusion equations A priori error estimates

a b s t r a c t In this paper, we develop a mass-conservative characteristic finite element scheme for optimal control problems governed by linear singularly perturbed, convection–diffusion equations. We discuss the case that the velocity fields are non-divergence-free. The space discretization of the state variable is done by piecewise linear continuous functions, whereas the control variable is approximated by piecewise constant functions due to the limitation of the regularity. The scheme preserves the mass balance for the original state equation. We derive a priori error estimates for both the control and state approximations. Some numerical examples are presented to show the efficiency of the proposed scheme. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Optimal control problems governed by convection–diffusion equations arise in many scientific and engineering applications, such as atmospheric pollution control problems [1,2]. Efficient numerical methods are essential to successful applications of such optimal control problems. Although several well-established techniques have been proposed and analyzed, for example, the streamline diffusion finite element method [3], the residual-free bubble method [4], the discontinuous Galerkin method [5], the edge-stabilization Galerkin method [6], and the characteristic method [7–13], it is still a hard and challenging work for the theoretical analysis and numerical simulation of convection-dominated diffusion equations. It is also not an easy thing to use these techniques straightly for such optimal control problems in many cases. To our best knowledge, there are only a few published results on optimal control problems governed by steady convection–diffusion equations; see [14,15] of SUPG method, [16] of the standard finite element discretizations with stabilization, [17] of symmetric stabilization method, [18] of edge-stabilization method, and [19] of the application of RT mixed DG scheme. For the approximation of constrained optimal control problems governed by time-dependent convection–diffusion equations, it is much more complicated and there are nearly no published papers. ⇑ Corresponding author. E-mail addresses: [email protected] (H. Fu), [email protected] (H. Rui). 0045-7825/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2012.05.019

In [20], a class of quadratic optimal control problems governed by linear advection-dominated diffusion equations are discussed with the application of the conventional modified methods of characteristic (MMOC), where the divergence-free velocity fields are dealt with. It is well known that mass conservation is an important property for the approximation scheme. However, in the framework of MMOC it is not trivial to maintain this property. Based on the idea proposed in [21], in this work we present a new characteristic finite element scheme which preserves the mass balance to the quadratic optimal control problems governed by singularly perturbed, two dimensional convection–diffusion state equations. We here do not assume the velocity is incompressible. We obtain a priori error estimates for both the control and state approximations. The main differences of the present work with [20] can be summarized as follows. First, a fully different mass-conservative characteristic (MCC) finite element scheme is adopted. Second, we have deleted the divergence-free limitation on the velocity fields of singularly perturbed convection–diffusion equations. Third, the Robin boundary condition is considered. Forth, the error estimates are obtained in the framework of eweighted energy norm for the state variable, and L2 -norm for the control variable. The rest of the paper is organized as follows. In Section 2, we revisit the optimal control problems and then derive the optimality conditions. In Section 3, we construct a mass-conservative numerical scheme with the characteristic finite element approximation. In Section 4 and Section 5, we establish a e-uniform priori error estimates for the model problem and auxiliary lemmas,

83

H. Fu, H. Rui / Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 82–92

respectively. In Section 6, some numerical examples are presented to observe the convergence behavior of the MCC scheme, we also check the mass error at the final time T for the state variable. We conclude with some further comments in Section 7. Throughout this work, we employ the usual notion for Lebesgue and Sobolev spaces, see [22] for details. In addition, we use C to denote a generic constant which is independent of the discrete parameters and may have different values in different circumstances.

!1p

Dtkf i kpX

Z

yðx; tÞ ¼

Z

y0 ðxÞ þ

X

Z

t

Z

0

ðf þ BuÞðx; tÞ þ

Z 0

X

t

Z

gðx; tÞ:

ð2:5Þ

@X

In the next section, we shall prove that the mass conservative relation (2.5) is also preserved for the MCC method in the discretization form. For a given control set K, the above convex optimal control problem can be restated as follows:

16p<1

i¼1

and the standard modification for p ¼ 1. Let

n o p l ð0; T; XÞ :¼ f : kf klp ð0;T;XÞ < 1 ;

ð2:4aÞ ð2:4bÞ

Substituting w ¼ 1 2 V in (2.4a) and integrating the time from 0 to t ð6 TÞ, we can easily derive the mass balance identity

X

;

8w 2 V;

yðx; 0Þ ¼ y0 ðxÞ:

Let X and XU be bounded open convex polygons in R2 , with Lipschitz boundaries @ X and @ XU . Let ð0; T be the time interval. Partition ð0; T by ti ¼ iDt for 0 6 i 6 N T with Dt :¼ T=N T . Let f i ¼ f ð; ti Þ. We define the discrete time-dependent norms

kf klp ð0;T;XÞ ¼

ðyt ; wÞ  ðv  rw; yÞ þ aðy; wÞ ¼ ðf þ Bu; wÞ þ hg; wi;

2. Optimal control problems and optimality conditions

NT X

Here we extent our discussion to the case with non-divergence-free velocity field, i.e., we delete the assumption r  v ¼ 0. Let ð; Þ and ð; ÞU denote the inner product of X and XU , respectively, and let h; i be the inner product of @ X. Denote aðv ; wÞ ¼ ðerv ; rwÞ. Then a possible weak formula for the state equation reads as follows: For given f ; u and y0 , find y 2 H1 ð0; T; L2 ðXÞÞ \ W such that

1 6 p 6 1:

1 2

Z

T

  kyðuÞ  zd k2 þ akuk2U dt;

Now we give a brief description of the mathematical model of the optimal control problems governed by singularly perturbed convection-dominated diffusion equations. To fix the idea, we shall take the state space W ¼ L2 ð0; T; VÞ with V ¼ H1 ðXÞ, the control space X ¼ L2 ð0; T; UÞ with U ¼ L2 ðXU Þ, and the observation space Y ¼ L2 ð0; T; HÞ with H ¼ L2 ðXÞ. Let B be a linear continuous operator from U to H, and let K be a closed convex set in X defined as follows:

min

K ¼ fv 2 X : n1 6 v 6 n2 ; a:e: in XU  ð0; Tg;

It is well known (see, e.g., [23]) that the control problem (2.6) has a unique solution ðy; uÞ 2 H1 ð0; T; L2 ðXÞÞ \ W  X, and that a pair ðy; uÞ is the solution of (2.6) if and only if there is a co-state p 2 H1 ð0; T; L2 ðXÞÞ \ W such that the triplet ðy; p; uÞ satisfies the following optimality conditions:

where the bounds n1 ; n2 2 R fulfill n1 < n2 . We consider the following distributed optimal control problem:

min u2K

1 2

Z

T

0



 ky  zd k2 þ akuk2U dt

ð2:1Þ

u2K

ðyt ðuÞ; wÞ  ðv  rw; yðuÞÞ þ aðyðuÞ; wÞ ¼ ðf þ Bu; wÞ þ hg; wi;

on @ X  ð0; TÞ;



ðyt ; wÞ  ðv  rw; yÞ þ aðy; wÞ ¼ ðf þ Bu; wÞ þ hg; wi;

8w 2 V;

yð0Þ ¼ y0 ; ð2:7Þ ð2:2aÞ



ð2:2bÞ ð2:2cÞ

ðpt ; qÞ  ðv  rp; qÞ þ aðq; pÞ ¼ ðy  zd ; qÞ;

8q 2 V;

pðTÞ ¼ 0; Z 0

yðx; 0Þ ¼ y0 ðxÞ in X:

8w 2 V;

yðuÞðx; 0Þ ¼ y0 ðxÞ:

combined with the following Robin boundary condition and initial condition

ðery  yv Þ  n ¼ g

ð2:6Þ

where yðuÞ 2 W subject to

subject to the convection–diffusion equation

yt þ r  ðv y  eryÞ ¼ f þ Bu in X  ð0; TÞ

0

T

ðau þ B p; v  uÞU dt P 0;

8v 2 K  X ¼ L2 ð0; T; UÞ;

ð2:8Þ

ð2:9Þ

where B is the adjoint operator of B.

We denote by L the convection–diffusion operator

Ly :¼ yt þ r  ðv y  eryÞ ¼ yt þ v  ry þ ðr  v Þy  eDy:

3. A mass-conservative characteristic finite element method

Here x ¼ ðx1 ; x2 Þ; yðx; tÞ is called state, which is e-dependent; uðx; tÞ is called control, which is usually constrained; a is a positive constant; 0 < e  1 is a small constant that scales the diffusion and characterizes the convection-dominance of (2.2); we assume the velocity v ¼ ðV 1 ðx; tÞ; V 2 ðx; tÞÞ satisfy the following conditions, which makes the following argument simple and clear. Hypothesis 2.1. The velocity v satisfies

In this section, we will give the mass-conservative characteristic (MCC) finite element approximation for the control problem (2.6). The approximation scheme is also applicable to the control problems with more general convex objective functionals. Here for simplicity we consider the n-simplex conforming Lagrange element, which is most widely used in practical computations. Let z ¼ Gðx ; t ; tÞ be an approximate characteristic curve passing through point x at time t , which is defined by

v 2 Cð0; T; W 1;1 ðXÞÞ2 ;

Gðx ; t ; tÞ :¼ x  v ðx ; t Þðt   tÞ:

v ¼ 0 on @ X:

In a previous paper [20], we have discussed the model control problem (2.1) with the restriction that the velocity v is incompressible, i.e.,

r  v ¼ 0;

8ðx; tÞ 2 X  ð0; TÞ:

ð2:3Þ

ð3:1Þ

We denote by  x be the foot at time ti1 of the characteristic curve with head x at time t i , and x represents the head of the characteristic curve with foot x at time t i1 , namely,

x ¼ Gðx; t i ; ti1 Þ;

x ¼ Gðx; t i ; t i1 Þ:

ð3:2Þ

84

H. Fu, H. Rui / Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 82–92

Remark 3.1. [10] Under Hypothesis 2.1 and Dt < 1= kvkCð0;T;W 1;1 ðXÞÞ2 ; G is a differential homeomorphism of X onto itself. Remark 3.1 guarantees that the tangent to the characteristic curve can not cross a boundary to an undefined location. The conventional characteristic methods, such as MMOC, usually combine the time derivative term and the convection term, i.e., the material derivative term

L0 y ¼ yt þ v  ry

ð3:3Þ

in the governing equations to carry out the temporal discretization in a Lagrange coordinate. These methods symmetrize the governing equation and stabilize their numerical approximations. They generate accurate numerical solutions and significantly reduce the numerical diffusion and grid-orientation effect present in upwind methods, even if large time steps and coarse spatial meshes are used. Here our idea is to directly apply the characteristic approximation to the term

L1 y ¼ yt þ v  ry þ ðr  v Þy:

ð3:4Þ

We set by J i the Jacobian determinant of the transformation from x to G that

J i ¼ det

    @Gðx; ti ; t i1 Þ ¼ 1  ðr  v i ÞDt þ O Dt 2 ¼ 1 þ OðDtÞ: @x ð3:5Þ

Then for sufficiently smooth state y, a first-order approximation to the term L1 y can be derived following from (3.5) that

L1 y ¼ ðyt þ v  ryÞjt¼ti i

Z X

ynh ¼

Z X

y0h þ

Z  n  Z X Dt f i þ Buih þ i¼1

ð3:9Þ

@X

Remark 3.3. Theorem 3.2 can be easily proved by choosing wh ¼ 1 in (3.8), multiplying Dt, and summing over i from 1 to n. For the detailed proof, we refer the reader to [21]. We say here the mass balance relation (3.9) is the discrete version of (2.5) and it is in the sense of discrete state and control variables. Since (3.8) has a unique solution yih for every uih , and (3.7) is convex in yih and is strictly convex in uih , we then conclude that the control problem (3.7) has a unique solution ðY ih ; U ih Þ; i ¼ 1; 2;    ; N T , and that a pair ðY ih ; U ih Þ 2 V h  K h ; i ¼ 1; 2;    ; N T , is the solution of (3.7) if and only if there is a co-state P hi1 2 V h ; i ¼ i h h h 1; 2;    ; N T , such that the triplet ðY ih ; Pi1 h ; Uh Þ 2 V  V  K ; i ¼ 1; 2;    ; N T , satisfies the following optimality conditions:

 8  i i1 i   Y h Y h J > > þ a Y ih ; wh ¼ ðf i þ BU ih ; wh Þ þ hg i ; wh i; ; w > h Dt < > 8wh 2 V h ; i ¼ 1; 2; . . . ; NT ; > > : 0 Y h ðxÞ ¼ yh0 ðxÞ;  8  i1 i   P h P h i1 > > þ a q ; q ; P ¼ ðY ih  zid ; qh Þ; > h h h Dt <



 U

ð3:10Þ

ð3:11Þ

> 8qh 2 V h ; i ¼ NT ; . . . ; 2; 1; > > : NT Ph ðxÞ ¼ 0; i aU ih þ B Pi1 h ; v h  Uh

i1 þ OðDtÞ þ ðr  v Þy i

X

 gi :

P 0; 8v h 2 K h  K \ U h ; i ¼ 1;2;...;N T ; ð3:12Þ

i1 y i1  J i i1  y yi  y þ þ OðDtÞ Dt Dt i1  J i yi  y þ OðDtÞ; ¼ Dt

¼

 i ¼ P i ð  where P h xÞ, and x is defined by (3.2). h

ð3:6Þ

i

 ¼ yðxÞ. In the following we shall point that where y ¼ yð; t i Þ and y the new scheme is mass balance. Let T h and T hU be regular triangulations of X and XU ,  ¼ [ hs   respectively, X s2T , XU ¼ [sU 2T hU sU . Let h ¼ maxs2T h hs ; hU ¼ maxsU 2T hU hsU , where hs and hsU denote the diameter of the element s and sU , respectively.  Þ, Associate with T h is a finite dimensional subspace V h of CðX such that vjs are first-order polynomials for all v 2 V h and s 2 T h . It is clear that V h  V. Associate with T hU is another finite dimensional subspace U h of L2 ðXU Þ, such that vjsU are constants for all v 2 U h and sU 2 T hU . Here there is no requirement for the continuity. A fully discrete approximate scheme of optimal control problem (2.6) is to find ðyih ; uih Þ 2 V h  K h ; i ¼ 1; 2;    ; N T , such that

  T 1X Dt kyih  zd ðx; t i Þk2 þ akuih k2U ui 2K h 2 i¼1 h N

min

ð3:7Þ

Remark 3.4. The MMOC which approximates the material derivative term in (3.3) leads to another discrete optimality conditions:

  8      Y ih Y hi1 > > þ r  vi Y ih ; wh þ a Y ih ; wh ; w > h D t > > > >   < ¼ f i þ BU ih ; wh þ hg i ; wh i; > > > > > 8wh 2 V h ; i ¼ 1; 2; . . . ; NT ; > > : 0 Y h ðxÞ ¼ yh0 ðxÞ;

ð3:13Þ

  8  i1       P h P ih b Ji > > þ r  vi qh ; Pi1 þ a qh ; Pi1 ¼ ðY ih  zid ; qh Þ; ; q > h h h Dt < > 8q 2 V h ; i ¼ NT ; . . . ; 2; 1; > > : NTh Ph ðxÞ ¼ 0 ð3:14Þ combined with the variational inequality (3.12), where bJ i is the reciprocal of Ji . Comparing to (3.13)–(3.14), the present MCC scheme (3.10)–(3.11) is simpler, symmetric and mass-conservative.

subject to i i1 yih  y h J ; wh Dt

!

  þ a yih ; wh

¼ ðf i þ Buih ; wh Þ þ hg i ; wh i;

4. A priori error estimates

8wh 2 V h ;

y0h ðxÞ ¼ yh0 ðxÞ;

ð3:8Þ

where K h is a closed convex set in U h , and yh0 2 V h is an approximation of y0 which will be specified later on. For ease of exposition, we also assume K h  K \ U h . yih

Theorem 3.2 (Mass balance). Let ði ¼ 1; 2; . . . ; N T Þ be the solution of (3.8). Under Hypothesis 2.1 and the condition Dt < 1= kvkCð0;T;W 1;1 ðXÞÞ2 , the following conservation holds

In this section, we are able to derive some a priori error esti2 mates in l ð0; T; L2 ðXU ÞÞ-norm for the control variable and in le ð0; T; H1 ðXÞÞ-norm for the state variable. Here le ð0; T; H1 ðXÞÞ is introduced by using the following e-weighted energy norm

kf kle ð0;T;H1 ðXÞÞ :¼

kf k2l1 ð0;T;L2 ðXÞÞ þ

NT X i¼1

where kf i k2a ¼ aðf i ; f i Þ.

!1=2

Dtkf i k2a

;

85

H. Fu, H. Rui / Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 82–92

To obtain a priori error estimates for the scheme (3.10)–(3.12), we introduce two auxiliary variables ðY ih ðuÞ; P ih ðuÞÞ 2 V h  V h ; i ¼ 1; 2;    ; N T , associate with the control u as follows:

Next, we give an estimate for the control approximation, it is proved to rely on the co-state approximation.

 8 i   Y h ðuÞY i1 ðuÞJ i > h > þ a Y ih ðuÞ; wh ; w h > Dt <

Lemma 4.2. Let ðy; p; uÞ and ðY h ; P h ; U h Þ be the solutions of (2.7)– (2.9) and (3.10)–(3.12), respectively. Assume that u 2 Cð0; T; H1 ðXU ÞÞ; p 2 Cð0; T; H1 ðXÞÞ \ H1 ð0; T; L2 ðXÞÞ; K h  K \ U h , and there exists Ph ui 2 K h ð1 6 i 6 N T Þ such that

> > > :

¼ ðf i þ Bui ; wh Þ þ hg i ; wh i;

ð4:1Þ

8wh 2 V h ;

Y 0h ðuÞ ¼ yh0 ðxÞ;



8    i i1 > < Ph ðuÞPh ðuÞ ; q þ a q ; Pi1 ðuÞ ¼ ðY i ðuÞ  zi ; q Þ; h h h h h d Dt > : NT Ph ðuÞ ¼ 0:

8qh 2 V h ;

For simplicity of illustration, in the following we use the notation

gi ¼ yi  Y ih ðuÞ; i ¼ 0; 1; . . . ; NT ;

fi ¼ Pih  P ih ðuÞ;

ni ¼ pi  Pih ðuÞ;



2

U

6 ChU ;

ð4:11Þ

kui  Ph ui kL2 ðXU Þ 6 ChU kui kH1 ðXU Þ :   ku  U h kl2 ð0;T;L2 ðXU ÞÞ 6 C hU þ Dt þ kp  P h ðuÞkl1 ð0;T;L2 ðXÞÞ ;

It is clear that h0 ¼ 0 and fNT ¼ 0.

where Ph ðuÞ is defined in (4.2).

aku  U h k2l2 ð0;T;L2 ðXU ÞÞ ¼

NT   X Dt aui ; ui  U ih

U

i¼1

Lemma 4.1. Let ðY h ; Ph Þ and ðY h ðuÞ; P h ðuÞÞ be the solutions of (3.10)–(3.11) and (4.1)–(4.2), respectively. There hold the following estimates:

!

i1

h h Dt

;wh

i¼1



ð4:4Þ

þ

Noticing h 2 V , select wh ¼ h as a test function in (4.5). The 2 inequality aða  bÞ P 12 ða2  b Þ and a direct calculation shows that

hi  hi1 i ;h Dt

!

 1  i 2 P kh k  khi1 k2 ; 2Dt

ð4:6Þ

NT   X i i Dt aU ih þ B Pi1 h ðuÞ; u  U h

i i i Dt B ðPi1 h ðuÞ  p Þ; u  U h

6

þ

NT   X Dt B ðPhi1 ðuÞ  pi Þ; ui  U ih :

Note that Ph ui 2 K , we then conclude from (3.7) that

aku  U h k2l2 ð0;T;L2 ðXU ÞÞ 6

NT   X i i Dt aU ih þ B P i1 h ; U h  Ph u NT X

i i Dt aU ih þ B Pi1 h ; Ph u  u

since h0 ¼ 0. Thus (4.3) follows from (4.9) and the discrete Gronwall’s lemma. Similarly, we can prove from Eqs. (3.11) and (4.2) that

Then (4.4) follows from (4.10) and (4.3).

ð4:10Þ h

U

  i1 i i Dt B ðPi1 h  P h ðuÞÞ; u  U h

NT X

U

NT     X i i i Dt B Pi1 h ðuÞ  p ; u  U h

U

i¼1

6

NT X



Dt aðU ih  ui Þ; Ph ui  ui

i¼1

þ

 U

NT X   Dt aui þ B pi ; Ph ui  ui U i¼1

þ

NT X   Dt B ðpi1  pi Þ; Ph ui  ui U i¼1

NT   X i i þ Dt B ðpi1  Pi1 h ðuÞÞ; u  Ph u i¼1

kfkle ð0;T;H1 ðXÞÞ 6 CkY h  Y h ðuÞkl2 ð0;T;L2 ðXÞÞ :



i¼1

i¼1

ð4:9Þ



i¼1

þ

i¼1

U

i¼1

 1  i 2 kh k  khi1 k2 þ khi k2a 6 Ckui  U ih k2U þ Ckhi k2 þ Ckhi1 k2 : 2 Dt ð4:8Þ

i¼1

U

i¼1 h

þ

N N N X X X 2Dtkhi k2a 6 C Dtðkhi k2 þ khi1 k2 Þ þ C Dtkui  U ih k2U

U

U

i¼1

Then it follows from Cauchy–Schwarz inequality, the continuous property of B and the fact Ji ¼ 1 þ OðDtÞ that

khN k2 þ



NT   X i i Dt aU ih þ B P i1 h ðuÞ; U h  u

ð4:7Þ

Multiply both sides of (4.8) by 2Dt and sum over i from 1 to N, we obtain

U



NT X

þ

khi1 k2 6 ð1 þ C DtÞkhi1 k2 :

U

i¼1

ð4:5Þ i

U

i¼1

!   hi1  hi1  J i i i i ; wh : þ a h ;wh ¼ ðBðU h  u Þ;wh Þ  Dt

h

i



NT   X ¼ Dt aui þ B pi ; ui  U ih

Proof. We first prove (4.3). It follows from Eqs. (4.1) and (3.10) that i



T  sumNi¼1 Dt aU ih ; ui  U ih

ð4:3Þ

kP h  Ph ðuÞkle ð0;T;H1 ðXÞÞ 6 Cku  U h kl2 ð0;T;L2 ðXU ÞÞ :

ð4:13Þ

Proof. It follows from (2.9) that

i ¼ NT ; . . . ; 1; 0:

kY h  Y h ðuÞkle ð0;T;H1 ðXÞÞ 6 Cku  U h kl2 ð0;T;L2 ðXU ÞÞ ;

ð4:12Þ

Then there exists a constant C such that

ð4:2Þ

hi ¼ Y ih  Y ih ðuÞ;

aui þ B pi ; Ph ui  ui

þ

U

NT   X i1 i i Dt B ðPi1 h ðuÞ  P h Þ; u  Ph u i¼1

U

86

H. Fu, H. Rui / Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 82–92

þ

NT   X i1 i i Dt B ðPi1 h  P h ðuÞÞ; u  U h

The following lemma describes the estimates between the exact solution ðy; pÞ and the auxiliary solution ðY h ðuÞ; Ph ðuÞÞ.

U

i¼1 NT   X i1 þ Dt B ðPi1 Þ; ui  U ih h ðuÞ  p

U

i¼1

þ

NT   X Dt B ðpi1  pi Þ; ui  U ih : U

i¼1 0

ð4:14Þ

NT

Noticing h ¼ f ¼ 0, it follows from (3.10) and (3.11) to (4.1) and (4.2) that for the sixth term on the right-hand side of (4.14) NT     X i1 i i Dt B Pi1 h  P h ðuÞ ; u  U h

¼

e

ð4:15Þ

e

ð4:16Þ

!

NT NT   X X hi  hi1  J i i1 Dt Dta hi ; fi1 ;f  Dt i¼1 i¼1 ! N NT NT i1   T X X X f  fi i Dta hi ; fi1 ¼  Dtðhi ; hi Þ ¼  Dt ;h  Dt i¼1 i¼1 i¼1

¼

¼ khk2l2 ð0;T;L2 ðXÞÞ 6 0:

where C depends on some spatial and temporal derivatives of y; p and zd . Proof. First, we give an estimate for the difference g between y and Y h ðuÞ. From Eq. (2.7), we can easily see that the exact solution y satisfies for i P 1

! i1  J i yi  y ;wh þ aðyi ;wh Þ ¼ ðf i þ Bui ; wh Þ  ðri ;wh Þ; 8wh 2 V h ; Dt

Then we obtain from Lemma 4.1, Cauchy-Schawarz inequality, and the above inequality that

aku  U h k2l2 ð0;T;L2 ðXU ÞÞ 6

NT X

  Dt aui þ B pi ; Ph ui  ui U þ CðdÞku

i¼1

P

2 h ukl2 ð0;T;L2 ðXU ÞÞ

!

2

NT   X i1 i i Dt Pi1 h  P h ðuÞ; Bðu  U h Þ i¼1

! 2 h 2 ky  Y h ðuÞkle ð0;T;H1 ðXÞÞ 6 C h þ minfh; Dtg þ pffiffiffi þ he1=2 þ Dt ;

h 2 kp  Ph ðuÞkle ð0;T;H1 ðXÞÞ 6 C h þ minfh; Dtg þ pffiffiffi þ he1=2 þ Dt ;

U

i¼1

Lemma 4.4. Let ðy; pÞ and ðY h ðuÞ; P h ðuÞÞ be the solutions of (2.7)– 1 (2.8) and (4.1)–(4.2), respectively. Assume that y; p 2 l ð0; T; 2 1 2 2 2 1 2 H ðXÞÞ \ H ð0; T; H ðXÞÞ \ H ð0; T; L ðXÞÞ; zd 2 H ð0; T; L ðXÞÞ. Then the following estimates hold

2

þ CðdÞðDtÞ

kpt k2L2 ð0;T;L2 ðXÞÞ

þ CðdÞkp  Ph ðuÞk2l2 ð0;T;L2 ðXÞÞ þ CdkP h ðuÞ  P h k2l2 ð0;T;L2 ðXÞÞ þ Cdku  U h k2l2 ð0;T;L2 ðXU ÞÞ   2 6 ChU kuk2Cð0;T;H1 ðXU ÞÞ þ kpk2Cð0;T;H1 ðXÞ

ð4:17Þ where

ri ¼

i1  J i @yi yi  y þ r  ðv i yi Þ  : @t Dt

By a standard backward difference error analysis [7,9] and the results in [21], we have

kri k2 6 C Dtkyss k2L2 ðti1 ;ti ;L2 ðXÞÞ þ C Dt kyt k2L2 ðti1 ;ti ;L2 ðXÞÞ þ CðDtÞ2 kyi1 k2H1 ðXÞ :

þ CðDtÞ2 kpt k2L2 ð0;T;L2 ðXÞÞ þ Ckp 

P h ðuÞk2l2 ð0;T;L2 ðXÞÞ

þ Cdku 

U h k2l2 ð0;T;L2 ðXU ÞÞ ;

ð4:18Þ Subtracting Eq. (4.1) from Eq. (4.17), we obtain an error equation on g ¼ y  Y h ðuÞ:

where and after, d is an arbitrary small positive number and CðdÞ depends on 1=d. We finish the proof by choosing Cd ¼ a2. h Remark 4.3. Since the control variable is approximated by piecewise constant functions, we take

Ph v jsU ¼

1

Z

j s U j sU

v ; 8sU 2 T

h U;

gi  g i1  Ji Dt

!

; wh

  þ a gi ; wh ¼ ðri ; wh Þ;

8wh 2 V h ;

i P 1; ð4:19Þ

where g ¼ y  Y h ðuÞ can be decomposed as g ¼ ðy  HyÞ þ ðHy Y h ðuÞÞ ¼ l þ m, where HyðtÞ 2 V h is defined to be the elliptic projection of yðtÞ 2 V such that

8wh 2 V h ;

8t 2 ð0; T:

where jsU j is the measure of sU . It is easy to see that Ph ui 2 K h , and it follows from [22] that for u 2 Cð0; T; H1 ðXU ÞÞ

aðyðtÞ  HyðtÞ; wh Þ ¼ 0;

ku  Ph ukCð0;T;L2 ðXU ÞÞ 6 ChU kukCð0;T;H1 ðXU ÞÞ :

8 < ky  HykHk ðXÞ 6 Ch2k kykH2 ðXÞ ;   : kðy  HyÞ k k 6 Ch2k kyk 2 þ ky k 2 t H ðXÞ ; t H ðXÞ H ðXÞ

ð4:21Þ

  2 ky  Hykle ð0;T;H1 ðXÞÞ 6 C h þ he1=2 kykl1 ð0;T;H2 ðXÞÞ :

ð4:22Þ

Moreover, if u 2 Cð0; T; H1 ðXU ÞÞ and p 2 Cð0; T; H1 ðXÞÞ, we have NT X   Dt aui þ B pi ; Ph ui  ui U i¼1 NT X X Z   ¼ Dt aui þ B pi  Ph ðaui þ B pi Þ ðPh ui  ui Þ i¼1

sU 2T hU

sU

6 kau þ B p  Ph ðau þ B pÞkCð0;T;L2 ðXU ÞÞ kPh u  ukCð0;T;L2 ðXU ÞÞ   2 6 ChU kukCð0;T;H1 ðXU ÞÞ þ kpkCð0;T;H1 ðXÞÞ kukCð0;T;H1 ðXU ÞÞ   2 2 6 ChU kuk2Cð0;T;H1 ðXU ÞÞ þ kpk2Cð0;T;H1 ðXÞÞ 6 ChU ; which proof the assumptions (4.11) and (4.12).

ð4:20Þ

It is well known that the following estimates hold for k ¼ 0; 1 [22]:

Since the estimate for l is known, we need only to derive the estimate for m. We choose wh ¼ mi and make use of (4.20) to rewrite Eq. (4.19) in terms of l and m





  þ a mi ; mi ! mi1  mi1  Ji i ¼ ; m  ðri ; mi Þ  Dt

mi  mi1 Dt

; mi

li  l i1  Ji Dt

! ; mi :

ð4:23Þ

87

H. Fu, H. Rui / Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 82–92

Multiply both sides of (4.23) by Dt and sum over 1 6 i 6 N. The same manner as (4.6)–(4.8) yields N N N   X X 1 N 2 X km k þ Dtkmi k2a 6 C Dt kmi k2 þ kmi1 k2 þ C Dtkri k2 2 i¼1 i¼1 i¼1

þ

N X  i1  J i ; mi Þj; jðli  l



pi1  p i Dt

ð4:24Þ

i¼1

where we choose yh0 to be the elliptic projection of y0 which satisfies (4.20), i.e., m0 ¼ 0. The estimate of the last term on the right-hand side of (4.24) presents major difficulties. For clarity of exposition, the proof is presented in Lemma 5.1, there we obtain N N X X  i1  J i ; mi Þj 6 C Dtkmi k2 þ C minfh2 ; ðDtÞ2 gkyk2l1 ð0;T;H2 ðXÞÞ jðli  l i¼1

It is clear that the results (4.21) and (4.22) are still valid with p instead of y. By choosing qh ¼ pi1 and utilizing the projection (4.20), Eq. (4.28) can be rewritten in the following form

; pi1



  þ a pi1 ; pi1

¼ ðvi1 ; pi1 Þ  

i1

þ y

i



i1

y ;p

qi1  q i



Dt   i1 þ zid  zi1 : d ;p

NT NT X X 1 M 2 Dtkpi1 k2a 6 C Dtðkpi1 k2 þ kpi k2 Þ kp k þ 2 i¼Mþ1 i¼Mþ1

4

N 1X

2

h

4

e

kyk2l1 ð0;T;H2 ðXÞÞ

Dtkm

ð4:25Þ

2

þ Ch

4

2



h

e

!

2

2

2

kyss kL2 ð0;T;L2 ðXÞÞ þ kyt kL2 ð0;T;L2 ðXÞÞ þ kykl2 ð0;T;H1 ðXÞÞ



2 kyt kL2 ð0;T;H2 ðXÞÞ :

Combining (4.26) with the well known estimate (4.22) for l, we finish the proof of (4.15). Next, we give the estimate for the difference h between p and P h ðuÞ. From Eq. (2.8), we can easily see that the exact solution p satisfies for i 6 N T

i    pi1  p i1 ; qh Þ; 8qh 2 V h ; ; qh þ a qh ; pi1 ¼ ðyi1  zi1 d ; qh Þ  ðv Dt

vi1

NT X

 i ; pi1 Þj 6 C jðqi1  q

i¼Mþ1

4

4

þ Ch kpt k2L2 ð0;T;H2 ðXÞÞ þ

! kpk2l1 ð0;T;H2 ðXÞÞ

e

NT 1 X Dtkpi1 k2a : 2 i¼Mþ1

Collecting the estimates (4.30) and (4.31), and the proved results (4.15) together we conclude NT X

kpM k2 þ

Dtkpi1 k2a

i¼Mþ1

X

6 CðDtÞ2

v ¼y;p

X

4

kv

X

kv ss k2L2 ð0;T;L2 ðXÞÞ þ

2 t kL2 ð0;T;H2 ðXÞÞ þ C

kv k2l1 ð0;T;H2 ðXÞÞ þ C

NT X

!

X v ¼y;zd

kv t k2L2 ð0;T;L2 ðXÞÞ þ kyk2l2 ð0;T;H1 ðXÞÞ 2

4

minfh ;ðDtÞ2 g þ h þ

h

4

!

e

Dtðkpi1 k2 þ kpi k2 Þ:

i¼Mþ1

Thus, it follows from the discrete Gronwall’s lemma for sufficient small Dt that kpkle ð0;T;H1 ðXÞÞ 6 C Dt

¼ ðvi1 ; qh Þ þ ðyi  Y ih ðuÞ; qh Þ þ ðyi1  yi ; qh Þ

8q h 2 V h ;

h

ð4:31Þ

v ¼y;p

  þ a qh ; ni1

þ ðzid  zi1 d ; qh Þ;

4

2

þ C minfh ; ðDtÞ g þ h þ



Thus, Eqs. (4.2) and (4.27) can be differenced to show that

!

ð4:30Þ

Dtkpi1 k2 2

v ¼y;p

i @pi1 pi1  p ¼ /i  : @s Dt

NT X i¼Mþ1

þ Ch

can be written as

ni1  ni ; qh Dt

 i ; pi1 Þj: jðqi1  q

The estimate of the last term on the right-hand side of (4.30) also presents major difficulties. For convenience of explanation, its proof is presented in Lemma 5.2, there we obtain

ð4:27Þ

v

Dtkzdi1  zid k

i¼Mþ1

ð4:26Þ

where

NT X

þ

kyk2l1 ð0;T;H2 ðXÞÞ

e

i1

NT X i¼Mþ1

4

Application of the discrete Gronwall’s lemma for sufficient small Dt yields that   kmkle ð0;T;H1 ðXÞÞ 6 C Dt kyss kL2 ð0;T;L2 ðXÞÞ þ kyt kL2 ð0;T;L2 ðXÞÞ þ kykl2 ð0;T;H1 ðXÞÞ ! 2 h 2 þ Ch kyt kL2 ð0;T;H2 ðXÞÞ þ C minfh; Dtg þ pffiffiffi kykl1 ð0;T;H2 ðXÞÞ :



Dtkyi  yi1 k2

i¼Mþ1

þC

i¼1

þ CðDtÞ

NT X

þC

N N X X k þ Dtkmi k2a 6 C Dtðkmi k2 þ kmi1 k2 Þ

þ C minfh ;ðDtÞ2 g þ

Dtkyi  Y ih ðuÞk2

i¼Mþ1

i¼1

N 2

i¼1

NT X

þC

i 2 ka :

Dtkvi1 k2

i¼Mþ1

We now combine (4.24) with the estimates (4.18) and (4.25) to obtain that km

NT X

þC

þ Ch kyt k2L2 ð0;T;H2 ðXÞÞ þ C

ð4:29Þ

Multiply both sides of (4.29) by Dt and sum over i from N T to M þ 1. We conclude by the same manner as (4.6)–(4.8) that

i¼1

þ

   ; pi1 þ yi  Y ih ðuÞ; pi1

i 6 NT :

ð4:28Þ

Similarly, We decompose the error n ¼ p  Ph ðuÞ as n ¼ ðp  HpÞ þ ðHp  P h ðuÞÞ ¼ q þ p, where Hp 2 V h is defined to be the elliptic projection of p 2 V which also satisfies Eq. (4.20).

X

kv ss kL2 ð0;T;L2 ðXÞÞ þ

v ¼y;p

X 2 kv t kL2 ð0;T;H2 ðXÞÞ þ Ch

X

! kv

v ¼y;zd

2 t kL2 ð0;T;L2 ðXÞÞ

þ kykl2 ð0;T;H1 ðXÞÞ

v ¼y;p

! 2 X h 2 þ C minfh; Dtg þ h þ pffiffiffi kv kl1 ð0;T;H2 ðXÞÞ :

e

v ¼y;p

ð4:32Þ

88

H. Fu, H. Rui / Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 82–92

Incorporating (4.32) with the well known estimate (4.22) for q, we finish the proof of (4.16). Thus Lemma 4.4 is proved. h Combing the bounds given by Lemmas 4.1, and 4.2, we have the following main result.

Case 2. Cr 6 c. This case requires much more attention. Note that Remark 3.1 guarantees that G is one-to-one mapping from X to itself, thus we have

i Þ  i1  J i ; mi Þ ¼ ðli ; mi Þ  ðli1 ; m ðli  l i  mi Þ: ¼ ðli  li1 ; mi Þ  ðli1 ; m

Theorem 4.5. Suppose that fy; p; ug and fY h ; Ph ; U h g are the solutions of (2.7)–(2.9) and (3.10)–(3.12), respectively. Assume that all conditions of Lemmas 4.1, and 4.2 em are valid. Then

ky  Y h kle ð0;T;H1 ðXÞÞ þ kp  P h kle ð0;T;H1 ðXÞÞ þ ku  U h kl2 ð0;T;L2 ðXU ÞÞ ! 2 h 2 1=2 6 C hU þ h þ minfh; Dtg þ pffiffiffi þ he þ Dt ;

e

We apply the estimate (4.21) to bound the first term on the righthand side of (5.4) by

 i  ðl  li1 ; mi Þ ¼

Z "Z

where the constant C depends on some spatial and temporal derivatives of y; p; zd and u, but is independent of hU , h and Dt. Proof. It follows from (4.13) and (4.16) that

ti

t i1

X

ð4:33Þ

#

lt dt mi dx 4

6 C Dtkmi k2 þ Ch kyt k2L2 ðti1 ;ti ;H2 ðXÞÞ :

! 2 h 6 C hU þ h þ minfh; Dtg þ pffiffiffi þ he1=2 þ Dt : 2

e

ð4:34Þ

  @x   ¼ 1 þ OðDtÞ: @z

Z  i1 i     mi Þ ¼  ðl ; m 

  d i m ðx þ sðx  xÞÞds dx d s 0 X  Z Z 1

  i1 ¼  l rmi ðx þ sðx  xÞÞds  ðx  xÞdx

li1

Z

1

0

X

Z Z 6 C Dtkli1 k

Moreover, it follows from Lemmas 4.1, 4.4 and (4.34) that

ky  Y h kle ð0;T;H1 ðXÞÞ þ kp  P h kle ð0;T;H1 ðXÞÞ

X

6 kY h  Y h ðuÞkle ð0;T;H1 ðXÞÞ þ kPh  Ph ðuÞkle ð0;T;H1 ðXÞÞ

1

6 C Dth kyi1 kH2 ðXÞ

þ ky  Y h ðuÞkle ð0;T;H1 ðXÞÞ þ kp  Ph ðuÞkle ð0;T;H1 ðXÞÞ 2

h 2 6 Cku  U h kl2 ð0;T;L2 ðXU ÞÞ þ C h þ minfh; Dtg þ pffiffiffi þ he1=2 þ Dt

!

1=2 jrmi ðx þ sðx  xÞÞj2 dsdx

0

2

Z

1=2 jrmi ðzÞj2 dz

X 4

6

1 h Dtkmi k2a þ C Dt kyk2l1 ð0;T;H2 ðXÞÞ : 2 e

e

!

h 2 6 C hU þ h þ minfh; Dtg þ pffiffiffi þ he1=2 þ Dt :

ð4:35Þ

e

Thus we obtain Theorem 4.5 from (4.34) and (4.35).

ð5:5Þ

Let z ¼ x þ sðx  xÞ, note that

Then we bound the second term on the right-hand side of (5.4) as follows:

  ju  U h kl2 ð0;T;L2 ðXU ÞÞ 6 C hU þ Dt þ kp  Ph ðuÞkl1 ð0;T;L2 ðXÞÞ

2

ð5:4Þ

h

ð5:6Þ Gathering Eqs. (5.3), (5.5), and (5.6) and summing over i from 1 to N, we finish the proof. h 1

5. Auxiliary lemmas

Lemma 5.2. Assume that p 2 l ð0; T; H2 ðXÞÞ \ H1 ð0; T; H2 ðXÞÞ. Let Hp 2 V h be the elliptic projection of p 2 V which satisfies (4.20) and q ¼ p  Hp. Then the estimate (4.31) holds.

In this section, we will prove the two auxiliary lemmas used to prove the error bounds in (4.25) and (4.31), respectively.

Proof. This proof of lemma is basically similar to that of Lemma 5.1. For a given constant 0 < c < 1, when Cr P c, we have

1

Lemma 5.1. Assume that y 2 l ð0; T; H2 ðXÞÞ \ H1 ð0; T; H2 ðXÞÞ. Let Hy 2 V h be the elliptic projection of y 2 V which satisfies (4.20) and l ¼ y  Hy. Then the estimate (4.25) holds.

 i ; pi1 Þ ¼ ðqi1 ; pi1 Þ  ðqi ; p  i1  J i Þ ðqi1  q



jv ðx; tÞjDt : Cr :¼ max ðx;tÞ2X½0;T h

ð5:1Þ

For a given constant 0 < c < 1, the proof is divided into two cases. Case 1. Cr P c which implies

h6

ð5:7Þ For Cr 6 c. This case requires much more attention. From Remark 3.1 we have that

Proof. Denote the Courant number by

C

 i1  ðq  q  i ; pi1 Þ 6 C Dtkpi1 k2 þ C Dt minfh2 ; ðDtÞ2 gkpk21 l ð0;T;H2 ðXÞÞ :

max jv ðx; tÞjDt 6 C Dt:

ð5:2Þ

c ðx;tÞ2X½0;T

Using the estimate (4.7), the facts (3.5) and (5.2), we have

 i1 Þ ¼ ðqi1  qi ; pi1 Þ þ ðqi ; pi1  p  i1  p  i1  J i Þ: þ ðqi ; p

ð5:8Þ

The first and second terms on the right-hand side of (5.8) could be estimated as Case 2 in Lemma 5.1,

 i1  ðq  qi ; pi1 Þ þ ðqi ; pi1  p  i1 Þ 4

6 C Dtkpi1 k2 þ Ch kpt k2L2 ðti1 ;ti ;H2 ðXÞÞ

     i  i1 i i  i ðl  l  J ; m Þ 6 Ckmi k kli k þ kli1  J k

4

1 h þ Dtkpi1 k2a þ C Dt kpk2l1 ð0;T;H2 ðXÞÞ : 2 e

2

6 Ch kmi kkykl1 ð0;T;H2 ðXÞÞ 2

6 C Dtkmi k2 þ C Dt minfh ; ðDtÞ2 gkyk2l1 ð0;T;H2 ðXÞÞ : ð5:3Þ

ð5:9Þ

For the third term on the right-hand side of (5.8), noticing Ji ¼ 1 þ OðDtÞ, we have

89

H. Fu, H. Rui / Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 82–92

   i  i1  i1 i   i1 k ðq ; p  p  J Þ 6 C Dtkqi kkp

Table 2 The results obtained by the MMOC scheme for Example 6.1.

2

6 C Dth kpi kH2 ðXÞ kpi1 k i1 2

6 C Dtkp

k þ C Dth

4

kpk2l1 ð0;T;H2 ðXÞÞ :

ð5:10Þ

Combining the estimates (5.7)–(5.10) and summing over i from N T to M þ 1, we finish the proof. h 6. Numerical experiments In this section, we conduct some numerical examples to check the efficiency of the MCC finite element scheme (3.10)–(3.12). We also compare the numerical results with those of obtained by the MMOC finite element scheme (3.13), and (3.14). Consider the problem:

min

1 2

Z

T



 ky  zd k2 þ ku  u0 k2 dt;

ð6:1Þ

0

1=h

8

16

32

64

Ey ðhÞ EOC Ey ðhÞ

4.2661E02 – 7.2211E02

1.3578E02 1.6516 3.0868E02

4.5767E03 1.5689 1.4609E02

1.8639E03 1.2960 7.1875E03

EOC Ep ðhÞ EOC Ep ðhÞ

– 9.2168E03 – 1.0769E02

1.2261 3.8561E03 1.2571 4.5923E03

1.0793 1.4048E03 1.4568 1.5973E03

1.0233 5.6902E04 1.3038 6.1767E04

EOC Eu ðhÞ EOC

– 2.0066E02 –

1.2296 9.5224E03 1.0754

1.5236 4.6293E03 1.0405

1.3707 2.2784E03 1.0228

MassErr(T) EOC

2.1328E01 –

9.7067E02 1.1357

4.7704E02 1.0249

2.3450E02 1.0245

Table 3 The results obtained by the MCC scheme for Example 6.2.

subject to

1=h

10

20

40

80

yt þ r  ðv y  eryÞ ¼ f þ Bu in X  ð0; TÞ;

ð6:2Þ

Ey ðhÞ EOC Ey ðhÞ

4.0986E02 – 4.1842E02

2.0719E02 0.9842 2.1392E02

1.0512E02 0.9789 1.0965E02

5.1748E03 1.0225 5.4205E03

with 0 6 u 6 0:5 and B ¼ I. The co-state p satisfies the following equation:

EOC Ep ðhÞ EOC Ep ðhÞ

– 3.1028E01 – 3.1114E01

0.9679 1.5601E01 0.9919 1.5641E01

0.9642 7.7889E02 1.0021 7.8080E02

1.0164 3.8935E02 1.0004 3.9027E02

EOC Eu ðhÞ EOC

– 1.0505E01 –

0.9922 6.2221E02 0.7556

1.0023 3.5453E02 0.8115

1.0005 2.0552E02 0.7866

MassErr(T) EOC

1.9962E01 –

9.3116E02 1.1002

4.8531E02 0.9401

2.4180E02 1.0051

ðery  yv Þ  n ¼ g

on @ X  ð0; TÞ;

 pt  v  rp  erp ¼ y  zd in X  ð0; TÞ; erp  n ¼ 0 on @ X  ð0; TÞ:

ð6:3Þ

For simplicity we shall use the same mesh partition for T h and T hU , and it all runs Dt ¼ h in the numerical tests. To solve the optimal control problems, we use the C++ software package: AFEpack, it is available at http://dsec.pku.edu.cn/~rli/. For constrained optimal control problems governed by convection–diffusion equations, people pay more attention on the state and the control. Therefore in the following numerical examples, we mostly center on the state variable y, which is approximated by piecewise linear elements; and the control variable u, which is discretized using piecewise constant elements. For an error functional EðhÞ we define the experimental order of convergence by

EOC ¼

log Eðh1 Þ  log Eðh2 Þ : log h1  log h2

In Tables 1–4 we investigate the error functionals

Ev ðhÞ :¼ kv  v h kl1 ð0;T;L2 ðXÞÞ ; for



Ev ðhÞ :¼ kv  v h kle ð0;T;H1 ðXÞÞ ;

v ¼ y; p;

Table 4 The results obtained by the MMOC scheme for Example 6.2. 1=h

10

20

40

80

Ey ðhÞ EOC Ey ðhÞ

5.1330E02 – 5.1965E02

2.5395E02 1.0153 2.5901E02

1.2851E02 0.9827 1.3208E02

6.4217E03 1.0009 6.6184E03

EOC Ep ðhÞ EOC Ep ðhÞ

– 5.5014E01 – 5.5139E01

1.0045 2.6862E01 1.0342 2.6916E01

0.9716 1.3228E01 1.0220 1.3254E01

0.9969 6.5684E02 1.0100 6.5807E02

EOC Eu ðhÞ EOC

– 1.0606E01 –

1.0346 6.2818E02 0.7556

1.0220 3.5826E02 0.8102

1.0101 2.0731E02 0.7892

MassErr(T) EOC

3.5748E01 –

1.7867E01 1.0006

9.4263E02 0.9225

4.8032E02 0.9727

Eu ðhÞ :¼ ku  uh kl2 ð0;T;L2 ðXÞÞ :

Table 1 The results obtained by the MCC scheme for Example 6.2. 1=h

8

16

32

64

Ey ðhÞ EOC Ey ðhÞ

4.2045E02 – 7.1551E02

1.1757E02 1.8384 2.9723E02

3.1227E03 1.9127 1.4107E02

9.3915E04 1.7334 6.9743E03

EOC Ep ðhÞ EOC Ep ðhÞ

– 6.8163E03 – 8.2250E03

1.2674 2.9972E03 1.1854 3.6630E03

1.0752 9.6527E04 1.6346 1.1377E03

1.0163 2.9752E04 1.6979 3.4132E04

EOC Eu ðhÞ EOC

– 1.9959E02 –

1.1670 9.4888E03 1.0727

1.6869 4.6153E03 1.0398

1.7369 2.2712E03 1.0230

MassErr(T) EOC

2.0389E02 –

6.9099E03 1.5611

3.8229E03 0.8540

1.8303E03 1.0626

Besides, the errors in mass balance of the two schemes for the state are compared, where the errors are carried out by taking the absolute relative error for mass balance at time T, i.e.

MassErrðTÞ ¼

j

R X

R yNh T  X yðTÞj R : j X yðTÞj

Example 6.1 (A body rotation problem). This example is an adaption of Example 1 from [21] for the convection–diffusion equation. We take the spatial domain X ¼ ð1; 1Þ2 and the time interval ½0; T ¼ ½0; 0:5. The velocity field is imposed as v ¼ ð1 þ sinðt  x1 Þ; 1 þ sinðt  x2 ÞÞ and it is plotted in Fig. 1. In this example v depends both on spacial variable x and temporal variable t, and we can check that r  v – 0. The corresponding

90

H. Fu, H. Rui / Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 82–92

analytical solutions for problems (6.2) and (6.3) with a constant diffusion coefficient e ¼ 4:0  103 ; f ¼ u, and zd ¼ y are given by

1 0.8

    1  cosðt  x1 Þ 1  cosðt  x2 Þ exp  ; yðx; tÞ ¼ exp 

0.6

e

0.4

e

pðx; tÞ ¼ 0;

0.2

u0 ðx; tÞ ¼ sinðptÞ sinðpx1 =2Þ sinðpx2 =2Þ;

0

uðx; tÞ ¼ maxf0; minð0:5; u0  pÞg:

−0.2

The boundary function g can be determined by inserting the state function y into Eq. (6.2). This example is a moving hump problem but with its height of the hump center not changing in the course of the time for the state equation. Its initial hump center is fixed at ð0; 0Þ. Fig. 2 show the approximate state solution and its contour-lines for the MCC approximation at T ¼ 0:5, at that time the hump center moves to ð0:5; 0:5Þ. The elevation plot of the approximate control and its corresponding contour-lines at T ¼ 0:5 are presented in Fig. 3.

−0.4 −0.6 −0.8 −1 −1

−0.5

0

0.5

1

Fig. 1. The velocity field at T ¼ 0:5.

1

1

0.5

0

Y

SOLUTION

0.5

0 -1

-1

-0.5

-0.5

-0.5 Y

0

0 0.5

0.5

X

-1 -1

1

1

-0.5

0

X

0.5

1

0.5

1

Fig. 2. The approximate solution (left) and contour-lines (right) for the state.

SOLUTION

0.4 0.2 0

Y

0.5

0

-1

-1 -0.5

-0.5 Y

0

0

X

-0.5

0.5

0.5 1

1

-0.5

0

X Fig. 3. The approximate solution (left) and contour-lines (right) for the control.

91

H. Fu, H. Rui / Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 82–92

1

that the corresponding analytical solutions of problems (6.2)–(6.3) are

0.9

yðx; tÞ ¼ sinðpt=2Þ exp 

0.8 0.7

! ðx1  1=2Þ2 þ ðx2  1=2Þ2 pffiffiffi ;

e

wðx; tÞ ¼ expð2tÞ cosðpx1 Þ cosðpx2 Þ;

0.6

0.4

pðx; tÞ ¼ wðx; tÞ  wðx; TÞ; ( 1:0 x1 þ x2 > 1:0; zðx; tÞ ¼ 0:0 x1 þ x2 6 1:0;

0.3

u0 ðx; tÞ ¼ 1  cosðpx1 =2Þ  cosðpx2 =2Þ þ z;

0.2

uðx; tÞ ¼ maxf0; minð0:5; u0  pÞg:

0.1

The boundary function g is determined by inserting the state function y into Eq. (6.2). This example is a hump changing its height in the course of the time for the state, and the optimal control has a strong jump (discontinuity), introduced by u0 . Although the jump does not ’’move’’ with time, it may expand or contract with time because the zero-set of the co-state p varies with time. In Figs. 5,6, we show the approximate solutions and their corresponding contourlines for the state and control at T ¼ 1, respectively. Because of the affections of both the effective searching for the characteristic points x and x, and the numerical integration error, the numerical solutions may have a slightly perturbation. It can be seen from the contour-lines of the approximate state solution, some extra disgusting information occur. In this numerical test, we set h ¼ 1=N with N ¼ 10; 20; 40; 80. In Tables 3,4 numerical results are presented for the MCC and MMOC schemes, respectively. Notice that the MCC scheme gives better agreement to the analytic solutions than the MMOC. In addition, it still can be readily seen that the MCC is a conservative scheme where the errors in mass are smaller than those of MMOC. The results displayed in Tables 1, 3 and Figs. 2,3, 2,3 show that the MCC scheme used in this paper has a good approximation to the quadratic optimal control problem governed by singularly perturbed, convection-dominated diffusion equations.

0.5

0

0.2

0.4

0.6

0.8

1

Fig. 4. The velocity field for Example 6.2.

In this numerical test, we set h ¼ 2=N with N ¼ 16; 32; 64; 128. In Tables 1,2 numerical results and convergence orders are presented for the MCC and MMOC schemes, respectively. They both show a first-order convergence for the control and even better than first-order convergence for the state and co-state. However, it can be seen that the MCC method has simulated the optimal control problems better than the MMOC. In addition, it is clear that the relative mass errors obtained by the MCC are smaller than those obtained by the MMOC. Example 6.2 (A hump changing its height). For the second example, we consider an optimization problem with internal layer in X ¼ ð0; 1Þ2 and T ¼ 1. We take the velocity field as v ¼ ð1  2x2 Þ sinðpx1 Þ sinðpx2 Þ; ð2x1  1Þ sinðpx1 Þ sinðpx2 ÞÞ, and it is plotted in Fig. 4. The small diffusion constant e ¼ 104 . It is easy to see that in this example we have that r  v – 0. Let the desired state zd and the external source function f be chosen such

1

1

0.8

SOLUTION

0.8 0.6 0.4

0.6

Y

0

0.4 0.2 0 0

0

0.2

0.2

0.2 0.4 Y

0.4 0.6

0.6

0.8

0.8 1

X

0

0

0.2

0.4

1 Fig. 5. The approximate solution (left) and contour-lines (right) for the state.

X

0.6

0.8

1

92

H. Fu, H. Rui / Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 82–92

0.4

0.2

0.6

Y

SOLUTION

0.8

0.4 0 0.2

0 1

0.8

0.4 0.6

X

0.6 0.4

0.2

0.2

Y

0.8 0

1

0.2

0.4

X

0.6

0.8

1

Fig. 6. The approximate solution (left) and contour-lines (right) for the control.

7. Concluding remarks In this paper, we apply a class of MCC scheme to optimal control problems governed by singularly perturbed convection–diffusion equations subject to non-divergence-free velocity fields and pointwise inequality constraints on the control variable. We derive some a priori error estimates in e-weighted energy norm for the state and L2 -norm for the control. The proposed scheme is symmetric and stable. Numerical experiments are given to confirm the theoretical convergence and to show its effeciency. We shall consider a posteriori error estimates and adaptive computation in the future. Acknowledgements The authors thank the editor and the anonymous referee for their valuable comments and suggestions on an earlier version of this paper. The first author was supported by the National Natural Science Foundation of China (Nos. 11126086 and 11101431) and the Fundamental Research Funds for the Central Universities (No. 12CX04083A); The second author was supported by the National Natural Science Foundation of China (No. 11171190). References [1] J. Zhu, Q. Zeng, A mathematical formulation for optimal control of air pollution, Sci. China Ser. D 46 (2003) 994–1002. [2] D. Parra-Guevara, Y.N. Skiba, Elements of the mathematical modelling in the control of pollutants emissions, Ecol. Model. 167 (2003) 263–275. [3] T.J.R. Hughes, A. Rooks, Streamline upwind/Petrov Galerkin formulations for the convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 54 (1982) 199–259. [4] F. Brezzi, A. Russo, Choosing bubbles for advection–diffusion problems, Math. Models Methods Appl. Sci. 4 (1994) 571–587. [5] C. Johnson, J. Pitkränta, An analysis of the discontinuous Galerkin method for scalar hyperbolic equation, Math. Comput. 46 (1986) 1–26. [6] E. Burman, P. Hansbo, Edge stabilization for Galerkin approximations of convection–diffusion–reaction problems, Comput. Methods Appl. Mech. Eng. 193 (2004) 1437–1453.

[7] J. Douglas, T.F. Russell, Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures, SIAM J. Numer. Anal. 19 (1982) 871– 885. [8] P. Houston, E. Süli, Adaptive Lagrange–Galerkin methods for unsteady convection–diffusion problems, Math. Comput. 70 (2000) 77–106. [9] R.E. Ewing, T.F. Russell, M.F. Wheeler, Convergence analysis of an approximation of miscible displacement in displacement in porous media by mixed finite elements and a modified method of characteristics, Comput. Methods Appl. Mech. Eng. 47 (1984) 73–92. [10] H. Rui, M. Tabata, A second order characteristic finite element scheme for convection–diffusion problems, Numer. Math. 92 (2002) 161–177. [11] J. Douglas, T.C. Huang, F. Pereira, The modified method of characteristics with adjusted advection, Numer. Math. 83 (1999) 353–369. [12] M.A. Celia, T.F. Russell, I. Herrera, R.E. Ewing, An Eulerian–Lagrangian localized adjoint method for the advection–diffusion equation, Adv. Water Res. 13 (1990) 187–206. [13] T. Arbogast, M.F. Wheeler, A characteristics-mixed finite element method for advection-dominated transport problems, SIAM J. Numer. Anal. 32 (1995) 404–424. [14] S.S. Collis, M. Heinkenschloss, Analysis of the streamline upwind/Petrov Galerkin method applied to the solution of optimal control problems, CAAM TR02-01, 2002. [15] M. Heinkenschloss, D. Leykekhman, Local error estimates for SUPG solutions of advection-dominated elliptic linear-quadratic optimal control problems, SIAM J. Numer. Anal. 47 (2010) 4607–4638. [16] R. Becker, B. Vexler, Optimal control of the convection–diffusion equation using stabilized finite element methods, Numer. Math. 106 (2007) 349– 367. [17] M. Braack, Optimal control in fluid mechanics by finite elements with symmetric stabilization, SIAM J. Control Optim. 48 (2009) 672–687. [18] N. Yan, Z. Zhou, A priori and a posteriori error analysis of edge stabilization Galerkin method for the optimal control problem governed by convectiondominated diffusion equation, J. Comput. Appl. Math. 223 (2009) 198– 217. [19] N. Yan, Z. Zhou, A RT mixed FEM/DG scheme for optimal control governed by convection diffusion equations, J. Sci. Comput. 41 (2009) 273–299. [20] H. Fu, H. Rui, A priori error estimates for optimal control problems governed by transient advection–diffusion equations, J. Sci. Comput. 38 (2009) 290– 315. [21] H. Rui, M. Tabata, A mass-conservative characteristic finite element scheme for convection–diffusion problems, J. Sci. Comput. 43 (2010) 416– 432. [22] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2002. [23] J.L. Lions, Optimal Control of Systems Governed by Partial Differential Equations, Springer-Verlag, Berlin, 1971.