# Adaptive Mixed Finite Element Methods for Parabolic Optimal Control Problems

Authors: Zuliang Lu
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

We will investigate the adaptive mixed finite element methods for
parabolic optimal control problems. The state and the costate are approximated by the lowest-order Raviart-Thomas mixed finite element spaces, and the control is approximated by piecewise constant elements. We derive a posteriori error estimates of the mixed finite element solutions for optimal control problems. Such a posteriori error estimates can be used to construct more efficient and reliable adaptive mixed finite element method for the optimal control problems. Next we introduce an adaptive algorithm to guide the mesh refinement. A numerical example is given to demonstrate our theoretical results.


---

## Body

## 1. Introduction

Optimal control problems are very important models in science and engineering numerical simulation. Finite element method of optimal control problems plays an important role in numerical methods for these problems. Let us mention two early papers devoted to linear optimal control problems by Falk [1] and Geveci [2]. Knowles was concerned with standard finite element approximation of parabolic time optimal control problems in [3]. In [4] Gunzburger and Hou investigated the finite element approximation of a class of constrained nonlinear optimal control problems. For quadratic optimal control problem governed by linear parabolic equation, Liu and Yan derived a posteriori error estimates for both the state and the control approximation in [5]. Systematic introductions of the finite element method for optimal control problems can be found in [6–10].Adaptive finite element approximation was the most important means of boosting the accuracy and efficiency of finite element discretization. The literature in this aspect was huge, see, for example, [11, 12]. Adaptive finite element method was widely used in engineering numerical simulation. There has been extensive studies on adaptive finite element approximation for optimal control problems. In [13], the authors have introduced some basic concept of adaptive finite element discretization for optimal control of partial differential equations. A posteriori error estimators for distributed elliptic optimal control problems were contained in Li et al. [14]. Recently an adaptive finite element method for the estimation of distributed parameter in elliptic equation was discussed by Feng et al. [15]. Note that all the above works aimed at standard finite element method.In many control problems, the objective functional contains the gradient of the state variables. Thus, the accuracy of the gradient is important in numerical discretization of the coupled state equations. When the objective functional contains the gradient of the state variable, mixed finite element methods should be used for discretization of the state equation with which both the scalar variable and its flux variable can be approximated in the same accuracy. In [16–20] we have done some primary works on a priori error estimates and superconvergence for linear optimal control problems by mixed finite element methods. We considered a posteriori error estimates of mixed finite element methods for quadratic and general optimal control problems in [21–23].In [24], the authors discussed the mixed finite element approximation for general optimal control problems governed by parabolic equation. And then, they derived a posteriori error estimates of mixed finite element solution. In this paper, we study the adaptive mixed finite element methods for the parabolic optimal control problems. We construct the mixed finite element discretization for the original problems and derive a useful posteriori error indicators. Furthermore, we provide an adaptive algorithm to guide the multimesh refinement. Finally, a numerical experiment shows that this algorithm works very well with the adaptive multimesh discretization.The plan of this paper is as follows. In the next section, we construct the mixed finite element discretization for the parabolic optimal control problems. Then, we derive a posteriori error estimates for the mixed finite element solutions in Section3. Next, we introduce an adaptive algorithm to guide the mesh refinement in Section 4. Finally, a numerical example is given to demonstrate our theoretical results in Section 5.


## 2. Mixed Methods of Optimal Control Problems

In this section, we investigate the mixed finite element approximation for parabolic optimal control problems. We adopt the standard notationWm,p(Ω) for Sobolev spaces on Ω with a norm ∥·∥m,p given by ∥v∥m,pp=∑|α|≤m∥Dαv∥Lp(Ω)p, a seminorm |·|m,p given by |v|m,pp=∑|α|=m∥Dαv∥Lp(Ω)p. We set W0m,p(Ω)={v∈Wm,p(Ω):v|∂Ω=0}. For p=2, we denote Hm(Ω)=Wm,2(Ω),H0m(Ω)=W0m,2(Ω), and ∥·∥m=∥·∥m,2,∥·∥=∥·∥0,2.We denote byLs(0,T;Wm,p(Ω)) the Banach space of all Ls integrable functions from J into Wm,p(Ω) with norm ∥v∥Ls(J;Wm,p(Ω))=(∫0T∥v∥Wm,p(Ω)sdt)1/sfors∈[1,∞) and the standard modification for s=∞. The details can be found in [25].The parabolic optimal control problems that we are interested in are as follows:(2.1)min⁡u∈K⊂U{∫0T(g1(p)+g2(y)+j(u))dt},yt(x,t)+div⁡p(x,t)=f+u(x,t),x∈Ω,p(x,t)=-A(x)∇y(x,t),x∈Ω,y(x,t)=0,x∈∂Ω,t∈J,y(x,0)=y0(x),x∈Ω,
where the bounded open set Ω⊂ℝ2 is a convex polygon with the boundary ∂Ω. Let K be a closed convex set in U=L2(J;L2(Ω)), f∈L2(J;L2(Ω)), J=[0,T], y0(x)∈H01(Ω). Furthermore, we assume that the coefficient matrix A(x)=(ai,j(x))2×2∈L∞(Ω;ℝ2×2) is a symmetric 2×2-matrix and there is a constant c>0 satisfying for any vector X∈ℝ2, XtAX≥c∥X∥ℝ22. j′ is positive, g1′, g2′, and j′ are locally Lipschitz on L2(Ω)2, W, U, and that there is a c>0 such that (j′(u)-j′(ũ),u-ũ)≥c∥u-ũ∥0, forallu,ũ∈U.Now we will describe the mixed finite element discretization of parabolic optimal control problems (2.1). Let (2.2)V=H(div⁡;Ω)={v∈(L2(Ω))2,div⁡v∈L2(Ω)},W=L2(Ω).
The Hilbert space V is equipped with the following norm: (2.3)‖v‖div⁡=‖v‖H(div⁡;Ω)=(‖v‖0,Ω2+‖div⁡v‖0,Ω2)1/2.We recast (2.1) as the following weak form: find (p,y,u)∈V×W×K such that(2.4)min⁡u∈K⊂U{∫0T(g1(p)+g2(y)+j(u))dt},(2.5)(A-1p,v)-(y,div⁡v)=0,∀v∈V,(2.6)(yt,w)+(div⁡p,w)=(f+u,w),∀w∈W,(2.7)y(x,0)=y0(x),∀x∈Ω.Similar to [26], the optimal control problems (2.4)–(2.7) have a unique solution (p,y,u), and a triplet (p,y,u) is the solution of (2.4)–(2.7) if and only if there is a costate (q,z)∈V×W such that (p,y,q,z,u) satisfies the following optimality conditions:(2.8)(A-1p,v)-(y,div⁡v)=0,∀v∈V,(2.9)(yt,w)+(div⁡p,w)=(f+u,w),∀w∈W,(2.10)y(x,0)=y0(x),∀x∈Ω,(2.11)(A-1q,v)-(z,div⁡v)=-(g1′(p),v),∀v∈V,(2.12)-(zt,w)+(div⁡q,w)=(g2′(y),w),∀w∈W,(2.13)z(x,T)=0,∀x∈Ω,(2.14)∫0T(j′(u)+z,ũ-u)Udt≥0,∀ũ∈K,
where (·,·)U is the inner product of U. In the rest of the paper, we will simply write the product as (·,·) whenever no confusion should be caused.Let𝒯h be regular triangulation of Ω. They are assumed to satisfy the angle condition which means that there is a positive constant C such that, for all τ∈𝒯h, C-1hτ2≤|τ|≤Chτ2, where |τ| is the area of τ, hτ is the diameter of τ and h=maxhτ. In addition C or c denotes a general positive constant independent of h.LetVh×Wh⊂V×W denote the Raviart-Thomas space [27] of the lowest order associated with the triangulation 𝒯h of Ω. Pk denotes the space of polynomials of total degree at most k. Let V(τ)={v∈P02(τ)+x·P0(τ)}, W(τ)=P0(τ). We define (2.15)Vh∶={vh∈V:∀τ∈Th,vh|τ∈V(τ)},Wh∶={wh∈W:∀τ∈Th,wh|τ∈W(τ)},Kh∶={ũh∈K:∀τ∈Th,ũh|τ=constant}.The mixed finite element discretization of (2.4)–(2.7) is as follows: compute (ph,yh,uh)∈Vh×Wh×Kh such that(2.16)min⁡uh∈Kh{∫0T(g1(ph)+g2(yh)+j(uh))dt},(A-1ph,vh)-(yh,div⁡vh)=0,∀vh∈Vh,(yht,wh)+(div⁡ph,wh)=(f+uh,wh),∀wh∈Wh,yh(x,0)=y0h(x),∀x∈Ω,
where y0h(x)∈Wh is an approximation of y0. The optimal control problem (2.16) again has a unique solution (ph,yh,uh), and a triplet (ph,yh,uh) is the solution of (2.16) if and only if there is a costate (qh,zh)∈Vh×Wh such that (ph,yh,qh,zh,uh) satisfies the following optimality conditions:(2.17)(A-1ph,vh)-(yh,div⁡vh)=0,∀vh∈Vh,(yht,wh)+(div⁡ph,wh)=(f+uh,wh),∀wh∈Wh,yh(x,0)=y0(x),∀x∈Ω,(A-1qh,vh)-(zh,div⁡vh)=-(g1′(ph),vh),∀vh∈Vh,-(zht,wh)+(div⁡qh,wh)=(g2′(yh),wh),∀wh∈Wh,zh(x,T)=0,∀x∈Ω,(j′(uh)+zh,ũh-uh)≥0,∀ũh∈Kh.We now consider the fully discrete approximation for the semidiscrete problem. LetΔt>0, N=T/Δt∈ℤ, and ti=iΔt, i∈ℤ. Also, let (2.18)ψi=ψi(x)=ψ(x,ti),dtψi=ψi-ψi-1Δt.
For i=1,2,…,N, construct the finite element spaces Vhi∈V, Whi∈W (similar as Vh). Similarly, construct the finite element spaces Khi∈Kh with the mesh 𝒯hi. Let hτi denote the maximum diameter of the element τi in 𝒯hi. Define mesh functions τ(·) and mesh size functions hτ(·) such that τ(t)|t∈(ti-1,ti]=τi, hτ(t)|t∈(ti-1,ti]=hτi. For ease of exposition, we will denote τ(t) and hτ(t) by τ and hτ, respectively.The following fully discrete approximation scheme is to find(phi,yhi,uhi)∈Vhi×Whi×Khi,i=1,2,…,N, such that(2.19)min⁡uhi∈Khi{∑i=1N∫ti-1ti(g1(phi)+g2(yhi)+j(uhi))},(2.20)(A-1phi,vh)-(yhi,div⁡vh)=0,∀vh∈Vhi,(2.21)(dtyhi,wh)+(div⁡phi,wh)=(f(x,ti)+uhi,wh),∀wh∈Whi,(2.22)yh0(x,0)=y0h(x),∀x∈Ω.
It follows that the optimal control problems (2.19)–(2.22) have a unique solution (phi,yhi,uhi),i=1,2,…,N, and that a triplet (phi,yhi,uhi)∈Vhi×Whi×Khi,  i=1,2,…,N, is the solution of (2.19)–(2.22) if and only if there is a costate (qhi-1,zhi-1)∈Vhi×Whi such that (phi,yhi,qhi-1,zhi-1,uhi)∈(Vhi×Whi)2×Khi satisfies the following optimality conditions:(2.23)(A-1phi,vh)-(yhi,div⁡vh)=0,∀vh∈Vhi,(2.24)(dtyhi,wh)+(div⁡phi,wh)=(f(x,ti)+uhi,wh),∀wh∈Whi,(2.25)yh0(x,0)=y0h(x),∀x∈Ω,(2.26)(A-1qhi-1,vh)-(zhi-1,div⁡vh)=-(g1′(phi),vh),∀vh∈Vhi,(2.27)-(dtzhi,wh)+(div⁡qhi-1,wh)=(g2′(yhi),wh),∀wh∈Whi,(2.28)zhN(x,T)=0,∀x∈Ω,(2.29)(j′(uhi)+zhi-1,ũh-uhi)≥0,∀ũh∈Khi.Fori=1,2,…,N, let (2.30)Yh|(ti-1,ti]=((ti-t)yhi-1+(t-ti-1)yhi)Δt,Zh|(ti-1,ti]=((ti-t)zhi-1+(t-ti-1)zhi)Δt,Ph|(ti-1,ti]=((ti-t)phi-1+(t-ti-1)phi)Δt,Qh|(ti-1,ti]=((ti-t)qhi-1+(t-ti-1)qhi)Δt,Uh|(ti-1,ti]=uhi.For any functionw∈C(J;L2(Ω)), let (2.31)ŵ(x,t)|t∈(ti-1,ti]=w(x,ti),w̃(x,t)|t∈(ti-1,ti]=w(x,ti-1).Then the optimality conditions (2.23)–(2.29) can be restated as.(2.32)(A-1P̂h,vh)-(Ŷh,div⁡vh)=0,∀vh∈Vhi,(2.33)(Yht,wh)+(div⁡P̂h,wh)=(f̂+Uh,wh),∀wh∈Whi,(2.34)Yh(x,0)=y0h(x),∀x∈Ω,(2.35)(A-1Q̃h,vh)-(Z̃h,div⁡vh)=-(g1′(P̂h),vh),∀vh∈Vhi,(2.36)-(Zht,wh)+(div⁡Q̃h,wh)=(g2′(Ŷh),wh),∀wh∈Whi,(2.37)Zh(x,T)=0,∀x∈Ω,(2.38)(j′(Uh)+Z̃h,ũh-Uh)≥0,∀ũh∈Khi.In the rest of the paper, we will use some intermediate variables. For any control functionUh∈K, we first define the state solution (p(Uh),y(Uh),q(Uh),z(Uh)) which satisfies(2.39)(A-1p(Uh),v)-(y(Uh),div⁡v)=0,∀v∈V,(2.40)(yt(Uh),w)+(div⁡p(Uh),w)=(f+Uh,w),∀w∈W,(2.41)y(Uh)(x,0)=y0(x),∀x∈Ω,(2.42)(A-1q(Uh),v)-(z(Uh),div⁡v)=-(g1′(p(Uh)),v),∀v∈V,(2.43)-(zt(Uh),w)+(div⁡q(Uh),w)=(g2′(y(Uh)),w),∀w∈W,(2.44)z(Uh)(x,T)=0,∀x∈Ω.


## 3. A Posteriori Error Estimates

In this section we study a posteriori error estimates of the mixed finite element approximation for the parabolic optimal control problems. Fixed givenu∈K, let ℳ1, ℳ2 be the inverse operators of the state equation (2.6), such that p(u)=ℳ1u and y(u)=ℳ2u are the solutions of the state equations (2.6). Similarly, for given Uh∈Kh, Ph(Uh)=ℳ1hUh, Yh(Uh)=ℳ2hUh are the solutions of the discrete state equation (2.33). Let(3.1)F(u)=g1(M1U)+g2(M2U)+j(u),Fh(Uh)=g1(M1hUh)+g2(M2hUh)+j(Uh).
It can be shown that(3.2)(F′(u),v)=(j′(u)+z,v),(F′(Uh),v)=(j′(Uh)+z(Uh),v),(Fh′(Uh),v)=(j′(Uh)+Z̃h,v).It is clear thatF and Fh are well defined and continuous on K and Khi. Also the functional Fh can be naturally extended on K. Then (2.4) and (2.19) can be represented as(3.3)min⁡u∈K{F(u)},(3.4)min⁡Uh∈Khi{Fh(Uh)}.In many application,F(·) is uniform convex near the solution u. The convexity of F(·) is closely related to the second order sufficient conditions of the optimal control problems, which are assumed in many studies on numerical methods of the problem. For instance, in many applications, u→g1(ℳ1U) and u→g2(ℳ2U) are convex. Then there is a constant c>0, independent of h, such that(3.5)∫0T(F′(u)-F′(Uh),u-Uh)U≥c‖u-Uh‖L2(J;L2(Ω))2.Theorem 3.1.
Letu and Uh be the solutions of (3.3) and (3.4), respectively. Assume that Khi⊂K. In addition, assume that (Fh′(Uh))|τ∈Hs(τ),forallτ∈𝒯h,(s=0,1), and there is a vh∈Khi such that
(3.6)|(Fh′(Uh),vh-u)|≤C∑τ∈Thhτ‖Fh′(Uh)‖Hs(τ)‖u-Uh‖L2(τ)s.
Then one has
(3.7)‖u-Uh‖L2(J;L2(Ω))2≤Cη12+C‖z(Uh)-Z̃h‖L2(J;L2(Ω))2,
where
(3.8)η12=∫0T∑τ∈Thhτ1+s‖j′(Uh)+Z̃h‖H1(τ)1+s.Proof.
It follows from (3.3) and (3.4) that
(3.9)∫0T(F′(u),u-v)≤0,∀v∈K,∫0T(Fh′(Uh),Uh-vh)≤0,∀vh∈Khi⊂K.
Then it follows from assumptions (3.5), (3.6) and Schwartz inequality that
(3.10)c‖u-Uh‖L2(J;L2(Ω))2≤∫0T(F′(u)-F′(Uh),u-Uh)≤∫0T{(Fh′(Uh),vh-u)+(Fh′(Uh)-F′(Uh),u-Uh)}≤C∫0T{∑τ∈Thhτ1+s‖Fh′(Uh)‖Hs(τ)1+s+‖Fh′(Uh)-F′(Uh)‖L2(Ω)2}+c2‖u-Uh‖L2(J;L2(Ω))2.
It is not difficult to show
(3.11)Fh′(Uh)=j′(Uh)+Z̃h,F′(Uh)=j′(Uh)+z(Uh),
where z(Uh) is defined in (2.39)–(2.44). Thanks to (3.11), it is easy to derive
(3.12)‖Fh′(Uh)-F′(Uh)‖L2(Ω)≤C‖Z̃h-z(Uh)‖L2(Ω).
Then by estimates (3.10) and (3.12) we can prove the requested result (3.7).In order to estimate the a posteriori error of the mixed finite element approximation solution, we will use the following dual equations:(3.13)-φt-div⁡(A∇φ)=G,x∈Ω,t∈(0,T],φ|∂Ω=0,t∈J,φ(x,T)=0,x∈Ω,(3.14)ψt-div⁡(A*∇ψ)=G,x∈Ω,t∈(0,T],ψ|∂Ω=0,t∈J,ψ(x,0)=0,x∈Ω.The following well-known stability results are presented in [28].Lemma 3.2.
Letφ and ψ be the solutions of (3.13), and (3.14), respectively. Then, for v=φ or v=ψ,
(3.15)‖v‖L∞(J;L2(Ω))≤C‖G‖L2(J;L2(Ω)),‖∇v‖L2(J;L2(Ω))≤C‖G‖L2(J;L2(Ω)),‖D2v‖L2(J;L2(Ω))≤C‖G‖L2(J;L2(Ω)),‖vt‖L2(J;L2(Ω))≤C‖G‖L2(J;L2(Ω)),
where D2v=∂2v/∂xi∂xj,1≤i,j≤2.Now, we are able to derive the main result.Theorem 3.3.
Let(Yh,Ph,Zh,Qh,Uh) and (y(Uh),p(Uh),z(Uh),q(Uh),Uh) be the solutions of  (2.32)–(2.38) and (2.39)–(2.44). Then,
(3.16)‖Yh-y(Uh)‖L2(J;L2(Ω))2+‖Ph-p(Uh)‖L2(J;L2(Ω))2≤Cη22,
where
(3.17)η22=‖Yht+div⁡P̂h-f̂-Uh‖L2(J;L2(Ω))2+‖(Yh-y(Uh))(x,0)‖L2(Ω)2+‖f̂-f‖L2(J;L2(Ω))2+‖(Ŷh-Yh)t‖L2(J;L2(Ω))2+‖P̂h-Ph‖L2(J;H(div⁡;Ω))2.Proof.
Lettingφ be the solution of (3.13) with G=Yh-y(Uh), we infer
(3.18)‖Yh-y(Uh)‖L2(J;L2(Ω))2=∫0T(Yh-y(Uh),F)dt=∫0T(Yh-y(Uh),-φt-div⁡(A∇φ))dt=∫0T(((Yh-y(Uh))t,φ)-(Ph-p(Uh),∇φ))dt+‖(Yh-y(Uh))(x,0)‖L2(Ω)2=∫0T((Yht-yt(Uh),φ)+(div⁡(Ph-p(Uh)),φ))dt+‖(Yh-y(Uh))(x,0)‖L2(Ω)2.
Then it follows from (2.39)-(2.40) that
(3.19)‖Yh-y(Uh)‖L2(J;L2(Ω))2=∫0T((Yht,φ)+(div⁡P̂h,φ)-(yt(Uh),φ)-(div⁡p(Uh),φ)+(div⁡(Ph-P̂h),φ))dt+‖(Yh-y(Uh))(x,0)‖L2(Ω)2=∫0T((Yht+div⁡P̂h-f̂-Uh,φ)+(f̂-f,φ)+(div⁡(Ph-P̂h),φ))dt+‖(Yh-y(Uh))(x,0)‖L2(Ω)2=∫0T((Yht+div⁡P̂h-f̂-Uh,φ)+(f̂-f,φ)+(div⁡(Ph-P̂h),φ))dt+‖(Yh-y(Uh))(x,0)‖L2(Ω)2≤C(δ)‖Yht+div⁡P̂h-f̂-Uh‖L2(J;L2(Ω))2+C(δ)‖f̂-f‖L2(J;L2(Ω))2+C(δ)‖P̂h-Ph‖L2(J;H(div⁡;Ω))2+C(δ)‖Ŷh-Yh‖L2(J;L2(Ω))2+‖(Yh-y(Uh))(x,0)‖L2(Ω)2+Cδ‖φ‖L2(J;L2(Ω))2,
where and after, δ is an arbitrary positive number, C(δ) is the constant dependent on δ-1.
Now, we are in the position of estimating the error∥Ph-p(Uh)∥L2(J;L2(Ω))2. First, we derive from (2.32)-(2.33) and (2.39)-(2.40) the following useful error equations:
(3.20)(A-1(P̂h-p(Uh)),vh)-(Ŷh-y(uh),div⁡vh)=0,(3.21)((Ŷh-y(Uh))t,wh)+(div⁡(P̂h-p(Uh)),wh)=(f̂-f,wh)-((Yh-Ŷh)t,wh),
where vh∈Vh, wh∈Wh. Choose vh=P̂h-p(Uh) and wh=Ŷh-y(Uh) as the test functions and add the two relations of (3.20)-(3.21) such that
(3.22)(A-1(P̂h-p(Uh)),P̂h-p(Uh))+((Ŷh-y(Uh))t,Ŷh-y(Uh))=(f̂-f,Ŷh-y(Uh))-((Yh-Ŷh)t,Ŷh-y(Uh)).
Then, using ϵ-Cauchy inequality, we can find an estimate as follows:
(3.23)c‖P̂h-p(Uh)‖L2(Ω)2+((Ŷh-y(Uh))t,Ŷh-y(Uh))≤C‖Ŷh-y(Uh)‖L2(Ω)2+C‖f̂-f‖L2(Ω)2+C‖(Ŷh-Yh)t‖L2(Ω)2.
Note that
(3.24)((Ŷh-y(Uh))t,Ŷh-y(Uh))=12∂∂t‖Ŷh-y(Uh)‖L2(Ω)2,
then, using the assumption on A, we verify that
(3.25)c‖P̂h-p(Uh)‖L2(Ω)2+12∂∂t‖Ŷh-y(Uh)‖L2(Ω)2≤C‖Ŷh-y(Uh)‖L2(Ω)2+C‖f̂-f‖L2(Ω)2+C‖(Ŷh-Yh)t‖L2(Ω)2.
Integrating (3.25) in time and since Ŷh(0)-y(Uh)(0)=0, applying Gronwall's lemma, we can easily obtain the following error estimate:
(3.26)‖P̂h-p(Uh)‖L2(J;L2(Ω))+‖Ŷh-y(Uh)‖L∞(J;L2(Ω))≤C‖f̂-f‖L2(J;L2(Ω))+C‖(Ŷh-Yh)t‖L2(J;L2(Ω)).
Using the triangle inequality and (3.26), we deduce that
(3.27)‖Ph-p(Uh)‖L2(J;L2(Ω))≤C(‖f̂-f‖L2(J;L2(Ω))+‖(Ŷh-Yh)t‖L2(J;L2(Ω))+‖P̂h-Ph‖L2(J;L2(Ω))).
Then letting δ be small enough, it follows from (3.18)–(3.27) that
(3.28)‖Yh-y(Uh)‖L2(J;L2(Ω))2+‖Ph-p(Uh)‖L2(J;L2(Ω))2≤Cη22.
This proves (3.16).Similarly, lettingψ be the solution of (3.14) with G=Zh-z(Uh), with the aid of (2.26)-(2.27), (2.42)-(2.43), we can conclude the following.Theorem 3.4.
Let(Yh,Ph,Zh,Qh,Uh) and (y(Uh),p(Uh),z(Uh),q(Uh),Uh) be the solutions of (2.32)–(2.38) and (2.39)–(2.44). Then,
(3.29)‖Zh-z(Uh)‖L2(J;L2(Ω))2+‖Qh-q(Uh)‖L2(J;L2(Ω))2≤C(η22+η32),
where
(3.30)η32=‖g2′(Ŷh)+Zht-div⁡(Q̃h)‖L2(J;L2(Ω))2+‖Z̃h-Zh‖L2(J;L2(Ω))2+‖(Z̃h-Zh)t‖L2(J;L2(Ω))2+‖Yh-Ŷh‖L2(J;L2(Ω))2+‖Q̃h-Qh‖L2(J;H(div⁡;Ω))2.Let(p,y,q,z,u) and (Ph,Yh,Qh,Zh,Uh) be the solutions of (2.8)–(2.14) and (2.32)–(2.38), respectively. We decompose the errors as follows: (3.31)p-Ph=p-p(Uh)+p(Uh)-Ph∶=ϵ1+ɛ1,y-Yh=y-y(Uh)+y(Uh)-Yh∶=r1+e1,q-Qh=q-q(Uh)+q(Uh)-Qh∶=ϵ2+ɛ2,z-Zh=z-z(Uh)+z(Uh)-Zh∶=r2+e2.From (2.8)–(2.13) and (2.39)–(2.44), we derive the error equations:(3.32)(A-1ϵ1,v)-(r1,div⁡v)=0,∀v∈V,(3.33)(r1t,w)+(div⁡ϵ1,w)=(u-Uh,w),∀w∈W,(3.34)(A-1ϵ2,v)-(r2,div⁡v)=-(g1′(p)-g1′(p(Uh))v),∀v∈V,(3.35)(r2t,w)+(div⁡ϵ2,w)=(g2′(y)-g2′(y(Uh)),w),∀w∈W.Theorem 3.5.
There is a constantC>0, independent of h, such that
(3.36)‖ϵ1‖L2(J;L2(Ω))+‖r1‖L2(J;L2(Ω))≤C‖u-Uh‖L2(J;L2(Ω)),(3.37)‖ϵ2‖L2(J;L2(Ω))+‖r2‖L2(J;L2(Ω))≤C‖u-Uh‖L2(J;L2(Ω)).Proof.
Part I
Choosev=ϵ1 and w=r1 as the test functions and add the two relations of (3.32)-(3.33), then we have
(3.38)(A-1ϵ1,ϵ1)+(r1t,r1)=(u-Uh,r1).
Then, using the Cauchy inequality, we can find an estimate as follows:
(3.39)(A-1ϵ1,ϵ1)+(r1t,r1)≤C(‖r1‖L2(Ω)2+‖u-Uh‖L2(Ω)2).
Note that
(3.40)(r1t,r1)=12∂∂t‖r1‖L2(Ω)2,
then, using the assumption on A, we can obtain that
(3.41)‖ϵ1‖L2(Ω)2+12∂∂t‖r1‖L2(Ω)2≤C(‖r1‖L2(Ω)2+‖u-Uh‖L2(Ω)2).
Integrating (3.41) in time and since r1(0)=0, applying the Gronwall's lemma, we can easily obtain the following error estimate
(3.42)‖ϵ1‖L2(J;L2(Ω))2+‖r1‖L2(J;L2(Ω))2≤C‖u-Uh‖L2(J;L2(Ω))2.
This implies (3.36).
Part II
Similarly, choosev=ϵ2 and w=r2 as the test functions and add the two relations of (3.34)-(3.35), then we can obtain that
(3.43)(A-1ϵ2,ϵ2)+(r2t,r2)=(g2′(y)-g2′(y(Uh)),r2)-(g1′(p)-g1′(p(Uh)),ϵ2).
Then, using the Cauchy inequality, we can find an estimate as follows:
(3.44)(A-1ϵ2,ϵ2)+(r2t,r2)≤C(‖r1‖L2(Ω)2+‖r2‖L2(Ω)2+‖ϵ1‖L2(Ω)2)+c2‖ϵ2‖L2(Ω)2.
Note that
(3.45)(r2t,r2)=12∂∂t‖r2‖L2(Ω)2,
then, using the assumption on A, we verify that
(3.46)‖ϵ2‖L2(Ω)2+12∂∂t‖r2‖L2(Ω)2≤C(‖r1‖L2(Ω)2+‖r2‖L2(Ω)2+‖ϵ1‖L2(Ω)2).
Integrating (3.46) in time and since r2(T)=0, applying the Gronwall's lemma, we can easily obtain the following error estimate
(3.47)‖ϵ2‖L2(J;L2(Ω))2+‖r2‖L2(J;L2(Ω))2≤C‖u-Uh‖L2(J;L2(Ω))2.
Then (3.37) follows from (3.47) and the previous statements immediately.Collecting Theorems3.1–3.5, we can derive the following results.Theorem 3.6.
Let(p,y,q,z,u) and (Ph,Yh,Qh,Zh,Uh) be the solutions of (2.8)–(2.14) and (2.32)–(2.38), respectively. In addition, assume that (j′(Uh)+Z̃h)|τ∈Hs(τ),forallτ∈𝒯h,(s=0,1), and that there is a vh∈Kh such that
(3.48)|(j′(Uh)+Z̃h,vh-u)|≤C∑τ∈Thhτ‖j′(Uh)+Z̃h‖Hs(τ)‖u-Uh‖L2(τ)s.
Then one has that, forallt∈(0,T],
(3.49)‖u-Uh‖L2(J;L2(Ω))2+‖y-Yh‖L2(J;L2(Ω))2+‖p-Ph‖L2(J;L2(Ω))2+‖z-Zh‖L2(J;L2(Ω))2+‖q-Qh‖L2(J;L2(Ω))2≤C∑i=13ηi2,
where η1, η2, and η3 are defined in Theorems 3.1–3.4.


## 4. An Adaptive Algorithm

In the section, we introduce an adaptive algorithm to guide the mesh refine process. A posteriori error estimates which have been derived in Section3 are used as an error indicator to guide the mesh refinement in adaptive finite element method.Now, we discuss the adaptive mesh refinement strategy. The general idea is to refine the mesh such that the error indicator likeη is equally distributed over the computational mesh. Assume that an a posteriori error estimator η has the form of η2=∑τiητi2, where τi is the finite elements. At each iteration, an average quantity of all ητi2 is calculated, and each ητi2 is then compared with this quantity. The element τi is to be refined if ητi2 is larger than this quantity. As ητi2 represents the total approximation error over τi, this strategy makes sure that higher density of nodes is distributed over the area where the error is higher.Based on this principle, we define an adaptive algorithm of the optimal control problems (2.1) as follows.Starting from initial triangulations𝒯h0 of Ω, we construct a sequence of refined triangulation 𝒯hj as follows. Given 𝒯hj, we compute the solutions (Ph,Yh,Qh,Zh,Uh) of the system (2.32)–(2.38) and their error estimator (4.1)ητ2=∫0T∑τ∈Thhτ1+s‖j′(Uh)+Z̃h‖H1(τ)1+s+‖Yht+div⁡P̂h-f̂-Uh‖L2(J;L2(τ))2+‖(Yh-y(Uh))(x,0)‖L2(τ)2+‖f̂-f‖L2(J;L2(τ))2+‖(Ŷh-Yh)t‖L2(J;L2(τ))2+‖P̂h-Ph‖L2(J;H(div⁡;τ))2+‖g2′(Ŷh)+Zht-div⁡(Q̃h)‖L2(J;L2(τ))2+‖Z̃h-Zh‖L2(J;L2(τ))2+‖(Z̃h-Zh)t‖L2(J;L2(τ))2+‖Yh-Ŷh‖L2(J;L2(τ))2+‖Q̃h-Qh‖L2(J;H(div⁡;τ))2,Ej=∑τ∈Thητ2.Then we adopt the following mesh refinement strategy. All the trianglesτ∈𝒯hj satisfying ητ2≥(αEj/n) are divided into four new triangles in 𝒯hj+1 by joining the midpoints of the edges, where n is the numbers of the elements of 𝒯hj, α is a given constant. In order to maintain the new triangulation 𝒯hj+1 to be regular and conformal, some additional triangles need to be divided into two or four new triangles depending on whether they have one or more neighbor which have refined. Then we obtain the new mesh 𝒯hj+1. The above procedure will continue until Ej≤tol, where tol is a given tolerance error.


## 5. Numerical Example

The purpose of this section is to illustrate our theoretical results by numerical example. Our numerical example is the following optimal control problem:(5.1)min⁡u∈K⊂U{12∫0T(‖p-pd‖2+‖y-yd‖2+‖u-u0‖2)dt},(5.2)yt+div⁡p=f+u,p=-∇y,x∈Ω,y|∂Ω=0,y(x,0)=0,(5.3)-zt+div⁡q=y-yd,q=-(∇z+p-pd),x∈Ω,z|∂Ω=0,z(x,T)=0.In our example, we choose the domainΩ=[0,1]×[0,1]. Let Ω be partitioned into 𝒯h as described Section 2. For the constrained optimization problem,(5.4)min⁡u∈KF(u),
where F(u) is a convex functional on U and K={u∈L2(Ω):u≥0a.e.inΩ×J}, the iterative scheme reads (n=0,1,2,…)(5.5)b(un+(1/2),v)=b(un,v)-ρn(F′(un),v),∀v∈U,un+1=PKb(un+(1/2)),
where b(·,·) is a symmetric and positive definite bilinear form such that there exist constants c0 and c1 satisfying(5.6)|b(u,v)|≤c1‖u‖U‖v‖U,∀u,v∈U,b(u,u)≥c0‖u‖U2,
and the projection operator PKbU→K is defined. For given w∈U find PKbw∈K such that(5.7)b(PKbw-w,PKbw-w)=min⁡u∈Kb(u-w,u-w).
The bilinear form b(·,·) provides suitable preconditioning for the projection algorithm. An application of (5.5) to the discretized parabolic optimal control problem yields the following algorithm: (5.8)b(un+(1/2),vh)=b(un,vh)-ρn∫0T(un+zn,vh),∀vh∈Uh,∫0T((pn,vh)-(yn,div⁡vh))=0,∀vh∈Vh,∫0T((ynt,wh)+(div⁡pn,wh))+(yn(0),w(0))=∫0T(f+un,wh),∀wh∈Wh,∫0T((qn,vh)-(zn,div⁡vh))=-∫0T(pn-pd,vh),∀vh∈Vh,∫0T(-(znt,wh)+(div⁡qn,wh))+(zn(T),wh(T))=∫0T(yn-yd,wh),∀wh∈Wh,un+1=PKb(un+(1/2)),un+(1/2),un∈Uh.
The main computational effort is to solve the four state and costate equations and to compute the projection PKbun+(1/2). In this paper we use a fast algebraic multigrid solver to solve the state and costate equations. Then it is clear that the key to saving computing time is how to compute PKbun+(1/2) efficiently. If one uses the C0 finite elements to approximate to the control, then one has to solve a global variational inequality, via, for example, semismooth Newton method. The computational load is not trivial. For the piecewise constant elements, Kh={uh:uh≥0} and b(u,v)=(u,v)U, then
(5.9)PKbun+(1/2)|τ=max⁡(0,avg(un+(1/2))|τ),
where avg(un+(1/2))|τ is the average of un+(1/2) over τ.In solving our discretized optimal control problem, we use the preconditioned projection gradient method withb(u,v)=(u,v)U and a fixed step size ρ=0.8. We now briefly describe the solution algorithm to be used for solving the numerical examples in this section.


### 5.1. Algorithm

(1)
Solve the discretized optimization problem with the projection gradient method on the current meshes, and calculate the error estimatorsηi.(2)
Adjust the meshes using the estimators, and update the solution on new meshes, as described.Now, we present a numerical example to illustrate our theoretical results.Example 5.1.
We choose the state function by(5.10)y(x1,x2)=sinπx1sinπx2sinπt
and the function f(x1,x2)=yt+divp-u with
(5.11)p(x1,x2)=-(πcos⁡πx1sinπx2sinπt,πsinπx1cos⁡πx2sinπt),q(x1,x2)=pd(x1,x2)=p(x1,x2).The costate function can be chosen as(5.12)z(x1,x2)=sinπx1sinπx2sinπt.
It follows from (5.2)-(5.3) that(5.13)yd(x1,x2)=y+zt-div⁡q.We assume that(5.14)λ={0.5,x1+x2>1.0,0.0,x1+x2≤1.0,u0(x1,x2)=1-sinπx12-sinπx22+λ.
Thus, the control function is given by(5.15)u(x1,x2)=max⁡(u0-z,0).In this example, the optimal control has a strong discontinuity, introduced byu0. The exact solution for the control u is plotted in Figure 1. The control function u is discretized by piecewise constant functions, whereas the state (y,p) and the costate (z,q) were approximation by the lowest-order Raviart-Thomas mixed finite elements. In Table 1, numerical results of u, y, and z on uniform and adaptive meshes are presented. It can be founded that the adaptive meshes generated using our error indicators can save substantial computational work, in comparison with the uniform meshes. At the same time, for the discontinuous control variable u, the accuracy has been improved obviously from the uniform meshes to the adaptive meshes in Table 1.Table 1
Numerical results on uniform and adaptive meshes.
On uniform meshOn adaptive meshuyzuyzNodes809780978097108913931393Sides239682396823968282538193819Elements158721587215872634824232423Dofs158721587215872634824232423TotalL2 error6.915e-034.065e-34.018e-36.527e-034.346e-34.323e-3Figure 1
The profile of the control solutionu at t=0.25.In Figure2, the adaptive mesh for u at t=0.25 is shown. In the computing, we use η1-η3 as the error indicators in the adaptive finite element method. It can be founded that the mesh adapts well to be neighborhood of the discontinuity, and a higher density of node points is indeed distributed along them.Figure 2
The adaptive meshes of the control solutionu at t=0.25.


## 5.1. Algorithm

(1)
Solve the discretized optimization problem with the projection gradient method on the current meshes, and calculate the error estimatorsηi.(2)
Adjust the meshes using the estimators, and update the solution on new meshes, as described.Now, we present a numerical example to illustrate our theoretical results.Example 5.1.
We choose the state function by(5.10)y(x1,x2)=sinπx1sinπx2sinπt
and the function f(x1,x2)=yt+divp-u with
(5.11)p(x1,x2)=-(πcos⁡πx1sinπx2sinπt,πsinπx1cos⁡πx2sinπt),q(x1,x2)=pd(x1,x2)=p(x1,x2).The costate function can be chosen as(5.12)z(x1,x2)=sinπx1sinπx2sinπt.
It follows from (5.2)-(5.3) that(5.13)yd(x1,x2)=y+zt-div⁡q.We assume that(5.14)λ={0.5,x1+x2>1.0,0.0,x1+x2≤1.0,u0(x1,x2)=1-sinπx12-sinπx22+λ.
Thus, the control function is given by(5.15)u(x1,x2)=max⁡(u0-z,0).In this example, the optimal control has a strong discontinuity, introduced byu0. The exact solution for the control u is plotted in Figure 1. The control function u is discretized by piecewise constant functions, whereas the state (y,p) and the costate (z,q) were approximation by the lowest-order Raviart-Thomas mixed finite elements. In Table 1, numerical results of u, y, and z on uniform and adaptive meshes are presented. It can be founded that the adaptive meshes generated using our error indicators can save substantial computational work, in comparison with the uniform meshes. At the same time, for the discontinuous control variable u, the accuracy has been improved obviously from the uniform meshes to the adaptive meshes in Table 1.Table 1
Numerical results on uniform and adaptive meshes.
On uniform meshOn adaptive meshuyzuyzNodes809780978097108913931393Sides239682396823968282538193819Elements158721587215872634824232423Dofs158721587215872634824232423TotalL2 error6.915e-034.065e-34.018e-36.527e-034.346e-34.323e-3Figure 1
The profile of the control solutionu at t=0.25.In Figure2, the adaptive mesh for u at t=0.25 is shown. In the computing, we use η1-η3 as the error indicators in the adaptive finite element method. It can be founded that the mesh adapts well to be neighborhood of the discontinuity, and a higher density of node points is indeed distributed along them.Figure 2
The adaptive meshes of the control solutionu at t=0.25.



---

*Source: 217493-2011-08-28.xml*