Adaptive Mixed Finite Element Methods for Parabolic Optimal Control Problems

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.

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 Section 3. 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 notation 𝑊^𝑚,𝑝(Ω) for Sobolev spaces on Ω with a norm ‖⋅‖_𝑚,𝑝 given by ‖𝑣‖^𝑝𝑚,_p=∑_{|𝛼|≤𝑚}‖𝐷^𝛼𝑣‖^𝑝𝐿_𝑝(Ω), a seminorm |⋅|_𝑚,𝑝 given by |𝑣|^{𝑝𝑚,𝑝}=∑_|𝛼|=𝑚‖𝐷^𝛼𝑣‖^𝑝𝐿_𝑝(Ω). We set 𝑊^0𝑚,𝑝(Ω)={𝑣∈𝑊^𝑚,𝑝(Ω)∶𝑣|_𝜕Ω=0}. For 𝑝=2, we denote 𝐻^𝑚(Ω)=𝑊^𝑚,2(Ω),𝐻^𝑚0(Ω)=𝑊^0𝑚,2(Ω), and ‖⋅‖_𝑚=‖⋅‖_𝑚,2,‖⋅‖=‖⋅‖_0,2.

We denote by 𝐿^𝑠(0,𝑇;𝑊^𝑚,𝑝(Ω)) the Banach space of all 𝐿^𝑠 integrable functions from 𝐽 into 𝑊^𝑚,𝑝(Ω) with norm ‖𝑣‖_𝐿^𝑠_(𝐽;𝑊^𝑚,𝑝_(Ω))=(∫^𝑇0‖𝑣‖^𝑠𝑊_{𝑚,𝑝(Ω)}𝑑𝑡)^1/𝑠for𝑠∈[1,∞) and the standard modification for 𝑠=∞. The details can be found in [25].

The parabolic optimal control problems that we are interested in are as follows:min_{𝑢∈𝐾⊂𝑈}^𝑇0𝑔₁𝑦(𝐩)+𝑔₂(𝑦)+𝑗(𝑢)𝑑𝑡_𝑡,𝑦(𝑥,𝑡)=0(𝑥,𝐩,𝑡(𝑥)𝑥+,∈𝑡d)𝜕i=vΩ𝐩,−(𝑡𝐴𝑥∈,(𝑡𝑥)𝐽),=∇𝑦𝑓(+𝑥𝑦,𝑢(𝑡𝑥)(,𝑥,0,)𝑡)𝑥=,∈𝑦₀𝑥Ω,∈Ω,(𝑥),𝑥∈Ω,(2.1) where the bounded open set Ω⊂ℝ² is a convex polygon with the boundary 𝜕Ω. Let 𝐾 be a closed convex set in 𝑈=𝐿²(𝐽;𝐿²(Ω)), 𝑓∈𝐿²(𝐽;𝐿²(Ω)), 𝐽=[0,𝑇], 𝑦₀(𝑥)∈𝐻¹⁰(Ω). Furthermore, we assume that the coefficient matrix 𝐴(𝑥)=(𝑎_𝑖,𝑗(𝑥))_2×2∈𝐿^∞(Ω;ℝ^2×2) is a symmetric 2×2-matrix and there is a constant 𝑐>0 satisfying for any vector 𝐗∈ℝ², 𝐗^𝑡𝐴𝐗≥𝑐‖𝐗‖^2ℝ₂. 𝑗^ is positive, 𝑔^1, 𝑔^2, and 𝑗^ are locally Lipschitz on 𝐿²(Ω)², 𝑊, 𝑈, and that there is a 𝑐>0 such that (𝑗^(𝑢)−𝑗^(𝑢)̃,𝑢−𝑢)̃≥𝑐‖𝑢−𝑢‖̃₀, forall𝑢,𝑢̃∈𝑈.

Now we will describe the mixed finite element discretization of parabolic optimal control problems (2.1). Let 𝐕=𝐻(div;Ω)=𝐯∈𝐿²(Ω)²,div𝐯∈𝐿²(Ω),𝑊=𝐿²(Ω).(2.2) The Hilbert space 𝐕 is equipped with the following norm: ‖𝐯‖_div=‖𝐯‖_𝐻(div;Ω)=‖𝐯‖^20,Ω+‖div𝐯‖^20,Ω^1/2.(2.3)

We recast (2.1) as the following weak form: find (𝐩,𝑦,𝑢)∈𝐕×𝑊×𝐾 such thatmin_{𝑢∈𝐾⊂𝑈}^𝑇0𝑔₁𝐴(𝐩)+𝑔₂(𝑦)+𝑗(𝑢)𝑑𝑡,(2.4)₋₁𝑦𝐩,𝐯−(𝑦,div𝐯)=0,∀𝐯∈𝐕,(2.5)_𝑡,𝑤+𝑦(d(𝑥iv,0𝐩),𝑤=)𝑦=₀(𝑓+𝑢,𝑤),∀𝑤∈𝑊,(2.6)(𝑥),∀𝑥∈Ω.(2.7)

Similar to [26], the optimal control problems (2.4)–(2.7) have a unique solution (𝐩,𝑦,𝑢), and a triplet (𝐩,𝑦,𝑢) is the solution of (2.4)–(2.7) if and only if there is a costate (𝐪,𝑧)∈𝐕×𝑊 such that (𝐩,𝑦,𝐪,𝑧,𝑢) satisfies the following optimality conditions:𝑦𝐴⁻¹𝐩,𝐯−(𝑦,div𝐯)=0,∀𝐯∈𝐕,(2.8)_𝑡,𝑤+𝑦((d𝑥iv,0𝐩),=𝑤)𝑦₀=(𝑓+𝑢,𝑤),∀𝑤∈𝑊,(2.9)𝐴(𝑥),∀𝑥∈Ω,(2.10)₋₁−𝑧𝐪,𝐯−(𝑧,div𝐯)=−𝑔^1(𝐩),𝐯,∀𝐯∈𝐕,(2.11)_𝑡,𝑤+(𝑧di(v𝑥𝐪,𝑇,𝑤))==0,𝑔^2∀(𝑥𝑦)∈,𝑤Ω,,∀𝑤∈𝑊,((22..1123))_𝑇0𝑗^(𝑢)+𝑧,𝑢̃−𝑢_𝑈𝑑𝑡≥0,∀𝑢̃∈𝐾,(2.14) where (⋅,⋅)_𝑈 is the inner product of 𝑈. In the rest of the paper, we will simply write the product as (⋅,⋅) whenever no confusion should be caused.

Let 𝒯_ℎ be regular triangulation of Ω. They are assumed to satisfy the angle condition which means that there is a positive constant 𝐶 such that, for all 𝜏∈𝒯_ℎ, 𝐶⁻¹ℎ^2𝜏≤|𝜏|≤𝐶ℎ^2𝜏, where |𝜏| is the area of 𝜏, ℎ_𝜏 is the diameter of 𝜏 and ℎ=maxℎ_𝜏. In addition 𝐶 or 𝑐 denotes a general positive constant independent of ℎ.

Let 𝐕_ℎ×𝑊_ℎ⊂𝐕×𝑊 denote the Raviart-Thomas space [27] of the lowest order associated with the triangulation 𝒯_ℎ of Ω. 𝑃_𝑘 denotes the space of polynomials of total degree at most 𝑘. Let 𝐕(𝜏)={𝐯∈𝑃²⁰(𝜏)+𝑥⋅𝑃₀(𝜏)}, 𝑊(𝜏)=𝑃₀(𝜏). We define 𝑊𝐕_ℎ∶=𝐯_ℎ∈𝐕∶∀𝜏∈𝒯_ℎ,𝐯_ℎ|_𝜏∈𝐕(𝜏)_ℎ,𝐾∶=𝑤_ℎ∈𝑊∶∀𝜏∈𝒯_ℎ,𝑤_ℎ|_𝜏∈𝑊(𝜏)_ℎ,∶=𝑢_ℎ̃∈𝐾∶∀𝜏∈𝒯_ℎ,𝑢_ℎ̃|_𝜏=constant.(2.15)

The mixed finite element discretization of (2.4)–(2.7) is as follows: compute (𝐩_ℎ,𝑦_ℎ,𝑢_ℎ)∈𝐕_ℎ×𝑊_ℎ×𝐾_ℎ such thatmin_{𝑢ℎ∈𝐾ℎ}𝐴^𝑇0𝑔₁𝐩_ℎ+𝑔₂𝑦_ℎ+𝑗𝑢_ℎ𝑑𝑡₋₁,𝑦_ℎ𝑡𝐩_ℎ,𝐯_ℎ−𝑦_ℎ,div𝐯_ℎ=0,∀𝐯_ℎ∈𝐕_ℎ,,𝑤_ℎ+d𝑦i_ℎv𝐩_ℎ,𝑤_ℎ=𝑓+𝑢_ℎ,𝑤_ℎ,∀𝑤_ℎ∈𝑊_ℎ,(𝑥,0)=𝑦^ℎ0(𝑥),∀𝑥∈Ω,(2.16) where 𝑦^ℎ0(𝑥)∈𝑊_ℎ is an approximation of 𝑦₀. The optimal control problem (2.16) again has a unique solution (𝐩_ℎ,𝑦_ℎ,𝑢_ℎ), and a triplet (𝐩_ℎ,𝑦_ℎ,𝑢_ℎ) is the solution of (2.16) if and only if there is a costate (𝐪_ℎ,𝑧_ℎ)∈𝐕_ℎ×𝑊_ℎ such that (𝐩_ℎ,𝑦_ℎ,𝐪_ℎ,𝑧_ℎ,𝑢_ℎ) satisfies the following optimality conditions:𝑦_ℎ𝑡𝐴⁻¹𝐩_ℎ,𝐯_ℎ−𝑦_ℎ,div𝐯_ℎ=0,∀𝐯_ℎ∈𝐕_ℎ,,𝑤_ℎ+d𝑦i_ℎv𝐩_ℎ,𝑤_ℎ=𝑓+𝑢_ℎ,𝑤_ℎ,∀𝑤_ℎ∈𝑊_ℎ,𝐴(𝑥,0)=𝑦₀(𝑥),∀𝑥∈Ω,₋₁−𝑧𝐪_ℎℎ𝑡,𝐯_ℎ−𝑧_ℎ,div𝐯_ℎ=−𝑔^1𝐩_ℎ,𝐯_ℎ,∀𝐯_ℎ∈𝐕_ℎ,,𝑤_ℎ+di𝑧v_ℎ𝐪_ℎ,𝑤_ℎ=𝑔^2𝑦_ℎ,𝑤_ℎ,∀𝑤_ℎ∈𝑊_ℎ,𝑗(𝑥,𝑇)=0,∀𝑥∈Ω,_𝑢_ℎ+𝑧_ℎ,𝑢_ℎ̃−𝑢_ℎ≥0,∀𝑢_ℎ̃∈𝐾_ℎ.(2.17)

We now consider the fully discrete approximation for the semidiscrete problem. Let Δ𝑡>0, 𝑁=𝑇/Δ𝑡∈ℤ, and 𝑡_𝑖=𝑖Δ𝑡, 𝑖∈ℤ. Also, let 𝜓^𝑖=𝜓^𝑖(𝑥)=𝜓𝑥,𝑡_𝑖,𝑑_𝑡𝜓^𝑖=𝜓^𝑖−𝜓^𝑖−1Δ𝑡.(2.18) For 𝑖=1,2,…,𝑁, construct the finite element spaces 𝐕^𝑖ℎ∈𝐕, 𝑊^𝑖ℎ∈𝑊 (similar as 𝐕_ℎ). Similarly, construct the finite element spaces 𝐾^𝑖ℎ∈𝐾_ℎ with the mesh 𝒯^𝑖ℎ. Let ℎ^𝑖𝜏 denote the maximum diameter of the element 𝜏^𝑖 in 𝒯^𝑖ℎ. Define mesh functions 𝜏(⋅) and mesh size functions ℎ_𝜏(⋅) such that 𝜏(𝑡)|_{𝑡∈(𝑡𝑖−1,𝑡𝑖]}=𝜏^𝑖, ℎ_𝜏(𝑡)|_{𝑡∈(𝑡𝑖−1,𝑡𝑖]}=ℎ_𝜏𝑖. For ease of exposition, we will denote 𝜏(𝑡) and ℎ_𝜏(𝑡) by 𝜏 and ℎ_𝜏, respectively.

The following fully discrete approximation scheme is to find (𝐩^𝑖ℎ,𝑦^𝑖ℎ,𝑢^𝑖ℎ)∈𝐕^𝑖ℎ×𝑊^𝑖ℎ×𝐾^𝑖ℎ,𝑖=1,2,…,𝑁, such thatmin_{𝑢𝑖ℎ∈𝐾𝑖ℎ}^𝑁_𝑖=1^𝑡𝑖_𝑡𝑖−1𝐴𝑔₁𝐩^𝑖ℎ+𝑔₂𝑦^𝑖ℎ+𝑗𝑢^𝑖ℎ,(2.19)₋₁𝑑𝐩^𝑖ℎ,𝐯_ℎ−𝑦^𝑖ℎ,div𝐯_ℎ=0,∀𝐯_ℎ∈𝐕^𝑖ℎ,(2.20)_𝑡𝑦^𝑖ℎ,𝑤_ℎ+div𝑦𝐩^𝑖ℎ,𝑤_ℎ=𝑓𝑥,𝑡_𝑖+𝑢^𝑖ℎ,𝑤_ℎ,∀𝑤_ℎ∈𝑊^𝑖ℎ,(2.21)_0ℎ(𝑥,0)=𝑦^ℎ0(𝑥),∀𝑥∈Ω.(2.22) It follows that the optimal control problems (2.19)–(2.22) have a unique solution (𝐩^𝑖ℎ,𝑦^𝑖ℎ,𝑢^𝑖ℎ),𝑖=1,2,…,𝑁, and that a triplet (𝐩^𝑖ℎ,𝑦^𝑖ℎ,𝑢^𝑖ℎ)∈𝐕^𝑖ℎ×𝑊^𝑖ℎ×𝐾^𝑖ℎ,  𝑖=1,2,…,𝑁, is the solution of (2.19)–(2.22) if and only if there is a costate (𝐪^ℎ𝑖−1,𝑧^ℎ𝑖−1)∈𝐕^𝑖ℎ×𝑊^𝑖ℎ such that (𝐩^𝑖ℎ,𝑦^𝑖ℎ,𝐪^ℎ𝑖−1,𝑧^ℎ𝑖−1,𝑢^𝑖ℎ)∈(𝐕^𝑖ℎ×𝑊^𝑖ℎ)²×𝐾^𝑖ℎ satisfies the following optimality conditions:𝐴𝑑⁻¹𝐩^𝑖ℎ,𝐯_ℎ−𝑦^𝑖ℎ,div𝐯_ℎ=0,∀𝐯_ℎ∈𝐕^𝑖ℎ,(2.23)_𝑡𝑦^𝑖ℎ,𝑤_ℎ+div𝑦𝐩^𝑖ℎ,𝑤_ℎ=𝑓𝑥,𝑡_𝑖+𝑢^𝑖ℎ,𝑤_ℎ,∀𝑤_ℎ∈𝑊^𝑖ℎ,(2.24)_0ℎ𝐴(𝑥,0)=𝑦^ℎ0(𝑥),∀𝑥∈Ω,(2.25)₋₁−𝑑𝐪^ℎ𝑖−1,𝐯_ℎ−𝑧^ℎ𝑖−1,div𝐯_ℎ=−𝑔^1𝐩^𝑖ℎ,𝐯_ℎ,∀𝐯_ℎ∈𝐕^𝑖ℎ,(2.26)_𝑡𝑧^𝑖ℎ,𝑤_ℎ+div𝐪𝑧^ℎ𝑖−1,𝑤_ℎ=𝑔^2𝑦^𝑖ℎ,𝑤_ℎ,∀𝑤_ℎ∈𝑊^𝑖ℎ,(2.27)_𝑁ℎ𝑗(𝑥,𝑇)=0,∀𝑥∈Ω,(2.28)_𝑢^𝑖ℎ+𝑧^ℎ𝑖−1,𝑢_ℎ̃−𝑢^𝑖ℎ≥0,∀𝑢_ℎ̃∈𝐾^𝑖ℎ.(2.29)

For 𝑖=1,2,…,𝑁, let 𝑌_ℎ|_{(𝑡𝑖−1,𝑡𝑖]}=𝑡^𝑖−𝑡𝑦^ℎ𝑖−1+𝑡−𝑡_𝑖−1𝑦^𝑖ℎ𝑍Δ𝑡_ℎ,|_{(𝑡𝑖−1,𝑡𝑖]}=𝑡^𝑖−𝑡𝑧^ℎ𝑖−1+𝑡−𝑡_𝑖−1𝑧^𝑖ℎ𝑃Δ𝑡_ℎ,|_{(𝑡𝑖−1,𝑡𝑖]}=𝑡^𝑖−𝑡𝐩^ℎ𝑖−1+𝑡−𝑡_𝑖−1𝐩^𝑖ℎ𝑄Δ𝑡_ℎ,|_{(𝑡𝑖−1,𝑡𝑖]}=𝑡^𝑖−𝑡𝐪^ℎ𝑖−1+𝑡−𝑡_𝑖−1𝐪^𝑖ℎ𝑈Δ𝑡_ℎ,|_{(𝑡𝑖−1,𝑡𝑖]}=𝑢^𝑖ℎ.(2.30)

For any function 𝑤∈𝐶(𝐽;𝐿²(Ω)), let 𝑤(𝑥,𝑡)|_{𝑡∈(𝑡𝑖−1,𝑡𝑖]}=𝑤𝑥,𝑡_𝑖,𝑤(𝑥,𝑡)|_{𝑡∈(𝑡𝑖−1,𝑡𝑖]}=𝑤𝑥,𝑡_𝑖−1.(2.31)

Then the optimality conditions (2.23)–(2.29) can be restated as.𝐴⁻¹𝑃_ℎ,𝐯_ℎ𝑌−𝑌^ℎ,div𝐯_ℎ=0,∀𝐯^ℎ∈𝐕^𝑖ℎ,(2.32)_ℎ𝑡,𝑤_ℎ+d𝑌iv𝑃_ℎ,𝑤_ℎ=𝑓+𝑈^ℎ,𝑤_ℎ,∀𝑤^ℎ∈𝑊^𝑖ℎ,(2.33)_ℎ𝐴(𝑥,0)=𝑦^ℎ0(𝑥),∀𝑥∈Ω,(2.34)₋₁𝑄_ℎ,𝐯_ℎ−𝑍−𝑍^ℎ,div𝐯_ℎ=−𝑔^1𝑃^ℎ,𝐯^ℎ,∀𝐯^ℎ∈𝐕^𝑖ℎ,(2.35)_ℎ𝑡,𝑤_ℎ+div𝑍𝑄_ℎ,𝑤_ℎ=𝑔^2𝑌^ℎ,𝑤^ℎ,∀𝑤^ℎ∈𝑊^𝑖ℎ,(2.36)_ℎ𝑗(𝑥,𝑇)=0,∀𝑥∈Ω,(2.37)_𝑈_ℎ+𝑍_ℎ,𝑢_ℎ̃−𝑈_ℎ≥0,∀𝑢^ℎ̃∈𝐾^𝑖ℎ.(2.38)

In the rest of the paper, we will use some intermediate variables. For any control function 𝑈_ℎ∈𝐾, we first define the state solution (𝐩(𝑈_ℎ),𝑦(𝑈_ℎ),𝐪(𝑈_ℎ),𝑧(𝑈_ℎ)) which satisfies𝑦𝐴⁻¹𝐩𝑈_ℎ,𝐯−𝑦𝑈_ℎ,div𝐯=0,∀𝐯∈𝐕,(2.39)_𝑡𝑈_ℎ,𝑤𝑦+𝑈div𝐩𝑈_ℎ,𝑤=𝑓+𝑈_ℎ,𝑤,∀𝑤∈𝑊,(2.40)_ℎ𝐴(𝑥,0)=𝑦₀(𝑥),∀𝑥∈Ω,(2.41)₋₁−𝑧𝐪𝑈_ℎ,𝐯−𝑧𝑈_ℎ,div𝐯=−𝑔^1𝐩𝑈_ℎ,𝐯,∀𝐯∈𝐕,(2.42)_𝑡𝑈_ℎ,𝑤+d𝑧iv𝑈𝐪𝑈_ℎ,𝑤=𝑔^2𝑦𝑈_ℎ,𝑤,∀𝑤∈𝑊,(2.43)_ℎ(𝑥,𝑇)=0,∀𝑥∈Ω.(2.44)

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 given 𝑢∈𝐾, let ℳ₁, ℳ₂ be the inverse operators of the state equation (2.6), such that 𝐩(𝑢)=ℳ₁𝑢 and 𝑦(𝑢)=ℳ₂𝑢 are the solutions of the state equations (2.6). Similarly, for given 𝑈_ℎ∈𝐾_ℎ, 𝑃_ℎ(𝑈_ℎ)=ℳ_1ℎ𝑈_ℎ, 𝑌_ℎ(𝑈_ℎ)=ℳ_2ℎ𝑈_ℎ are the solutions of the discrete state equation (2.33). Let𝐹𝐹(𝑢)=𝑔₁ℳ₁𝑈+𝑔₂ℳ₂𝑈+𝑗(𝑢),_ℎ𝑈_ℎ=𝑔₁ℳ_1ℎ𝑈_ℎ+𝑔₂ℳ_2ℎ𝑈_ℎ+𝑗𝑈_ℎ.(3.1) It can be shown that𝐹𝐹^(𝑢),𝑣=𝑗^(𝑢)+𝑧,𝑣_,𝐹𝑈_ℎ,𝑣=𝑗^𝑈_ℎ+𝑧𝑈_ℎ,𝑣_ℎ,𝑈_ℎ,𝑣=𝑗^𝑈_ℎ+𝑍_ℎ,𝑣.(3.2)

It is clear that 𝐹 and 𝐹_ℎ are well defined and continuous on 𝐾 and 𝐾^𝑖ℎ. Also the functional 𝐹_ℎ can be naturally extended on 𝐾. Then (2.4) and (2.19) can be represented asmmini_{𝑈ℎ∈𝐾𝑖ℎ}n_𝑢∈𝐾{𝐹(𝑢)},(3.3)𝐹_ℎ𝑈_ℎ.(3.4)

In many application, 𝐹(⋅) is uniform convex near the solution 𝑢. The convexity of 𝐹(⋅) 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, 𝑢→𝑔₁(ℳ₁𝑈) and 𝑢→𝑔₂(ℳ₂𝑈) are convex. Then there is a constant 𝑐>0, independent of ℎ, such that^𝑇0𝐹^(𝑢)−𝐹^𝑈_ℎ,𝑢−𝑈_ℎ_𝑈≥𝑐‖‖𝑢−𝑈_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}.(3.5)

Theorem 3.1. Let 𝑢 and 𝑈_ℎ be the solutions of (3.3) and (3.4), respectively. Assume that 𝐾^𝑖ℎ⊂𝐾. In addition, assume that (𝐹^ℎ(𝑈_ℎ))|_𝜏∈𝐻^𝑠(𝜏),forall𝜏∈𝒯_ℎ,(𝑠=0,1), and there is a 𝑣_ℎ∈𝐾^𝑖ℎ such that ||𝐹^ℎ𝑈_ℎ,𝑣_ℎ−𝑢||≤𝐶_{𝜏∈𝒯ℎ}ℎ_𝜏‖‖𝐹^ℎ𝑈_ℎ‖‖_𝐻^𝑠_(𝜏)‖‖𝑢−𝑈_ℎ‖‖^𝑠𝐿_2(𝜏).(3.6) Then one has ‖‖𝑢−𝑈_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}≤𝐶𝜂²¹+𝐶‖‖𝑧(𝑈_ℎ)−𝑍_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)},(3.7) where 𝜂²¹=^𝑇0_𝜏_∈𝒯ℎℎ^𝜏1+𝑠‖‖𝑗^(𝑈_ℎ)+𝑍_ℎ‖‖^𝐻1+𝑠_1(𝜏).(3.8)

Proof. It follows from (3.3) and (3.4) that ^𝑇0𝐹^(𝑢),𝑢−𝑣≤0,∀𝑣∈𝐾,_𝑇0𝐹^ℎ𝑈_ℎ,𝑈_ℎ−𝑣_ℎ≤0,∀𝑣_ℎ∈𝐾^𝑖ℎ⊂𝐾.(3.9) Then it follows from assumptions (3.5), (3.6) and Schwartz inequality that 𝑐‖‖𝑢−𝑈_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}≤≤^𝑇0_𝑇0𝐹^(𝑢)−𝐹^𝑈_ℎ,𝑢−𝑈_ℎ≤𝐶_𝑇0𝐹^ℎ𝑈_ℎ,𝑣_ℎ−𝑢+𝐹^ℎ𝑈_ℎ−𝐹^𝑈_ℎ,𝑢−𝑈_ℎ⎧⎪⎨⎪⎩_{𝜏∈𝒯ℎ}ℎ^𝜏1+𝑠‖‖𝐹^ℎ𝑈_ℎ‖‖^{𝐻1+𝑠𝑠}_(𝜏)+‖‖𝐹^ℎ𝑈_ℎ−𝐹^𝑈_ℎ‖‖^2𝐿2_(Ω)+𝑐⎫⎪⎬⎪⎭2‖‖𝑢−𝑈^ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}.(3.10) It is not difficult to show 𝐹^ℎ𝑈_ℎ=𝑗^𝑈_ℎ+𝑍_ℎ,𝐹^𝑈_ℎ=𝑗^𝑈_ℎ+𝑧𝑈_ℎ,(3.11) where 𝑧(𝑈_ℎ) is defined in (2.39)–(2.44). Thanks to (3.11), it is easy to derive ‖‖𝐹^ℎ𝑈_ℎ−𝐹^𝑈_ℎ‖‖_𝐿²_(Ω)≤𝐶‖‖𝑍_ℎ−𝑧𝑈_ℎ‖‖_𝐿²_(Ω).(3.12) 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:−𝜑_𝑡−div(𝐴∇𝜑𝜑|_𝜕Ω)=𝐺,𝑥∈Ω,𝑡∈(0,𝑇],𝜓𝜑(𝑥,𝑇)==00,,𝑡𝑥∈∈𝐽,Ω,(3.13)_𝑡−div𝐴^∗∇𝜓𝜓|_𝜕Ω=𝐺,𝑥∈Ω,𝑡∈(0,𝑇],𝜓(𝑥,0)==00,,𝑡𝑥∈∈𝐽,Ω.(3.14)

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 𝑣=𝜑 or 𝑣=𝜓, ‖‖∇𝑣𝑣‖‖_𝐿^∞_(𝐽;𝐿²_{(Ω))𝐿2(𝐽;𝐿2(Ω))}≤𝐶‖𝐺‖_𝐿²_(𝐽;𝐿²_(Ω)),‖‖𝐷₂≤𝐶‖𝐺‖_𝐿²_(𝐽;𝐿²_(Ω)),‖‖𝑣_𝑡𝑣‖‖_{𝐿2(𝐽;𝐿2(Ω))}≤𝐶‖𝐺‖_𝐿²_(𝐽;𝐿²_(Ω)),‖‖_{𝐿2(𝐽;𝐿2(Ω))}≤𝐶‖𝐺‖_𝐿²_(𝐽;𝐿²_(Ω)),(3.15) where 𝐷²𝑣=𝜕²𝑣/𝜕𝑥_𝑖𝜕𝑥_𝑗,1≤𝑖,𝑗≤2.

Now, we are able to derive the main result.

Theorem 3.3. Let (𝑌_ℎ,𝑃_ℎ,𝑍_ℎ,𝑄_ℎ,𝑈_ℎ) and (𝑦(𝑈_ℎ),𝐩(𝑈_ℎ),𝑧(𝑈_ℎ),𝐪(𝑈_ℎ),𝑈_ℎ) be the solutions of (2.32)–(2.38) and (2.39)–(2.44). Then, ‖‖𝑌_ℎ−𝑦𝑈_ℎ‖‖^2𝐿2_𝐽;𝐿²_(Ω)+‖‖𝑃_ℎ−𝐩𝑈_ℎ‖‖^2𝐿2_𝐽;𝐿²_(Ω)≤𝐶𝜂²²,(3.16) where 𝜂²²=‖‖𝑌_ℎ𝑡+div𝑃_ℎ−𝑓−𝑈_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖𝑌_ℎ−𝑦𝑈_ℎ(𝑥,0)‖‖^2𝐿2_(Ω)+‖‖𝑓−𝑓‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖‖𝑌_ℎ−𝑌_ℎ_𝑡‖‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖𝑃_ℎ−𝑃_ℎ‖‖^2𝐿_{2(𝐽;𝐻(div;Ω))}.(3.17)

Proof. Letting 𝜑 be the solution of (3.13) with 𝐺=𝑌_ℎ−𝑦(𝑈_ℎ), we infer ‖‖𝑌_ℎ−𝑦𝑈_ℎ‖‖^2𝐿2_𝐽;𝐿²_(Ω)==^𝑇0𝑌_ℎ−𝑦𝑈_ℎ,𝐹𝑑𝑡_𝑇0=𝑌_ℎ−𝑦𝑈_ℎ,−𝜑_𝑡−div(𝐴∇𝜑)𝑑𝑡_𝑇0+‖‖𝑌𝑌_ℎ−𝑦𝑈_ℎ_𝑡,𝜑−𝑃_ℎ−𝐩𝑈_ℎ,∇𝜑𝑑𝑡_ℎ−𝑦𝑈_ℎ(𝑥,0)‖‖^2𝐿2_(Ω)=+^𝑇0‖‖𝑌𝑌_ℎ𝑡−𝑦_𝑡𝑈_ℎ,𝜑+div𝑃_ℎ−𝐩𝑈_ℎ,𝜑𝑑𝑡_ℎ−𝑦𝑈_ℎ(𝑥,0)‖‖^2𝐿2_(Ω).(3.18) Then it follows from (2.39)-(2.40) that ‖‖𝑌_ℎ−𝑦𝑈_ℎ‖‖^2𝐿2_𝐽;𝐿²_(Ω)=+^𝑇0‖‖𝑌𝑌^ℎ𝑡,𝜑+div𝑃_ℎ,𝜑−𝑦_𝑡𝑈_ℎ,𝜑−div𝐩𝑈_ℎ,𝜑+div𝑃_ℎ−𝑃_ℎ,𝜑𝑑𝑡_ℎ−𝑦𝑈_ℎ(𝑥,0)‖‖^2𝐿2_(Ω)=+^𝑇0‖‖𝑌𝑌^ℎ𝑡+div𝑃_ℎ−𝑓−𝑈_ℎ,𝜑+𝑓−𝑓,𝜑+div𝑃_ℎ−𝑃_ℎ,𝜑𝑑𝑡_ℎ−𝑦𝑈_ℎ(𝑥,0)‖‖^2𝐿2_(Ω)=+^𝑇0‖‖𝑌𝑌^ℎ𝑡+div𝑃_ℎ−𝑓−𝑈_ℎ,𝜑+𝑓−𝑓,𝜑+div𝑃_ℎ−𝑃_ℎ,𝜑𝑑𝑡_ℎ−𝑦𝑈_ℎ(𝑥,0)‖‖^2𝐿2_(Ω)≤𝐶(𝛿)‖‖𝑌_ℎ𝑡+div𝑃_ℎ−𝑓−𝑈_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+𝐶(𝛿)‖‖𝑓−𝑓‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+𝐶(𝛿)‖‖𝑃_ℎ−𝑃_ℎ‖‖^2𝐿_{2(𝐽;𝐻(div;Ω))}+𝐶(𝛿)‖‖𝑌_ℎ−𝑌_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖𝑌_ℎ−𝑦𝑈_ℎ(𝑥,0)‖‖^2𝐿2_(Ω)+𝐶𝛿‖𝜑‖^2𝐿_{2𝐽;𝐿2(Ω)},(3.19) where and after, 𝛿 is an arbitrary positive number, 𝐶(𝛿) is the constant dependent on 𝛿⁻¹.
Now, we are in the position of estimating the error ‖𝑃_ℎ−𝐩(𝑈_ℎ)‖^2𝐿_{2(𝐽;𝐿2(Ω))}. First, we derive from (2.32)-(2.33) and (2.39)-(2.40) the following useful error equations: 𝑌_ℎ𝐴⁻¹𝑃^ℎ−𝐩𝑈_ℎ,𝐯_ℎ−𝑌^ℎ−𝑦𝑢_ℎ,div𝐯_ℎ=0,(3.20)−𝑦𝑈_ℎ_𝑡,𝑤_ℎ+div𝑃^ℎ−𝐩𝑈_ℎ,𝑤_ℎ=𝑓−𝑓,𝑤^ℎ−𝑌^ℎ−𝑌_ℎ_𝑡,𝑤_ℎ,(3.21) where 𝐯_ℎ∈𝐕_ℎ, 𝑤_ℎ∈𝑊_ℎ. Choose 𝑣_ℎ=𝑃_ℎ−𝐩(𝑈_ℎ) and 𝑤_ℎ=𝑌_ℎ−𝑦(𝑈_ℎ) as the test functions and add the two relations of (3.20)-(3.21) such that 𝐴=⁻¹𝑃𝑓_ℎ−−𝑓𝐩,𝑌𝑈_ℎℎ,𝑃_ℎ−𝐩𝑈_ℎ+𝑌_ℎ−𝑦𝑈_ℎ_𝑡,𝑌_ℎ−𝑦𝑈_ℎ−𝑦𝑈_ℎ−𝑌_ℎ−𝑌_ℎ_𝑡,𝑌_ℎ−𝑦𝑈_ℎ.(3.22) Then, using 𝜖-Cauchy inequality, we can find an estimate as follows: 𝑐‖‖𝑃_ℎ≤−𝐶𝐩‖‖𝑌𝑈_ℎℎ‖‖^2𝐿2_(Ω)+𝑌_ℎ−𝑦𝑈_ℎ_𝑡,𝑌_ℎ−𝑦𝑈_ℎ−𝑦𝑈_ℎ‖‖^2𝐿2_(Ω)+𝐶‖‖𝑓−𝑓‖‖^2𝐿_2(Ω)+𝐶‖‖‖𝑌_ℎ−𝑌_ℎ_𝑡‖‖‖^2𝐿_2(Ω).(3.23) Note that 𝑌^ℎ−𝑦𝑈_ℎ_𝑡,𝑌_ℎ−𝑦𝑈_ℎ=12𝜕𝜕𝑡‖‖𝑌^ℎ−𝑦𝑈_ℎ‖‖^2𝐿2_(Ω),(3.24) then, using the assumption on 𝐴, we verify that 𝑐‖‖𝑃_ℎ−𝐩(𝑈_ℎ)‖‖^2𝐿_2(Ω)+12𝜕𝜕𝑡‖‖𝑌^ℎ−𝑦𝑈_ℎ‖‖^2𝐿2_(Ω)≤𝐶‖‖𝑌_ℎ−𝑦𝑈_ℎ‖‖^2𝐿2_(Ω)+𝐶‖‖𝑓−𝑓‖‖^2𝐿_2(Ω)+𝐶‖‖‖𝑌_ℎ−𝑌_ℎ_𝑡‖‖‖^2𝐿_2(Ω).(3.25) Integrating (3.25) in time and since 𝑌_ℎ(0)−𝑦(𝑈_ℎ)(0)=0, applying Gronwall's lemma, we can easily obtain the following error estimate: ‖‖𝑃_ℎ−𝐩𝑈_ℎ‖‖_𝐿²_(𝐽;𝐿²_(Ω))+‖‖𝑌_ℎ−𝑦𝑈_ℎ‖‖_𝐿^∞_(𝐽;𝐿²_(Ω))≤𝐶‖‖𝑓−𝑓‖‖_{𝐿2(𝐽;𝐿2(Ω))}+𝐶‖‖‖𝑌_ℎ−𝑌_ℎ_𝑡‖‖‖_{𝐿2(𝐽;𝐿2(Ω))}.(3.26) Using the triangle inequality and (3.26), we deduce that ‖‖𝑃_ℎ−𝐩𝑈_ℎ‖‖_𝐿²_(𝐽;𝐿²_(Ω))≤𝐶‖‖𝑓−𝑓‖‖_{𝐿2(𝐽;𝐿2(Ω))}+‖‖‖𝑌_ℎ−𝑌_ℎ_𝑡‖‖‖_{𝐿2(𝐽;𝐿2(Ω))}+‖‖𝑃_ℎ−𝑃_ℎ‖‖_{𝐿2(𝐽;𝐿2(Ω))}.(3.27) Then letting 𝛿 be small enough, it follows from (3.18)–(3.27) that ‖‖𝑌_ℎ−𝑦𝑈_ℎ‖‖^2𝐿2_𝐽;𝐿²_(Ω)+‖‖𝑃_ℎ−𝐩𝑈_ℎ‖‖^2𝐿2_𝐽;𝐿²_(Ω)≤𝐶𝜂²².(3.28) This proves (3.16).

Similarly, letting 𝜓 be the solution of (3.14) with 𝐺=𝑍_ℎ−𝑧(𝑈_ℎ), with the aid of (2.26)-(2.27), (2.42)-(2.43), we can conclude the following.

Theorem 3.4. Let (𝑌_ℎ,𝑃_ℎ,𝑍_ℎ,𝑄_ℎ,𝑈_ℎ) and (𝑦(𝑈_ℎ),𝐩(𝑈_ℎ),𝑧(𝑈_ℎ),𝐪(𝑈_ℎ),𝑈_ℎ) be the solutions of (2.32)–(2.38) and (2.39)–(2.44). Then, ‖‖𝑍_ℎ−𝑧𝑈_ℎ‖‖^2𝐿2_𝐽;𝐿²_(Ω)+‖‖𝑄_ℎ−𝐪𝑈_ℎ‖‖^2𝐿2_𝐽;𝐿²_(Ω)≤𝐶𝜂²²+𝜂²³,(3.29) where 𝜂²³=‖‖𝑔^2𝑌^ℎ+𝑍^ℎ𝑡−div𝑄_ℎ‖‖^2𝐿_{2(𝐽;𝐿2(Ω))}+‖‖𝑍_ℎ−𝑍_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖‖𝑍_ℎ−𝑍_ℎ_𝑡‖‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖𝑌_ℎ−𝑌_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖𝑄_ℎ−𝑄_ℎ‖‖^2𝐿_{2(𝐽;𝐻(div;Ω))}.(3.30)

Let (𝐩,𝑦,𝐪,𝑧,𝑢) and (𝑃_ℎ,𝑌_ℎ,𝑄_ℎ,𝑍_ℎ,𝑈_ℎ) be the solutions of (2.8)–(2.14) and (2.32)–(2.38), respectively. We decompose the errors as follows: 𝐩𝑦−−𝑃𝑌_ℎℎ=𝐩−𝐩𝑈_ℎ+𝐩𝑈_ℎ−𝑃_ℎ∶=𝜖₁+𝜀₁,𝐪−𝑄_ℎ=𝑦−𝑦𝑈_ℎ+𝑦𝑈_ℎ−𝑌_ℎ∶=𝑟₁+𝑒₁,𝑧−𝑍_ℎ=𝐪−𝐪𝑈_ℎ+𝐪𝑈_ℎ−𝑄_ℎ∶=𝜖₂+𝜀₂,=𝑧−𝑧𝑈_ℎ+𝑧𝑈_ℎ−𝑍_ℎ∶=𝑟₂+𝑒₂.(3.31)

From (2.8)–(2.13) and (2.39)–(2.44), we derive the error equations:𝑟𝐴⁻¹𝜖₁,𝐯−𝑟₁,div𝐯=0,∀𝐯∈𝐕,(3.32)_1𝑡𝐴,𝑤+div𝜖₁,𝑤=𝑢−𝑈_ℎ,𝑤,∀𝑤∈𝑊,(3.33)₋₁𝑟𝜖₂,𝐯−𝑟₂,div𝐯=−𝑔^1(𝐩)−𝑔^1𝐩𝑈_ℎ𝐯,∀𝐯∈𝐕,(3.34)_2𝑡,𝑤+div𝜖₂,𝑤=𝑔^2(𝑦)−𝑔^2𝑦𝑈_ℎ,𝑤,∀𝑤∈𝑊.(3.35)

Theorem 3.5. There is a constant 𝐶>0, independent of ℎ, such that ‖‖‖‖𝜖𝜖₁‖‖_{𝐿2(𝐽;𝐿2(Ω))}+‖‖𝑟₁‖‖_{𝐿2(𝐽;𝐿2(Ω))}≤𝐶‖‖𝑢−𝑈_ℎ‖‖_{𝐿2(𝐽;𝐿2(Ω))},(3.36)₂‖‖_{𝐿2(𝐽;𝐿2(Ω))}+‖‖𝑟₂‖‖_{𝐿2(𝐽;𝐿2(Ω))}≤𝐶‖‖𝑢−𝑈_ℎ‖‖_{𝐿2(𝐽;𝐿2(Ω))}.(3.37)

Proof. Part I
Choose 𝐯=𝜖₁ and 𝑤=𝑟₁ as the test functions and add the two relations of (3.32)-(3.33), then we have 𝐴⁻¹𝜖₁,𝜖₁+𝑟_1𝑡,𝑟₁=𝑢−𝑈_ℎ,𝑟₁.(3.38) Then, using the Cauchy inequality, we can find an estimate as follows: 𝐴⁻¹𝜖₁,𝜖₁+𝑟_1𝑡,𝑟₁≤𝐶‖‖𝑟₁‖‖^2𝐿_2(Ω)+‖‖𝑢−𝑈_ℎ‖‖^2𝐿_2(Ω).(3.39) Note that 𝑟_1𝑡,𝑟₁=12𝜕𝜕𝑡‖‖𝑟¹‖‖^2𝐿_2(Ω),(3.40) then, using the assumption on 𝐴, we can obtain that ‖‖𝜖₁‖‖^2𝐿_2(Ω)+12𝜕𝜕𝑡‖‖𝑟¹‖‖^2𝐿_2(Ω)≤𝐶‖‖𝑟₁‖‖^2𝐿_2(Ω)+‖‖𝑢−𝑈_ℎ‖‖^2𝐿_2(Ω).(3.41) Integrating (3.41) in time and since 𝑟₁(0)=0, applying the Gronwall's lemma, we can easily obtain the following error estimate ‖‖𝜖₁‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖𝑟₁‖‖^2𝐿_{2𝐽;𝐿2(Ω)}≤𝐶‖‖𝑢−𝑈_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}.(3.42) This implies (3.36).
Part II
Similarly, choose 𝑣=𝜖₂ and 𝑤=𝑟₂ as the test functions and add the two relations of (3.34)-(3.35), then we can obtain that 𝐴⁻¹𝜖₂,𝜖₂+𝑟_2𝑡,𝑟₂=𝑔^2(𝑦)−𝑔^2𝑦𝑈_ℎ,𝑟₂−𝑔^1(𝐩)−𝑔^1𝐩𝑈_ℎ,𝜖₂.(3.43) Then, using the Cauchy inequality, we can find an estimate as follows: 𝐴⁻¹𝜖₂,𝜖₂+𝑟_2𝑡,𝑟₂≤𝐶‖‖𝑟₁‖‖^2𝐿_2(Ω)+‖‖𝑟₂‖‖^2𝐿_2(Ω)+‖‖𝜖₁‖‖^2𝐿_2(Ω)+2𝑐‖‖𝜖²‖‖^2𝐿_2(Ω).(3.44) Note that 𝑟_2𝑡,𝑟₂=12𝜕𝜕𝑡‖‖𝑟²‖‖^2𝐿_2(Ω),(3.45) then, using the assumption on 𝐴, we verify that ‖‖𝜖₂‖‖^2𝐿_2(Ω)+12𝜕𝜕𝑡‖‖𝑟²‖‖^2𝐿_2(Ω)≤𝐶‖‖𝑟₁‖‖^2𝐿_2(Ω)+‖‖𝑟₂‖‖^2𝐿_2(Ω)+‖‖𝜖₁‖‖^2𝐿_2(Ω).(3.46) Integrating (3.46) in time and since 𝑟₂(𝑇)=0, applying the Gronwall's lemma, we can easily obtain the following error estimate ‖‖𝜖₂‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖𝑟₂‖‖^2𝐿_{2𝐽;𝐿2(Ω)}≤𝐶‖‖𝑢−𝑈_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}.(3.47) Then (3.37) follows from (3.47) and the previous statements immediately.

Collecting Theorems 3.1–3.5, we can derive the following results.

Theorem 3.6. Let (𝐩,𝑦,𝐪,𝑧,𝑢) and (𝑃_ℎ,𝑌_ℎ,𝑄_ℎ,𝑍_ℎ,𝑈_ℎ) be the solutions of (2.8)–(2.14) and (2.32)–(2.38), respectively. In addition, assume that (𝑗^(𝑈_ℎ)+𝑍_ℎ)|_𝜏∈𝐻^𝑠(𝜏),forall𝜏∈𝒯_ℎ,(𝑠=0,1), and that there is a 𝑣_ℎ∈𝐾_ℎ such that |||𝑗_𝑈_ℎ+𝑍_ℎ,𝑣_ℎ−𝑢|||≤𝐶_{𝜏∈𝒯ℎ}ℎ_𝜏‖‖𝑗^𝑈_ℎ+𝑍_ℎ‖‖_{𝐻𝑠(𝜏)}‖‖𝑢−𝑈_ℎ‖‖^𝑠𝐿_2(𝜏).(3.48) Then one has that, forall𝑡∈(0,𝑇], ‖‖𝑢−𝑈_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖𝑦−𝑌_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖𝐩−𝑃_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖𝑧−𝑍_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}+‖‖𝐪−𝑄_ℎ‖‖^2𝐿_{2𝐽;𝐿2(Ω)}≤𝐶³_𝑖=1𝜂^2𝑖,(3.49) where 𝜂₁, 𝜂₂, and 𝜂₃ 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 Section 3 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𝜏_𝑖, where 𝜏_𝑖 is the finite elements. At each iteration, an average quantity of all 𝜂^2𝜏_𝑖 is calculated, and each 𝜂^2𝜏_𝑖 is then compared with this quantity. The element 𝜏_𝑖 is to be refined if 𝜂^2𝜏_𝑖 is larger than this quantity. As 𝜂^2𝜏_𝑖 represents the total approximation error over 𝜏_𝑖, 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 𝒯_ℎ0 of Ω, we construct a sequence of refined triangulation 𝒯_ℎ𝑗 as follows. Given 𝒯_ℎ𝑗, we compute the solutions (𝑃_ℎ,𝑌_ℎ,𝑄_ℎ,𝑍_ℎ,𝑈_ℎ) of the system (2.32)–(2.38) and their error estimator 𝜂^2𝜏=^𝑇0_𝜏_∈𝒯ℎℎ^𝜏1+𝑠‖‖𝑗^𝑈_ℎ+𝑍_ℎ‖‖^𝐻1+𝑠_1(𝜏)+‖‖𝑌_ℎ𝑡+div𝑃_ℎ−𝑓−𝑈_ℎ‖‖^2𝐿_{2𝐽;𝐿2(𝜏)}+‖‖𝑌_ℎ−𝑦𝑈_ℎ(𝑥,0)‖‖^2𝐿2_(𝜏)+‖‖𝑓−𝑓‖‖^2𝐿_{2𝐽;𝐿2(𝜏)}+‖‖‖𝑌_ℎ−𝑌_ℎ_𝑡‖‖‖^2𝐿_{2𝐽;𝐿2(𝜏)}+‖‖𝑃_ℎ−𝑃_ℎ‖‖^2𝐿_{2(𝐽;𝐻(div;𝜏))}+‖‖𝑔^2𝑌^ℎ+𝑍^ℎ𝑡−div𝑄_ℎ‖‖^2𝐿_{2𝐽;𝐿2(𝜏)}+‖‖𝑍_ℎ−𝑍_ℎ‖‖^2𝐿_{2𝐽;𝐿2(𝜏)}+‖‖‖𝑍_ℎ−𝑍_ℎ_𝑡‖‖‖^2𝐿_{2𝐽;𝐿2(𝜏)}𝐸_𝑗+‖‖𝑌_ℎ−𝑌_ℎ‖‖^2𝐿_{2𝐽;𝐿2(𝜏)}+‖‖𝑄_ℎ−𝑄_ℎ‖‖^2𝐿_{2(𝐽;𝐻(div;𝜏))},=_{𝜏∈𝒯ℎ}𝜂^2𝜏.(4.1)

Then we adopt the following mesh refinement strategy. All the triangles 𝜏∈𝒯_ℎ𝑗 satisfying 𝜂^2𝜏≥(𝛼𝐸_𝑗/𝑛) are divided into four new triangles in 𝒯_ℎ𝑗+1 by joining the midpoints of the edges, where 𝑛 is the numbers of the elements of 𝒯_ℎ𝑗, 𝛼 is a given constant. In order to maintain the new triangulation 𝒯_ℎ𝑗+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 𝒯_ℎ𝑗+1. The above procedure will continue until 𝐸_𝑗≤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:min_{𝑢∈𝐾⊂𝑈}12^𝑇0‖‖𝐩−𝐩^𝑑𝑦‖‖²+‖‖𝑦−𝑦_𝑑‖‖²+‖‖𝑢−𝑢₀‖‖²𝑑𝑡,(5.1)_𝑡−𝑧_𝑡+div𝐩=𝑓+𝑢,𝐩=−∇𝑦,𝑥∈Ω,𝑦|_𝜕Ω=0,𝑦(𝑥,0)=0,(5.2)+div𝐪=𝑦−𝑦_𝑑,𝐪=−∇𝑧+𝐩−𝐩_𝑑,𝑥∈Ω,𝑧|_𝜕Ω=0,𝑧(𝑥,𝑇)=0.(5.3)

In our example, we choose the domain Ω=[0,1]×[0,1]. Let Ω be partitioned into 𝒯_ℎ as described Section 2. For the constrained optimization problem,min_𝑢∈𝐾𝐹(𝑢),(5.4) where 𝐹(𝑢) is a convex functional on 𝑈 and 𝐾={𝑢∈𝐿²(Ω)∶𝑢≥0𝑎.𝑒.inΩ×𝐽}, the iterative scheme reads (𝑛=0,1,2,…)𝑏𝑢_𝑛+(1/2),𝑣=𝑏𝑢𝑢_𝑛,𝑣−𝜌_𝑛𝐹^𝑢_𝑛,𝑣,∀𝑣∈𝑈,_𝑛+1=𝑃^𝑏𝐾𝑢_𝑛+(1/2),(5.5) where 𝑏(⋅,⋅) is a symmetric and positive definite bilinear form such that there exist constants 𝑐₀ and 𝑐₁ satisfying||𝑏(𝑢,𝑣)||≤𝑏𝑐(₁𝑢‖,𝑢𝑢‖)_𝑈≥‖𝑐𝑣₀‖_𝑈,∀𝑢,𝑣∈𝑈,‖𝑢‖^2𝑈,(5.6) and the projection operator 𝑃^𝑏𝐾𝑈→𝐾 is defined. For given 𝑤∈𝑈 find 𝑃^𝑏𝐾𝑤∈𝐾 such that𝑏𝑃^𝑏𝐾𝑤−𝑤,𝑃^𝑏𝐾𝑤−𝑤=min_𝑢∈𝐾𝑏(𝑢−𝑤,𝑢−𝑤).(5.7) The bilinear form 𝑏(⋅,⋅) provides suitable preconditioning for the projection algorithm. An application of (5.5) to the discretized parabolic optimal control problem yields the following algorithm: 𝑏𝑢_𝑛+(1/2),𝑣_ℎ=𝑏𝑢_𝑛,𝑣_ℎ−𝜌_𝑛^𝑇0𝑢_𝑛_𝑇0+𝑧_𝑛,𝑣_ℎ,∀𝑣_ℎ∈𝑈_ℎ,_𝑇0𝐩_𝑛,𝐯_ℎ−𝑦_𝑛,div𝐯_ℎ=0,∀𝐯_ℎ∈𝐕_ℎ,𝑦_𝑛𝑡,𝑤_ℎ𝑇0+div𝐩_𝑛,𝑤_ℎ+𝑦_𝑛(0),𝑤(0)=^𝑇0𝑓+𝑢_𝑛,𝑤_ℎ,∀𝑤_ℎ∈𝑊_ℎ,_𝑇0𝐪_𝑛,𝐯_ℎ−𝑧_𝑛,div𝐯_ℎ=−^𝑇0𝐩_𝑛−𝐩_𝑑,𝑣_ℎ,∀𝐯_ℎ∈𝐕_ℎ,−𝑧_𝑛𝑡,𝑤_ℎ+div𝐪𝑢_𝑛,𝑤_ℎ𝑛+1+𝑧_𝑛(𝑇),𝑤_ℎ(𝑇)=^𝑇0𝑦_𝑛−𝑦_𝑑,𝑤_ℎ,∀𝑤_ℎ∈𝑊_ℎ,=𝑃^𝑏𝐾𝑢_𝑛+(1/2),𝑢_𝑛+(1/2),𝑢_𝑛∈𝑈_ℎ.(5.8) The main computational effort is to solve the four state and costate equations and to compute the projection 𝑃^𝑏𝐾𝑢_𝑛+(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 𝑃^𝑏𝐾𝑢_𝑛+(1/2) efficiently. If one uses the 𝐶⁰ 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, 𝐾_ℎ={𝑢_ℎ∶𝑢_ℎ≥0} and 𝑏(𝑢,𝑣)=(𝑢,𝑣)_𝑈, then 𝑃^𝑏𝐾𝑢_𝑛+(1/2)|_𝜏=max0,avg𝑢_𝑛+(1/2)|_𝜏,(5.9) where avg(𝑢_𝑛+(1/2))|_𝜏 is the average of 𝑢_𝑛+(1/2) over 𝜏.

In solving our discretized optimal control problem, we use the preconditioned projection gradient method with 𝑏(𝑢,𝑣)=(𝑢,𝑣)_𝑈 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 𝜂_𝑖.(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 𝑦𝑥₁,𝑥₂=sin𝜋𝑥₁sin𝜋𝑥₂sin𝜋𝑡(5.10) and the function 𝑓(𝑥₁,𝑥₂)=𝑦_𝑡+div𝐩−𝑢 with 𝐩𝑥₁,𝑥₂=−𝐪𝜋𝑥cos𝜋𝑥₁sin𝜋𝑥₂sin𝜋𝑡,𝜋sin𝜋𝑥₁cos𝜋𝑥₂sin𝜋𝑡₁,,𝑥₂=𝐩_𝑑𝑥₁,𝑥₂=𝐩𝑥₁,𝑥₂.(5.11)

The costate function can be chosen as𝑧𝑥₁,𝑥₂=sin𝜋𝑥₁sin𝜋𝑥₂sin𝜋𝑡.(5.12) It follows from (5.2)-(5.3) that𝑦_𝑑𝑥₁,𝑥₂=𝑦+𝑧_𝑡−div𝐪.(5.13)

We assume that𝜆=0.5,𝑥¹0.0,𝑥₁+𝑥₂>1.0,𝑢+𝑥₂≤1.0,₀𝑥₁,𝑥₂=1−sin𝜋2𝑥¹−sin𝜋2𝑥²+𝜆.(5.14) Thus, the control function is given by𝑢𝑥₁,𝑥₂=max𝑢₀−𝑧,0.(5.15)

In this example, the optimal control has a strong discontinuity, introduced by 𝑢₀. The exact solution for the control 𝑢 is plotted in Figure 1. The control function 𝑢 is discretized by piecewise constant functions, whereas the state (𝑦,𝐩) and the costate (𝑧,𝐪) were approximation by the lowest-order Raviart-Thomas mixed finite elements. In Table 1, numerical results of 𝑢, 𝑦, and 𝑧 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 𝑢, 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.

Figure 1 

The profile of the control solution 𝑢 at 𝑡=0.25.

In Figure 2, the adaptive mesh for 𝑢 at 𝑡=0.25 is shown. In the computing, we use 𝜂₁−𝜂₃ 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 solution 𝑢 at 𝑡=0.25.

Acknowledgments

The authors express their thanks to the referees for their helpful suggestions, which led to improvement of the presentation. This work was supported by National Nature Science Foundation under Grant 10971074 and Hunan Provincial Innovation Foundation For Postgraduate CX2009B119.