On gradient-type optimization method utilizing mixed finite element approximation for optimal boundary control problem governed by bi-harmonic equation

Applied Mathematics and Computation 186 (2007) 1429–1440 www.elsevier.com/locate/amc

On gradient-type optimization method utilizing mixed finite element approximation for optimal boundary control problem governed by bi-harmonic equation Li Bingjie

a,b,*

, Liu Sanyang

a

a

b

Department of Mathematics, Xidian University, Xi’an, Shaan xi 710071, China College of Science, Air Force Engineering University, Xi’an, Shaan xi 710051, China

Abstract In this paper, a numerical method based on mixed finite element method for optimal boundary control problem governed by bi-harmonic equation is presented. The system of optimality equations consisting of state and costate function is derived. Based on the system, a gradient-type optimization method using a mixed finite approximation for solving the optimal boundary control is developed. In our algorithm, the trace of vortex function of costate equation on boundary plays a key role. Finally, a local error analysis at every iteration for this gradient-type optimization method is given.  2006 Elsevier Inc. All rights reserved. Keywords: Optimal control; Mixed finite element; Bi-harmonic equation; Error estimate

1. Introduction In this paper, we are interested in the optimal boundary control problem that the state equation is governed by bi-harmonic equation 1 m 2 2 min J ðu; wÞ ¼ kwðyÞ  wz ðyÞkL2 ðXÞ þ kuðyÞkL2 ðCÞ ; ð1:1Þ 2 2 ð1:2Þ D2 wðyÞ ¼ f ðyÞ; y 2 X  R2 ; ow ðyÞ ¼ uðyÞ; y 2 C; ð1:3Þ wðyÞ ¼ 0; on where D is Laplace operator, X a bounded polygonal domain in R2 and C its boundary. wz(y) 2 L2(X), f(y) 2 L2(X) are given functions, and m is the weight of the cost of the control function u(y). The normal vector n always points into the exterior of X. We assume for the sake of simplicity that X is convex and m is not zero. This optimal control system can be used to study the control or reduction techniques of noise emanating from *

Corresponding author. Address: Department of Mathematics, Xidian University, Xi’an, Shaan xi 710071, China. E-mail addresses: [email protected] (L. Bingjie), [email protected] (L. Sanyang).

0096-3003/$ - see front matter  2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.07.142

1430

L. Bingjie, L. Sanyang / Applied Mathematics and Computation 186 (2007) 1429–1440

vibrating structures and some control problems of elastic plate bending vibration [7,10,12,13,15]. The behavior of bending vibration can be controlled by the conditions imposed on the boundary. Here the function wz(y) is the desired target state of w(y), then the problem of finding the boundary control u to achieve the desired target state is cast into the least square sense by (1.1). Generally speaking, the habit for solving the numerical solution is standard finite element method (FEM) [11,16]. However, standard finite element computation require rather sophisticated finite element such as the 21-degree-of-freedoms or non-conforming element of Hermite type. Thus it is very expensive compute and further leads to some inconveniences for numerical algorithm. But we found that the mixed finite element method is available for solving this kind of problem. Some previous works have been made about mixed finite element method (for example, Ciarlet–Raviart, Herrmann–Miyoshi, Hellan–Herrmann–Johnson mixed methods) [1,6,7,9,10,14,17,18]. Although these methods are widely used in numerical computations for partial differential equations, it has not been fully utilized in optimal control. Especially, much less research is available for the optimal boundary control problems governed by bi-harmonic equation because the boundary control function may bring some of theoretical and numerical difficulties. Motivated by the above works, in this paper, we develop the mixed finite element method for solving the problem (1.1)–(1.3). In Section 2, we introduce the costate equation with clamped boundary conditions where the costate function is k(y), and develop the system of optimality equations consisting of state and costate function. From the optimality system we discover that the vortex function Dk is just equal to mu on C in the optimal sense (see the Section 2: (2.8)), which automatically generate a suitable gradient. In Section 3, we derive a gradient-type optimization scheme, and propose a mixed finite element method for the system of optimality equations. We point out that the mixed finite element approaches employed in our work can be generalized to other optimization algorithms, for example, conjugate gradient method [16]. Section 4 focuses the local truncated error estimate at every iteration. 2. Optimal control problem In this section, we describe the optimality system of the optimal boundary control problem (1.1)–(1.3). 1

Theorem 2.1. Assume that X is convex, then there exists a unique solution ðw ; u Þ 2 ðH 2 ðXÞ \ H 10 ðXÞÞ  H 2 ðCÞ to the optimal control problem (1.1)–(1.3) such that 8 2 D wðyÞ ¼ f ðyÞ; y 2 X  R2 ; > > > > ow > y 2 C; > < wðyÞ ¼ 0; on ðyÞ ¼ uðyÞ; 2 ð2:1Þ D kðyÞ ¼ ðwðyÞ  wz ðyÞÞ; y 2 X  R2 ; > > > ok > y 2 C; > kðyÞ ¼ 0; on ðyÞ ¼ 0; > : muðyÞ ¼ DkðyÞ; y 2 C: 1

Proof. For any u 2 H 2 ðCÞ, there exists the solution w 2 H 2 ðXÞ \ H 10 ðXÞ of the problem (1.2), (1.3). Let w = w(u), where w(u) denotes the solution of (1.2), (1.3) as a function of u, then the mapping u ! w(u) form 1 1 H 2 ðCÞ to H 2 ðXÞ \ H 10 ðXÞ is affine and continuous [2,3,8]. Let w 0 (u, du) be its first derivative at u 2 H 2 ðCÞ for w(u) in the direction du, where du is the variation for u. We have ( 2 y 2 X  R2 ; D w0 ðu; duÞ ¼ 0; ð2:2Þ 0 ¼ du; y 2 C: w0 ðu; duÞ ¼ 0; ow ðu;duÞ on Let J^ðuÞ ¼ J ðwðuÞ; uÞ, then the mapping u ! J^ðuÞ is twice Frechet differentiable, and its first derivative at u is given by Z Z J^0 ðu; duÞ ¼ ðw  wz Þw0 ðu; duÞ dX þ m udu dC; X

C

and its second derivative at u is given by Z 2 2 J^00 ðuÞðdu; duÞ ¼ kw0 ðu; duÞkL2 ðXÞ þ ðw  wz Þw00 ðuÞðdu; duÞ dX þ mkdukL2 ðCÞ : X

L. Bingjie, L. Sanyang / Applied Mathematics and Computation 186 (2007) 1429–1440

1431

Because the second derivative of w(u) is zero, obviously, the functional J^ is uniformly convex. Since the domain X is convex in R2, this implies existence of a unique solution u* to (1.1)–(1.3). Suppose that (w*, u*) is a solution of (1.1)–(1.3), w* = w(u*), the solution is characterized by the first order optimality condition J^0 ðu ; duÞ ¼ 0 for all du J^0 ðu ; duÞ ¼ ðw  wz ; w0 ðu ; duÞÞL2 ðXÞ þ mðu ; duÞL2 ðCÞ ¼ 0:

ð2:3Þ

Introducing a costate function k* as the unique solution to the problem (with clamped boundary conditions) ( D2 k ðyÞ ¼ ðw ðyÞ  wz ðyÞÞ; y 2 X  R2 ; ð2:4Þ  k ðyÞ ¼ 0; ok ðyÞ ¼ 0; y 2 C: on Then by the boundary value problems (2.2) and (2.4), we have Z Z 2  0  w ðu ; duÞD k dX ¼  w0 ðu ; duÞðw  wz Þ dX; X

Z

ð2:5Þ

X

k D2 w0 ðu ; duÞ dX ¼ 0:

ð2:6Þ

X

Utilizing the Green formula of bi-harmonic operator, we have  Z Z  Z 0  oDk   ow ðu ; duÞ 0  0   w ðu ; duÞ Dw ðu ; duÞDk dX  Dk dC ¼  w0 ðu ; duÞðw  wz Þ dX; on on X C X  Z Z   0  ok oDw ðu ; duÞ   k Dk Dw0 ðu ; duÞ dX  Dw0 ðu ; duÞ dC ¼ 0: on on X C By and the boundary conditions of (2.2) and (2.4), we have Z Z Dk du dC ¼ w0 ðu ; duÞðw  wz Þ dX: C

ð2:7Þ

X

By the above equality and (2.3), there have J^0 ðu ; duÞ ¼ ðDk ; duÞL2 ðCÞ þ mðu ; duÞL2 ðCÞ ¼ 0;

1

8du 2 H 2 ðCÞ:

ð2:8Þ

Obviously, by the fundamental theorem of variational, the last equality in (2.1) holds, and ðw ; u Þ 2 1 ðH 2 ðXÞ \ H 10 ðXÞÞ  H 2 ðCÞ. This completes the proof. h Remark 2.1. Theorem 2.1 implies that the solution to the optimality system (2.1) is equivalent to the solution of (1.1)–(1.3). It is well known that if the largest interior angle of the boundary C is less than 126.28 and w  wz 2 L2(X), then k 2 H4(X). If X is a convex polygonal domain then k 2 H3+s(X), where 0 < s 6 1 is a constant that depends on the largest interior angle of the boundary [1]. From above, we can see that for each w  wz 2 L2(X), the costate problem (2.4) with homogeneous boundary conditions admits only one solution satisfying k 2 H 3 ðXÞ \ H 20 ðXÞ [7,18]. Thus, from (2.1), we can conclude the regularity property [2,3] that if 1 wz 2 L2(X), f 2 L2(X), then ðw ; k ; u Þ 2 ðH 2 ðXÞ \ H 10 ðXÞÞ  ðH 3 ðXÞ \ H 20 ðXÞÞ  H 2 ðCÞ. 3. Mixed finite element method 3.1. Gradient-type algorithm in the continuous case The mixed finite element method is used to simultaneously approximate the flow function and the vortex function. Let by w and u denote the flow function and the vortex function of the state equation, respectively, where u = Dw. Let by k and w denote the flow function and the vortex function of the costate function, respectively, where w = Dk.

1432

L. Bingjie, L. Sanyang / Applied Mathematics and Computation 186 (2007) 1429–1440

Let c0, c1 be the following trace mapping:  ow c0 w ¼ wjC ; c1 w ¼  ; on C

3

1

then the mapping {c0, c1} is linear continuous and onto H 2 ðCÞ  H 2 ðCÞ from H2(X). The problem (1.2), (1.3) is equivalent to 8 DuðyÞ ¼ f ðyÞ; y 2 X  R2 ; > > < DwðyÞ ¼ uðyÞ; y 2 X  R2 ; > > : c0 wðyÞ ¼ 0; c1 wðyÞ ¼ uðyÞ; y 2 C:

ð3:1Þ

It is well known that there exists a linear operator S : ujC !  ow jC , where S is an isomorphism from 1 1 on H ðCÞ onto H 2 ðCÞ and that it induces asymmetric, continuous, H 2 ðCÞ-elliptic bilinear form via Z ow 1 ~Þ ¼ ðSujC ; r ~Þ ¼  ~ dC; 8~ SðujC ; r r r 2 H 2 ðCÞ: on C 12

1

~, then for any u 2 H 2 ðCÞ, f 2 L2(X), by the mixed finite Moreover, if we denote by r the extension to X of r element method [5] and the Green formula, we can get the following variational problem corresponding to (1.2), (1.3): Find ðu; wÞ 2 H 1 ðXÞ  H 10 ðXÞ such that (R R R ur dX  X rr  rw dX ¼  C ur dC; 8r 2 H 1 ðXÞ; X R R ð3:2Þ ru  rl dX ¼ X f l dX; 8l 2 H 10 ðXÞ: X Similarly, we consider the costate problem (2.4) with homogeneous boundary conditions. Because of 1 k 2 H 3 ðXÞ \ H 20 ðXÞ, we can see that w = Dk 2 H1 (X), Dw = w  wz 2 L2(X), and c0 w 2 H 2 ðCÞ, namely 1 1 c0 Dk 2 H 2 ðCÞ. This result is essential owing to the control u 2 H 2 ðCÞ. In the following we can see that the trace of w (or Dk) on C will play a key role, both theoretically and numerically. By the mixed finite element method, we can obtain the variational problem corresponding to the costate problem (2.4): Find ðw; kÞ 2 H 1 ðXÞ  H 10 ðXÞ such that (R R wr dX  X rr  rk dX ¼ 0; 8r 2 H 1 ðXÞ; RX R ð3:3Þ rw  rl dX ¼  X ðw  wz Þl dX; 8l 2 H 10 ðXÞ: X In the following, we firstly propose a gradient-type algorithm for the optimality system (2.1) in the continuous case. The optimal control u* is determined by the optimality condition J~0 ðu; duÞ ¼ 0:

ð3:4Þ ð0Þ

1 2

An arbitrary function u 2 H ðCÞ may be specified as an initial guess for the control input (in a general way, set u(0) = 0). Firstly, we set i = 0, and solve the problems, respectively ( 2 ð0Þ D w ðyÞ ¼ f ðyÞ; y 2 X  R2 ; ð3:5Þ ð0Þ wð0Þ ðyÞ ¼ 0; owon ðyÞ ¼ uð0Þ ; y 2 C; ( 2 ð0Þ D k ðyÞ ¼ ðwð0Þ ðyÞ  wz ðyÞÞ; y 2 X  R2 ; ð3:6Þ ð0Þ y 2 C; kð0Þ ðyÞ ¼ 0; okon ðyÞ ¼ 0; to get the gradient g(0) = mu(0) + c0Dk(0) at u(0). Then for i P 0, assuming that the iterative control u(i), and the gradient g(i) are known, we update u(i) via u(i+1) = u(i)  qig(i), where qi is the optimum iterative step at u(i) such that min J ðwðuðiÞ  qgðiÞ Þ; uðiÞ  qgðiÞ Þ ¼ J ðwðuðiÞ  qi gðiÞ Þ; uðiÞ  qi gðiÞ Þ: qP0

L. Bingjie, L. Sanyang / Applied Mathematics and Computation 186 (2007) 1429–1440

1433

Obviously, qi is characterized as ðmðuðiÞ  qi gðiÞ Þ þ c0 DkðiÞ ; gðiÞ ÞL2 ðCÞ ¼ 0, where kðiÞ is the solution of ( 2  ðiÞ ðyÞ ¼ f ðyÞ; Dw y 2 X  R2 ; ðiÞ

 ðiÞ ðyÞ ¼ 0; owon ðyÞ ¼ uðiÞ  qi gðiÞ ; y 2 C; w ( 2 D kðiÞ ðyÞ ¼ ð wðiÞ ðyÞ  wz ðyÞÞ; y 2 X  R2 ;  kðiÞ ðyÞ ¼ 0;

ðiÞ ok on

ðyÞ ¼ 0;

y 2 C:

~ ðiÞ and ~ kðiÞ are respectively the solutions of Assume that w ( 2 ðiÞ ~ ðyÞ ¼ 0; y 2 X  R2 ; Dw ~ ðiÞ ðyÞ ¼ 0; w (

o~ wðiÞ on

ðyÞ ¼ gðiÞ ;

kðiÞ ðyÞ ¼ ~ D2 ~ wðiÞ ðyÞ; ~ kðiÞ ðyÞ ¼ 0;

~ðiÞ ok on

ðyÞ ¼ 0;

ð3:7Þ

y 2 C;

y 2 X  R2 ;

ð3:8Þ

y 2 C;

w(i) and k(i) are respectively the solutions of ( 2 D wðiÞ ðyÞ ¼ f ðyÞ; y 2 X  R2 ; (

wðiÞ ðyÞ ¼ 0;

owðiÞ on

ðyÞ ¼ uðiÞ ; y 2 C;

D2 kðiÞ ðyÞ ¼ ðwðiÞ ðyÞ  wz ðyÞÞ; y 2 X  R2 ; kðiÞ ðyÞ ¼ 0;

okðiÞ on

ðyÞ ¼ 0;

y 2 C;

kðiÞ holds owing to the linear characteristics of the considered boundary value problems, then  kðiÞ ¼ kðiÞ  qi ~ and the the optimum iterative step qi is characterized as ðmðuðiÞ  qi gðiÞ Þ þ c0 DkðiÞ  qi c0 D~ kðiÞ ; gðiÞ ÞL2 ðCÞ ¼ 0: Obviously, g(i) = mu(i) + c0Dk(i) holds, Let g~ðiÞ ¼ mgðiÞ þ c0 D~kðiÞ , then we have qi ¼

ðgðiÞ ; gðiÞ ÞL2 ðCÞ ð~ gðiÞ ; gðiÞ ÞL2 ðCÞ

:

ð3:9Þ

It is easy to prove that gðiþ1Þ ¼ gðiÞ  qi g~ðiÞ holds. From above, the gradient-type algorithm can be given as follows. Algorithm 1. Gradient-type algorithm in the continuous case Step 1. Choose an initial boundary control u(0), and set e P 0. Step 2. Solve the problems (3.5) and (3.6), set g(0) = mu(0) + c0Dk(0). Set i = 0. kgðiÞ k2 2 Step 3. Solve the problems (3.7) and (3.8), and set g~ðiÞ ¼ mgðiÞ þ c0 D~kðiÞ , qi ¼ ð~gðiÞ ;gðiÞLÞ ðCÞ , u(i+1) = u(i)  qig(i), L2 ðCÞ gðiþ1Þ ¼ gðiÞ  qi g~ðiÞ . Step 4. If

kgðiþ1Þ k2 2

L ðCÞ

kgð0Þ k2 2

6 e, then take u* = u(i+1); else, go to Step 5.

L ðCÞ

Step 5. Set i:¼i + 1, and go to Step 3. 3.2. Mixed finite element algorithm based on gradient-type method In the following we derive the mixed finite element algorithm based on gradient-type method. Firstly, we consider the variational problem corresponding to (3.5), (3.6): Find ðuð0Þ ; wð0Þ Þ 2 H 1 ðXÞ  H 10 ðXÞ such that (R R R uð0Þ r dX  X rr  rwð0Þ dX ¼  C uð0Þ r dC; 8r 2 H 1 ðXÞ; X R R ð3:10Þ ruð0Þ  rl dX ¼ X f l dX; 8l 2 H 10 ðXÞ; X

1434

L. Bingjie, L. Sanyang / Applied Mathematics and Computation 186 (2007) 1429–1440

and find ðwð0Þ ; kð0Þ Þ 2 H 1 ðXÞ  H 10 ðXÞ such that ( R ð0Þ R w r dX  X rr  rkð0Þ dX ¼ 0; X R R rwð0Þ  rl dX ¼  X ðwð0Þ  wz Þl dX; X

8r 2 H 1 ðXÞ; 8l 2 H 10 ðXÞ:

ð3:11Þ

~ ðiÞ Þ 2 H 1 ðXÞ  H 10 ðXÞ Then we consider the variational problems corresponding to (3.7), (3.8): Find ð~ uðiÞ ; w such that (R R R ~ ðiÞ r dX  X rr  r~ wðiÞ dX ¼  C gðiÞ r dC; 8r 2 H 1 ðXÞ; u X R ð3:12Þ r~ uðiÞ  rl dX ¼ 0; 8l 2 H 10 ðXÞ; X ~ ðiÞ ; ~ and find ðw kðiÞ Þ 2 H 1 ðXÞ  H 10 ðXÞ such that (R R ~ ðiÞ r dX  rr  r~ w kðiÞ dX ¼ 0; 8r 2 H 1 ðXÞ; X X R R ðiÞ ~ ðiÞ  rl dX ¼  w ~ l dX; rw 8l 2 H 10 ðXÞ: X X

ð3:13Þ

The classical theory of mixed finite element method guarantees the existence of a unique of the above problems, respectively [4,5,19]. But for the convenience of the readers, we sketch a depiction. We consider an abstract framework which is applicable to the problems above. Let X and M be two Hilbert spaces with inner products denoted by (Æ, Æ)X and (Æ, Æ)M, norms denoted by kÆkX and kÆkM, and dual spaces denoted by X 0 and M 0 , respectively. Let a(Æ, Æ) and b(Æ, Æ) be bounded bilinear forms on X · X, X · M, respectively, and g, f be two bounded linear functionals in X 0 and M 0 respectively, with hÆ, Æi denoting the duality. Find (u, w) 2 X · M such that aðu; rÞ þ bðr; wÞ ¼ hg; ri; 8r 2 X ; bðu; lÞ ¼ hf ; li;

8l 2 M:

ð3:14Þ

Let B denote the linear operator defined by hBr; liM 0 M ¼ bðr; lÞ; 8r 2 X ;

l 2 M:

The kernel of B: Ker B = {u0 2 Xj b(u0, l) = 0, "l 2 M}. We recall that a sufficient condition for problem (3.14) to be well-posed is the coercivity of the bilinear form a(Æ, Æ) on the kernel of B, and the surjectivity of the linear operator B : X ! M, on the other hand, the following conditions hold 2

9a0 > 0; s:t: aðu0 ; u0 Þ P a0 ku0 kX ; 8u0 2 Ker B; bðr; lÞ P b0 klkM ; 8l 2 M: 9b0 > 0; s:t: sup r6¼0 krkX

ð3:15Þ ð3:16Þ

A define element approximation of (3.14) consisting in choosing define dimensional spaces Xh  X, Mh  M and in solving the discrete problem: Find (uh, wh) 2 Xh · Mh such that aðuh ; rh Þ þ bðrh ; wh Þ ¼ hg; rh i; 8rh 2 X h ; bðuh ; lh Þ ¼ hf ; lh i;

8lh 2 M h :

ð3:17Þ

It is standard to obtain that (3.17) is unique solvable if there exist two positive constants ah,0 and bh,0 such that 2

aðuh;0 ; uh;0 Þ P ah;0 kuh;0 kX ; bðrh ; lh Þ P bh;0 klh kM ; sup rh 6¼0 krh kX

8uh;0 2 Ker Bh ; 8lh 2 M h ;

ð3:18Þ ð3:19Þ

where Ker Bh = {uh,0 2 Xhjb(uh,0, lh) = 0, "lh 2 Mh}. Moreover if both ah,0 and bh,0 are uniformly bounded away from zero with respect to h, the following error estimate holds: ku  uh kX þ kw  wh kM 6 C inf fku  rh kX þ kw  lh kM g: rh ;lh

ð3:20Þ

L. Bingjie, L. Sanyang / Applied Mathematics and Computation 186 (2007) 1429–1440

1435

In general, finding finite element spaces satisfying simultaneously both (3.18) and (3.19) is not a trivial task, because the kernel of Bh is not a subspace of kernel of B, so that (3.18) is not a consequence of (3.15) (unless the bilinear form a(Æ, Æ) is coercive on the whole space X). In many instance, condition (3.18) and (3.19) impose contradictory requirement on the choice of the discrete spaces Xh and Mh, and only quite special choices are available. Thus the appropriate choices of Xh and Mh are essential. Indeed condition (3.19) essentially requires that the space Xh is sufficiently ‘‘rich’’ with respect to Mh. If Xh is too large, the Ker Bh may consequently grow, and condition (3.18) might fail [4]. 1 1 In the problems (3.10)–(3.13), obviously, because R R of X ¼ H ðXÞ; R M ¼ H 0 ðXÞ, the bilinear forms are characterized by aðu; rÞ ¼ X ur dX and bðr; lÞ ¼  X rr  rl dX ¼ X lDr dX ¼ hDr; liM 0 M , and the continuous linear operator B is Br = Dr, and the kernel of B is Ker B = {u0 2 XjDu0 = 0}. Some previous works about the appropriate choices of Xh and Mh are available [5,17,18]. In the following we are to use some results of these literatures. We consider the discretizations of the problems (3.10)–(3.13). Let sh = {K} be a quasi-uniform partition of X with h the maximum diameter of the partition such that: sh finite, K  X, "K 2 sh, [K2sh K ¼ X, and for any K, K 0 2 sh, K 5 K 0 , at most, K and K 0 have a side or a vertex in common, Let Pm be the space of polynomials of two variables of degree less than or equal to m. We know that if the flow function and vortex function are approximated by piecewise quadratic polynomials the mixed finite element scheme can achieve the approximation order for associated error estimates of the bi-harmonic equation [18] and thus we usually introduce the following finite dimensional spaces X h :¼ fu 2 C 0 ðXÞ; ujK 2 P m ; m P 2; 8K 2 sh g; M h :¼ X h \ H 10 ðXÞ: ð0Þ

ð0Þ

Consider the discrete variational problems used to approximate (3.10)–(3.13): Find ðuh ; wh Þ 2 X h  M h , such that ( R ð0Þ R R ð0Þ ð0Þ u rh dX  X rrh  rwh dX ¼  C uh rh dC; 8rh 2 X h ; X h ð3:21Þ R R ð0Þ ruh  rlh dX ¼ X f lh dX; 8lh 2 M h ; X ð0Þ

ð0Þ

and find ðwh ; kh Þ 2 X h  M h , such that ( R ð0Þ R ð0Þ w rh dX  X rrh  rkh dX ¼ 0; 8rh 2 X h ; X h R R ð0Þ ð0Þ rwh  rlh dX ¼  X ðwh  wz Þlh dX; 8lh 2 M h : X

ð3:22Þ ðiÞ

ðiÞ

~ h Þ 2 X h  M h , such Then we consider the variational problems corresponding to (3.7), (3.8): Find ð~ uh ; w that ( R ðiÞ R R ðiÞ ðiÞ ~ r dX  X rrh  r~ wh dX ¼  C gh rh dC; 8rh 2 X h ; u X h h ð3:23Þ R ðiÞ r~ uh  rlh dX ¼ 0; 8lh 2 M h ; X ðiÞ ~ ðiÞ ; ~ and find ðw h kh Þ 2 X h  M h , such that ( R ðiÞ R ðiÞ ~ r dX  rrh  r~ w kh dX ¼ 0; 8rh 2 X h ; X h h X R ðiÞ R ~ ðiÞ  rl dX ¼  w ~ l dX; rw 8lh 2 M h : h h X X h h

ð3:24Þ

Now let us to consider the discrete variational problem (3.21). We explicitly introduce the finite element basis set N þN B

I X h ¼ spanfaj gj¼1

;

N

I M h ¼ spanfaj gj¼1 ;

so that we have NI interior and NB boundary degrees of freedom respectively. Let X Ih ; X Bh  X h satisfy X Ih  ð0Þ ð0Þ X Bh X h where X Ih ¼ fu 2 X h ; ujC ¼ 0g, then uh and wh can be written ð0Þ

uh ¼

NI X j¼1

ð0Þ

uh;j aj þ

NX I þN B j¼N I þ1

ð0Þ

uh;j aj ;

ð0Þ

wh ¼

NI X j¼1

ð0Þ

wh;j aj :

1436

L. Bingjie, L. Sanyang / Applied Mathematics and Computation 186 (2007) 1429–1440

Of course, we have to choose basis set with local support for approximating boundary functions. Let bj = ajjC, NI + 1 6 j 6 NI + NB, denote the piecewise polynomial functions obtained by restricting the basis functions aj to the boundary C. We note here that fbj gN I þ16j6N I þN B span a subspace of dimension NB. Because ð0Þ the basis functions {bj} have a small support, the boundary function uh can be written ð0Þ

uh ¼

NB X

ð0Þ

uh;j bN I þj :

j¼1

Consider the linear system (3.21) in block matrix forms with the notations: Ah ¼ ðakj ÞðN I þN B ÞðN I þN B Þ , akj = a(ak, aj), Bh ¼ ðbkj ÞN I ðN I þN B Þ , bkj = b(aj, ak), T

ðbk ; bj ÞL2 ðCÞ , f ¼ ðf1 ; f2 ; . . . ; fN I Þ , fj ¼ ðf ; aj ÞL2 ðXÞ , ð0Þ ð0Þ ð0Þ ð0Þ ð0Þ ð0Þ ð0Þ ðuh;1 ; uh;2 ; . . . ; uh;N I þN B ÞT , wh ¼ ðwh;1 ; wh;2 ; . . . ; wh;N I ÞT , the finite element version of (3.21) 0 1 ! ! 0 ð0Þ T Uh Ah Bh B ð0Þ C ¼ @ Uh uh A: ð0Þ Bh 0 wh f

Uh ¼ ðukj ÞN B N B , ð0Þ ð0Þ ð0Þ T j = 1, 2, . . . , NI, ¼ ðuh;1 ; uh;2 ; . . . ; uh;N B Þ , and denote the transpose of Bh by BTh . We can ð0Þ uh

ð0Þ

ð0Þ

ð0Þ

ukj ¼ ð0Þ

Uh ¼ obtain

ð0Þ T

Obviously, we can decompose the matrix Bh, and the vector Uh : Bh ¼ ðBI BB Þ, Uh ¼ ðUI ; UB Þ . The matrix Ah is replaced by the diagonal matrix diag(AI, AB). Then the system (3.21) can be written in the form of a block system 0 10 ð0Þ 1 0 1 UI 0 AI 0 BTI B CB ð0Þ C ð0Þ C C B ð3:25Þ @ 0 AB BTB AB @ UB A ¼ @ Uh uh A: ð0Þ f BI BB 0 w h

In a similar way, the linear system (3.22) in block matrix forms is 0 10 ð0Þ 1 0 1 WI 0 AI 0 BTI B CB ð0Þ C C C B 0 @ 0 AB BTB AB A; @ WB A ¼ @ ð0Þ ð0Þ ðw  wz Þ BI BB 0 Kh ð0Þ

ð0Þ

ð0Þ T

ð0Þ

ð3:26Þ

ð0Þ

T

where wð0Þ ¼ ðw1 ; w2 ; . . . ; wN I Þ , wj ¼ ðwh ; aj ÞL2 ðXÞ , j = 1, 2, . . . , NI, wz ¼ ðwz;1 ; wz;2 ; . . . ; wz;N I Þ , wz;j ¼ ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ð0Þ

ðwz ; aj ÞL2 ðXÞ , Kh ¼ ðkh;1 ; kh;2 ; . . . ; kh;N I ÞT , WI ¼ ðwh;1 ; wh;2 ; . . . ; wh;N I ÞT , WB ¼ ðwh;N I þ1 ; wh;N I þ2 ; . . . ; wh;N I þN B ÞT . PN B ð0Þ ð0Þ ð0Þ wh;N I þj bN I þj . The trace of wh on C can be express as c0 wh ¼ j¼1 In a similar way, the linear systems (3.23) and (3.24) in block matrix forms are, respectively 0 10 ðiÞ 1 0 1 ~I U 0 AI 0 BTI B C B C C B ðiÞ C B 0 AB BT CB ~ ðiÞ Uh gh C ð3:27Þ C¼B U B AB B @ @ A; @ A 0 BI BB 0 ~ ðiÞ w h 0 1 0 1 0 1 ~ ðiÞ W 0 AI 0 BTI I B C C B C C B B 0 AB BT CB ~ ðiÞ 0 C ð3:28Þ C¼B W B AB B @ A; @ @ A ~ ðiÞ w ~ ðiÞ BI BB 0 K h

ðiÞ ðiÞ ðiÞ ðiÞ ðiÞ ðiÞ ðiÞ ðiÞ T ~ ðiÞ ~ ðiÞ ~ ðiÞ ~ ðiÞ ~ h;2 ; . . . ; w ~ h;N I ÞT , where gh ¼ ðgh;1 ; gh;2 ; . . . ; gh;N B ÞT , w wh;1 ; w U uh;1 ; u I ¼ ð~ h ¼ ð~ h;N I Þ , h;2 ; . . . ; u ðiÞ T T ~ ðiÞ ~ ðiÞ ~ ðiÞ T ~ ðiÞ ~ ðiÞ ~ ðiÞ ~ ðiÞ ~ ðiÞ ~ ðiÞ ~ ðiÞ ð~ uh;N I þ1 ; u h;N I þN B Þ , and WI ¼ ðwh;1 ; wh;2 ; . . . ; wh;N I Þ , WB ¼ ðwh;N I þ1 ; wh;N I þ2 ; . . . ; wh;N I þN B Þ , h;N I þ2 ; . . . ; u ðiÞ ðiÞ ðiÞ ðiÞ ðiÞ ðiÞ ðiÞ ðiÞ T ~ ðiÞ on ~ ðiÞ ¼ ð~ ~2 ; . . . ; w ~ N I ÞT , w ~ j ¼ ð~ ð~ kh;1 ; ~ kh;2 ; . . . ; ~ kh;N I Þ , w w1 ; w wh ; aj ÞL2 ðXÞ , j = 1, 2, . . . , NI. The trace of w h P ~ ðiÞ ¼ N B w ~ ðiÞ b be express as c0 w . h j¼1 h;N I þj N I þj

~ ðiÞ U B ¼ ðiÞ ~ Kh ¼ C can

L. Bingjie, L. Sanyang / Applied Mathematics and Computation 186 (2007) 1429–1440

1437

All the foregoing works lead us to the following algorithm. Algorithm 2. Mixed finite element algorithm based on gradient type method ð0Þ

ð0Þ

ð0Þ

ð0Þ T

Step 1. Choose an initial control vector uh ¼ ðu1 ; u2 ; . . . ; uN B Þ , and set e P 0. ð0Þ Step 2. Solve the system of linear Eqs. (3.25) to determine wh , and solve the system of linear Eqs. (3.26) to ð0Þ ð0Þ ð0Þ ð0Þ ð0Þ determine Kh and Wh . Set gh ¼ muh  WB . Set i = 0. ~ ðiÞ Step 3. Solve the system of linear (3.27) to determine w h , and solve the system of linear (3.28) to determine ðiÞ ðiÞ ðiÞ ðiÞ ðiÞ ~ . Set ~ ~B , ~ and W g ¼ mg  W K h

h

qh;i ¼ Step 4. If

h

ðiÞ

ðiÞ

ðiÞ

ðiÞ

ðgh ÞT gh ð~ gh Þ T gh

ðiþ1Þ T ðiþ1Þ Þ gh ð0Þ ð0Þ

ðgh

ðgh ÞT gh

;

ðiþ1Þ

uh

ðiÞ

ðiÞ

¼ uh  qh;i gh ; ðiþ1Þ

6 e, then take uh ¼ uh

ðiþ1Þ

gh

ðiÞ

ðiÞ

¼ gh  qh;i ~gh :

; else, go to Step 5.

Step 5. Set i = i + 1, and go to Step 3. Remark 3.1. The same stiffness matrix in the linear (3.25)–(3.28) depends solely on the geometry of X in iterative process, which is to get some of convenience to ours algorithm and can save substantial computational work because the stiffness matrix is simply computed once in iterative process. 4. Error estimators The convergence and stability of gradient type algorithm are classical for quadratic objective function within optimization framework [16]. Hence we are only to consider the local convergence between Algorithms 1 and 2 at every iteration. In order to investigate the convergence of the proposed algorithm, at every iteration we evaluate the local ðiþ1Þ ðiþ1Þ error kuh  uðiþ1Þ kL2 ðCÞ (see the objective function (1.1)), where u(i+1), uh respectively come from AlgoPof ðiþ1Þ N B ðiþ1Þ rithm 1 and 2 after i + 1 iterations, and uh ¼ j¼1 uh;j bN I þj . For the sake of simplicity, we assume that ð0Þ uh ¼ uð0Þ ¼ 0. then some results of literatures [1,4,5,18] can directly be used. ð0Þ

ð0Þ

Theorem 4.1. Let u(0) and w(0) be the solutions of (3.10), and uh , wh be the solutions of (3.21), if the finite element spaces contain polynomials of degree m P 2 and w(0) 2 Hm+1(X), then the following error estimates hold ð0Þ

kwð0Þ  wh kH 1 ðXÞ 6 Chm kwð0Þ kH mþ1 ðXÞ ; ð0Þ

kw



ð0Þ wh kH 1 ðXÞ

ð0Þ

þ hku



ð0Þ uh kL2 ðXÞ

ð4:1Þ m

ð0Þ

6 Ch kw kH mþ1 ðXÞ ;

ð4:2Þ

where C is positive constant, independent of h. Proof. The proof of this theorem is straightforward (see [1,5,18]).

h

Remark 4.1. For polygonal domain X, the solution does not have the required in (4.1) and (4.2) regularity for m P 3 and the estimate does remain valid. for X a convex polygonal domain and m P 3 instead of (4.1) and (4.2) we have (see [1]) ð0Þ

kwð0Þ  wh kH 1 ðXÞ 6 Ch2þs kwð0Þ kH 3þs ðXÞ ; ð0Þ

ð0Þ

kwð0Þ  wh kH 1 ðXÞ þ hkuð0Þ  uh kL2 ðXÞ 6 Ch2þs kwð0Þ kH 3þs ðXÞ ; where 0 < s 6 1 depends on the maximal interior angle of the boundary C. For example, if X is a rectangular domain, then s = 1 and the error has optimal convergence rate [1]. For the sake of convenience of depiction, let X = H1(X), M ¼ H 10 ðXÞ, and the product space Q = X · M,with element of the form q = {w(0), k(0)}, which is equipped with the natural norm 2

2

1

kqkQ ¼ ðkwð0Þ kX þ kkð0Þ kM Þ2 ; then we have the following result.

1438

L. Bingjie, L. Sanyang / Applied Mathematics and Computation 186 (2007) 1429–1440 ð0Þ

ð0Þ

Theorem 4.2. Let w(0) and k(0) be the solutions of (3.11), and wh , kh be the solutions of (3.22), the following error estimate holds ð0Þ

ð0Þ

ð0Þ 

kðwð0Þ  wh ; kð0Þ  kh ÞkQ 6 Cf inf kðwð0Þ  rh ; kð0Þ  lh ÞkQ þ kwð0Þ  wh kM g; rh ;lh

ð4:3Þ

where C is positive constant, independent of h, and k  kM denotes the norm of the dual space of M. Proof. We define a bilinear operator A(q, p) on Q by Z Z Z Aðq; pÞ :¼ wð0Þ r dX  rr  rkð0Þ dX  rwð0Þ  rl dX: X

X

X

The variational (3.11) can be written in compact form as Aðq; pÞ ¼ F ðpÞ;

8p ¼ ðr; lÞ 2 Q;

where the linear functional F(Æ) defined by F ðpÞ ¼ F ðfr; lgÞ :¼ It is easy to prove that the following condition

ð4:4Þ

R

ð0Þ

X

ðw

 wz Þl dX.

jAðq; pÞj 6 C A kqkQ kpkQ ; and the stronger inf–sup type inequality for A(q, p) ( ) jAðq; pÞj inf sup P g: q2Q p2Q kqkQ kpkQ

ð4:5Þ

Let Qh = Xh · Mh  Q, Then approximations qh 2 Qh are determined by Aðqh ; ph Þ ¼ F ðph Þ;

8ph 2 Qh ;

ð0Þ ð0Þ fwh ; kh g,

R

ð0Þ ð0Þ frh ; lh g,

ð4:6Þ ð0Þ ðwh X

where qh ¼ ph ¼ F ðph Þ ¼  wz Þlh dX. This discretization is stable under the appropriate choices of Xh and Mh, in which a discrete analogue (4.5) is fulfilled by the same argument as used above ( ) jAðqh ; ph Þj inf sup P g: ð4:7Þ qh 2Qh ph 2Qh kqh kQ kp h kQ For any q 2 Q, we define the orthographic projection operator Ph such that Aðq  Ph q; ph Þ ¼ 0;

8ph 2 Qh ;

kq  Ph qkQ 6 C inf fkwð0Þ  qh 2Qh

ð0Þ wh kV

ð4:8Þ þ kkð0Þ 

ð0Þ kh kM g:

ð4:9Þ

By means of (4.4) and (4.6), we can get ð0Þ

Aðq  qh ; ph Þ ¼ ðwð0Þ  wh ; lh ÞL2 ðXÞ : By the estimate (4.7), we have gkPh q  qh kQ 6 sup

ph 2Qh

jAðPh q  qh ; ph Þj jAðq  qh ; ph Þj ð0Þ  ¼ sup 6 C A kwð0Þ  wh kM : kph kQ kph kQ ph 2Qh

We conclude the asserted estimate by noting the triangle inequality and (4.9). This completes the proof. ð0Þ

h

ð0Þ

Theorem 4.3. Let w(0) and k(0) be the solutions of (3.11), and wh , kh be the solutions of (3.22), if the finite element spaces contain polynomials of degree m P 2 and w(0) 2 Hm+1(X), k(0) 2 Hm+1(X), then the following error estimates hold: ð0Þ

kkð0Þ  kh kH 1 ðXÞ 6 Chm ðkwð0Þ kH mþ1 ðXÞ þ kkð0Þ kH mþ1 ðXÞ Þ; kkð0Þ 

ð0Þ kh kH 1 ðXÞ

þ hkwð0Þ 

ð0Þ wh kL2 ðXÞ

6 Chm ðkwð0Þ kH mþ1 ðXÞ þ kkð0Þ kH mþ1 ðXÞ Þ;

where C is positive constant, independent of h.

ð4:10Þ ð4:11Þ

L. Bingjie, L. Sanyang / Applied Mathematics and Computation 186 (2007) 1429–1440

1439

Proof. The proof of this theorem is straightforward by using the results of Theorem 4.1 and 4.2 (or see [4,5,18]). h ðiÞ

ðiÞ

Theorem 4.4. Assume that uðiÞ ¼ uh , gðiÞ ¼ gh hold in boundary C after i iterations in Algorithms 1 and 2 (namely the i iteration in Algorithm 2 is exact relative to Algorithm 1). Let u(i+1) be the control function after ðiþ1Þ i + 1 iterations in Algorithm 1, and uh be the discrete control function expressed by means of (3.28) in Algomþ1 ðiÞ ðiÞ ~ ~ 2 H ðXÞ, k 2 H mþ1 ðXÞ, then we obtain the local truncated error estimate as follows: rithm 2, and w ðiþ1Þ

kuðiþ1Þ  uh

3

kL2 ðCÞ 6 Chm2 ðk~ wðiÞ kH mþ1 ðXÞ þ k~ kðiÞ kH mþ1 ðXÞ Þ:

ð4:12Þ

Proof. By (3.12), (3.13) and (3.23), (3.24), similar to Theorems 4.1 and 4.3, we can conclude the error estimates as follows: ðiÞ k~ kðiÞ  ~ kh kH 1 ðXÞ 6 Chm ðk~ wðiÞ kH mþ1 ðXÞ þ k~ kðiÞ kH mþ1 ðXÞ Þ;

~ðiÞ

kk 

ðiÞ ~ kh kH 1 ðXÞ

~ ðiÞ

þ hkw 

~ ðiÞ k 2 w h L ðXÞ

m

ðiÞ

ð4:13Þ ~ðiÞ

6 Ch ðk~ w kH mþ1 ðXÞ þ kk kH mþ1 ðXÞ Þ:

ð4:14Þ

By step 3 in Algorithm 1 and step 3 in Algorithm 2, we have ðiÞ ~ ðiÞ  c w ~ ðiÞ k~ gðiÞ  g~h kL2 ðCÞ 6 kc0 w 0 h kL2 ðCÞ :

ð4:15Þ

Because of the linearity and continuity of the gradient operator, there exists a constant C such that ðiþ1Þ

kuðiþ1Þ  uh

ðiþ1Þ

kL2 ðCÞ 6 Ckgðiþ1Þ  gh

kL2 ðCÞ :

ð4:16Þ

Because the optimum iterative steps qi, qi,h are bounded, we have ðiþ1Þ

kgðiþ1Þ  gh

ðiÞ

ðiÞ

kL2 ðCÞ 6 kgðiÞ  gh kL2 ðCÞ þ maxfqi ; qi;h gk~ gðiÞ  g~h kL2 ðCÞ :

ð4:17Þ

Given a quasi-uniform subdivision of X, the following inverse inequality holds 1

ðiÞ

ðiÞ

2 ~ ðiÞ ~ ðiÞ  c w ~ ~ kc0 w 0 h kL2 ðCÞ 6 Ch kw  wh kL2 ðXÞ :

ð4:18Þ

We can conclude the asserted estimate (4.12) by noting (4.14)–(4.18) and the assumption conditions of this theorem. This completes the proof. h Remark 4.2. For a convex polygonal domain X and m P 3, similar to Remark 4.1, we have ðiþ1Þ

kuðiþ1Þ  uh

1

kL2 ðCÞ 6 Ch2þs ðk~ wðiÞ kH 3þs ðXÞ þ k~ kðiÞ kH 3þs ðXÞ Þ;

where 0 < s 6 1 depends on the maximal interior angle of the boundary C. If X is a rectangular domain, we get ðiþ1Þ

kuðiþ1Þ  uh

3 kL2 ðCÞ 6 Ch2 ðk~ wðiÞ kH 4 ðXÞ þ k~ kðiÞ kH 4 ðXÞ Þ:

Acknowledgement This research is supported by the National Science Foundation of China under Grant No. 69972036. References [1] A.B. Andreev, R.D. Lazarov, M.R. Rachrva, Postprocessing and higher order convergence of the mixed finite element approximations of bi-harmonic eigenvalue problems, J. Comput. Appl. Math 182 (2005) 333–349. [2] A. Borzi, K. Kunisch, D.Y. Kwak, Accuracy and convergence properties of the finite difference multi-grid solution of an optimal control optimality system, SIAM J. Control Optim. 41 (2003) 477–1497. [3] A. Borzi, K. Kunisch, A multigrid scheme for elliptic constrained optimal control problems, Comput. Optim. Appl. 31 (2005) 309– 333.

1440

L. Bingjie, L. Sanyang / Applied Mathematics and Computation 186 (2007) 1429–1440

[4] D. Boffi, C. Lovadina, Analysis of new augmented Lagrangian formulations for mixed finite element schemes, Numer. Math. 75 (1997) 405–419. [5] P.G. Ciarlet, P.A. Raviart, A mixed finite method for the bi-harmonic equation, in: C. de Boor (Ed.), mathematical Aspects of Finite Elements in Partial Differential Equations, Academic Press, 1974, pp. 125–145. [6] H.M. Gu, X.N. Wu, Alternating-direction iterative technique for mixed finite element methods for Stokes equations, Appl. Math. Comput. 162 (3) (2005) 1035–1047. [7] M.R. Hanisch, Two-level additive Schwarz preconditioners for fourth-order mixed methods, Electron. Trans. Numer. Anal. 22 (2006) 1–16. [8] N. Hinze, A variational discretization concept in control constrained optimization: the linear-quadratic case, Comput. Optim. Appl. 30 (2005) 45–61. [9] Y. Kim, H. Huh, New rectangular mixed finite element method for second-order elliptic problems, Appl. Math. Comput. 127 (2–3) (2002) 375–385. [10] K. Kulshreshtha, N. Nataraj, M. Jung, A parallel mixed finite element implementation of fourth-order plate bending problems in distributed memory environments, Appl. Math. Comput. 163 (1) (2005) 253–274. [11] R. Li, W.B. Liu, H.P. Ma, T. Tang, Adaptive finite element approximation for distributed elliptic optimal control problems, SIAM J. Control Optim. 41 (2003) 1321–1349. [12] B.J. Li, MRM boundary integral equation of two kinds of higher- order fundamental solution, Acta Math. Appl. Sinica 23 (2000) 534–542 (in Chinese). [13] B.J. Li, The multiple reciprocity method (MRM) boundary variation equation for boundary value problem of partial differential equation D2u  sDu + k2u = 0, Acta Math. Appl. Sinica 25 (2002) 660–666 (in Chinese). [14] S.M. Mefire, Mixed finite element and boundary element approximation in 3D magnetostatics for computation of the magnetic induction, Appl. Math. Comput. 125 (2–3) (2002) 399–421. [15] J. Ortega, E. Zuazua, Generic simplicity of the spectrum and stabilization for a plate equation, SIAM J. Control Optim. 39 (2001) 1585–1614. [16] A.M. Ramos, R. Glowinski, J. Periaux, Nash equilibria for the multiobjective control of linear partial differential equations, J. Optim. Theory Appl. 112 (2002) 457–498. [17] D.J. Silvester, M.D. Mihajlovic, A black-box multigrid preconditioner for the biharmonic equation, Numer. Math. 44 (2004) 151– 163. [18] Y.D. Yang, J.B. Gao, Remarks on the Ciarlet–Raviart mixed finite element, Electron. Trans. Numer. Anal. 4 (1996) 158–164. [19] D. Yang, Iterative schemes for mixed finite element methods with applications to elasticity and compressible flow problems, Numer. Math. 83 (2002) 177–200.