Immersed finite element method and its analysis for parabolic optimal control problems with interfaces

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.1 (1-22)

Applied Numerical Mathematics ••• (••••) •••–•••

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Immersed finite element method and its analysis for parabolic optimal control problems with interfaces Zhiyue Zhang a , Dong Liang b,∗ , Quanxiang Wang c a b c

School of Mathematical Sciences, Jiangsu Key Laboratory for NSLSCS, Nanjing Normal University, Nanjing 210023, China Department of Mathematics and Statistics, York University, Toronto, Ontario, M3J 1P3, Canada College of Engineering, Nanjing Agricultural University, Nanjing 210031, China

a r t i c l e

i n f o

Article history: Received 4 April 2019 Received in revised form 22 August 2019 Accepted 22 August 2019 Available online xxxx Keywords: PDE-constrained optimization Parabolic interface problems Variational discretization Immersed finite element

a b s t r a c t In this paper, we develop the immersed finite element method for parabolic optimal control problems with interfaces. By employing the definition of directional derivative of Lagrange function, first-order necessary optimality conditions in qualified form for parabolic optimal control problems with interfaces are established. The parabolic state equations are treated with the immersed finite element method using a simple uniform mesh which is independent of the interface. In the interface elements, functions are locally piecewise bilinear functions according to the subelements formed by the actual interface curves. By introducing the auxiliary functions which are the solutions of interface parabolic equations with non-homogeneous and homogeneous jump conditions, optimal error estimates are proved for the proposed schemes to the controls, states and adjoint states in both the semi-discrete case and the fully discrete case. Numerical experiments show the performance of the proposed scheme to solve the parabolic optimal control problems with interfaces and confirm the theoretical results. © 2019 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Optimal control problems governed by parabolic PDEs with interfaces arise in science and engineering, where the optimization or optimal control of a process is in a domain which is composed of several materials separated by curves or surfaces (called interfaces). Coefficients in parabolic PDEs may have a jump across the curve interface corresponding to different materials. Hence, it is challenging to design efficient numerical methods for the optimal control problems that involve parabolic equations with curve interfaces. One difficult task for solving the curve interface PDE problems is to design schemes on unfitted meshes near the interface. The immersed finite element method (IFEM) proposed in [16] is among a few methods based on linear finite element discretization and unfitted meshes, for example, uniform triangulations. The main advantage of the IFEM is that the basis functions in the approximate interface triangles satisfy the interface conditions, while the ones in the non-interface triangles are standard basis functions. Paper [15] gave optimal approximation capabilities of the immersed finite element space. And reference [6] proved optimal error estimates in L 2 and H 1 norms. Recently, [21] gave a posteriori error estimate of the fitted finite element method for linear parabolic interface problems. [14] used noncomforming spectral element method for solving parabolic interface problem with nonhomogeneous jump conditions. In [10] a parabolic problem with a time

*

Corresponding author. E-mail address: [email protected] (D. Liang).

https://doi.org/10.1016/j.apnum.2019.08.024 0168-9274/© 2019 IMACS. Published by Elsevier B.V. All rights reserved.

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.2 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

2

dependent interface was solved by the Crank-Nicolson type immersed finite element methods. A symmetric interior penalty finite element approach was developed to a parabolic interface problem with low regularity solutions in [22]. On optimal control problems with parabolic equations without interfaces, the developments of numerical methods have been made in many publications (see, e.g., [3,5,7,12,17–20,23,24]) and therein references. [7] gave a priori error estimate in L 2 norm for parabolic control problem with low regularity state equation constraints. A posteriori error estimate for both the state and control of the finite element approximation of quadratic optimal control problem governed by linear parabolic equation was obtained in [17]. The state of the art in the optimal a priori analysis of distributed parabolic optimal control problems can be found in [18]. In [19], for the optimal control of non-smooth semilinear parabolic equations, the equivalence of strong stationarity to purely primal conditions was obtained. Numerical arts for the optimal control problem with semi-linear parabolic equations constraints can be found in [3,5,20]. However, to the best of our knowledge, there are few works that concern numerical methods based on unfitted meshes for the parabolic optimal control problems with interfaces. The difficulty in analysis and computation is to discrete the problems near the interfaces where it only has the lower regularity of the state variable for such problems. Constrained by parabolic equation with interface, the solution of the problem has a lower regularity near the interface, which leads to numerical analysis and computation of the optimal control problem much challenging. In this paper, we propose the distributed control problem governed by parabolic interface problems by combining the variational discretization concept with the IFEM on uniform triangulations. By employing the definition of directional derivative of Lagrange function, first-order necessary optimality conditions in qualified form for parabolic optimal control problem with interface are established. In each interface element, functions are locally piecewise bilinear form according to the subelements formed by the actual curve interface replace its line approximation. By introducing the auxiliary functions which are the solution of interface parabolic equations with non-homogeneous and homogeneous jump conditions, optimal error estimates are proved for the proposed schemes to the control, state and adjoint state in semi-discrete case and fully discrete case. Our approach results in the scheme to solve the state y equation and the adjoint state p equation marching in opposite directions. A numerically implementable iteration can be used to solve single equations for the discrete controls. Numerical experiments show the performance of the proposed scheme to solve the parabolic optimal control problems with interfaces and confirm the theoretical analysis results. The remainder of this paper is organized as follows. In section 2, the model problem is introduced and optimality conditions and regularity results of the problem are analyzed. Section 3 presents the discretization of the optimal control problem based on the variational discretization approach and the immersed finite element method. And optimal error estimates are proved. In Section 4, the implementation details of the numerical method are described briefly. In Section 5, numerical experiments are provided. Some conclusions are addressed in Section 6. 2. Parabolic optimal control problem Consider the following parabolic interface problem,

∂t y (x, t ) − ∇ · (β(x)∇ y (x, t )) = f (x, t ) + u (x, t ), [ y] = 0, [β∂n y] = g (x, t ), y (x, t ) = ϕ (x, t ), y (x, 0) = y 0 ,

in (\) × (0, T ),

(2.1)

across  × (0, T ),

(2.2)

on ∂  × (0, T ),

(2.3)

x ∈ ,

(2.4)

where  ∈ R2 is a bounded domain separated by a closed interface  ∈ C 2 and  T = (0, T ] × , [ v ] = v + − v − denotes the jump of the function v (x) across an interface . We assume that the interface  separates the domain  into two sub-domains + and − , and − lies strictly inside , see Fig. 2 for an illustration. The vector n is the unit normal direction of  pointing to + . f ∈ H 1 (0, T ; L 2 ()), g ∈ L 2 (0, T ; H 1/2 ()), ϕ ∈ L 2 (0, T ; H 1/2 (∂ )) and y 0 ∈ H 1 (). The coefficient β(x) is a positive and piecewise constant, that is,

β(x) = β + , if x ∈ + ,

β(x) = β − , if x ∈ − .

(2.5)

For a time interval I = (0, T ) we use the following notation for the inner products and norms on L (), L (0, T ; L ()) and L 2 (0, T ; L 2 (s ))(s = ±): 2

( v , w ) := ( v , w )L 2 () ,  v  :=  v L 2 () ,

( v , w ) I := ( v , w ) L 2 ( I , L 2 ()) ,

 v  I :=  v L 2 ( I , L 2 ()) ,

2

2

( v , w ) I s := ( v , w ) L 2 ( I , L 2 (s ))

 v  I s :=  v L 2 ( I , L 2 (s )) .

The weak formulation of the state equation (2.1)-(2.4) is: find y ∈ L 2 (0, T ; H 1 ()) ∩ H 1 (0, T ; L 2 ()) ∩ { y |∂  = ϕ }) such that

∀ v ∈ L 2 (0, T ; H 01 ()), with y (0) = y 0 ,

(∂t y , v ) + a( y , v ) = ( f + u , v ) −  g , v  , where a( y , v ) =





s=± s

β ∇ y · ∇ vdx = s



s=± (β

s

∇ y , ∇ v )s ,  g , v  =



 g vd.

(2.6)

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.3 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

3

Problem 2.1. (P) Consider the optimal control problem of minimizing

J ( y, u) =

T

1

2

 y − yd  dt +

2

α

T u 2 dt

2

0

(2.7)

0

over all ( y , u ) ∈ L 2 (0, T ; H 1 ()) ∩ H 1 (0, T ; L 2 ()) ∩ { y |∂  = ϕ }) × L 2 (0, T ; L 2 ()) subject to the parabolic interface problem (2.6) and the control constraints

ua ≤ u ≤ ub ,

(2.8)

α is a fixed positive number and the set of admissible controls for (P) can be written as   U ad = u ∈ L 2 ( T ) : ua ≤ u ≤ u b .

The regularization parameter

We make the following smoothness assumption on the data of the problem, that is, yd ∈ L 2 (0, T ; L 2 ()), and ua , u b ∈ ˜ 2 ()) ∩ H 1 (0, T ; H˜ 1 ()) ∩ { y |∂  = ϕ }, H˜ 2 () = y ∈ H 1 () : y ∈ H˜ 2 (s ), s = L 2 (0, T ; L 2 ()). Denote H 2,1 ( T ) = L 2 (0, T ; H

   ± and H˜ 1 () = y ∈ L 2 () : y ∈ H˜ 1 (s ), s = ± equipped with the norms  y 2˜ 2  y 2˜ 1 H ()

=

 y 2H 1 (+ )

+  y 2H 1 (− )

respectively.

2, 1 H 0 ( T )

H () ˜1

=  y 2H 2 (+ ) +  y 2H 2 (− ) and

= L (0, T ; H˜ ()) ∩ H (0, T ; H ()) ∩ { y |∂  = 0}. 2

2

1

Since the problem is quadratic and convex, by applying an idea of the formal Lagrange technique see [12,23], we have the existence of solutions and the optimality conditions. Theorem 2.2. The problem (P) admits a unique optimal control u ∈ L 2 (0, T ; L 2 ()), with an associated state y ∈ L 2 (0, T ; H 1 ()) ∩ H 1 (0, T ; L 2 ()) ∩ { y |∂  = ϕ } and an adjoint state p ∈ L 2 (0, T ; H 01 ()) ∩ H 1 (0, T ; L 2 ()) that satisfy the state equation

(∂t y , v ) + a( y , v ) = ( f + u , v ) −  g , v  ,

∀ v ∈ H 01 () t ∈ (0, T ], with y (x, 0) = y 0 (x),

(2.9)

the adjoint equation

−(∂t p , v ) + a( v , p ) = ( y − yd , v ),

∀ v ∈ H 01 (),

t ∈ (0, T ], with p (x, T ) = 0,

(2.10)

and the variational inequality

(α u + p , w − u ) I ≥ 0,

∀ w ∈ U ad .

(2.11)

Moreover, the variational inequality is equivalent to

u = P[ u a , u b ] −

1

α

p ,

(2.12)

where P[ua ,ub ] ( v ) denotes the projection of v ∈ R onto the interval [ua , u b ]. The adjoint equation (2.10) is the weak form of the following interface problem

−∂t p (x, t ) − ∇ · (β(x)∇ p (x, t )) = y (x, t ) − yd (x, t ), [p] = 0, [β∂n p] = 0,

across  × (0, T ),

in (\) × (0, T ),

(2.13) (2.14)

p (x, t ) = 0,

on ∂  × (0, T ),

(2.15)

p (x, T ) = 0,

x ∈ ,

(2.16)

Proof. We define the Lagrangian functional for problem (P) 2,1

L ( y , u , p ) : H 2,1 ( T ) × U ad × H 0 ( T ) → R, as

(2.17)

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.4 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

4

T 

T 

L( y, u, p) = J ( y, u) −

∂t yp 1 dxdt − 0 \

T



0 \

T

+

β∇ y · ∇ p 1 dxdt



( f + u ) p 1 dxdt − 0 \

gp + 1 dsdt

(2.18)

0 

T 

T 

−

[ y ] p 2 dsdt − 0 

 ( y − ϕ ) p 3 dsdt −

0 ∂

( y (0) − y 0 ) p 1 (0)dx. 

In the following, we first clarify that the equation D y L ( y , u , p )h = 0, for all h leads to the adjoint equation. Using the definition of the directional derivative of (2.18) with respect to the state variable, set equal to zero and evaluated at the optimal pair. Upon integrating by parts with respect to the time t, invoking the Green’s formula, note that [ p 1 ] = 0, p 1 |∂  = 0 and finally regrouping the terms, we obtain

T  D y L ( y , u , p )h =

T +

⎛ ⎝

( y − yd )hdxdt −

+ ∇ · (β + ∇ p + 1 )h dx +

0

+

0

−

h(x, T ) p 1 (x, T )dx +

\

0 



T 





∂t p 1 · hdxdt 0 \

∂ p+ 1 + β+ − h ds − ∂n





⎞ + ∂ p β + 1 h+ ds⎠ dt ∂n

∂

⎛ ⎞ T   ∂ p− − − − − − 1 + ⎝ ∇ · (β ∇ p 1 )h dx − β h ds⎠ dt ∂ n− 

T  −

T  [h] p 2 dsdt −

0 

hp 3 dsdt 0 ∂

T  =

T 

 ( y − yd )hdxdt −

h(x, T ) p 1 (x, T )dx +

\

0 

T

(2.19)

0 \



+

T

 

∂ p1 β ∂n

∇ · (β∇ p 1 )hdxdt + 0 \

0 

T   +

β

+ + ∂ p1

∂n

∂t p 1 · hdxdt





T 

+

h dsdt + 

[h]

∂ p− 1 β− − ∂n

 − p 2 dsdt

0 

 − p 3 hdsdt .

0 ∂

First, we note that for all h ∈ C 0∞ ((\) T ) with homogeneous jump conditions [h] = 0 on the interface  the expressions h( T ) and h, vanish on  and ∂  T ,  T , respectively. We have

T 

T  ( y − yd )hdxdt +

T  ∂t p 1 · hdxdt +

0 \

0 

∇ · (β∇ p 1 )hdxdt = 0,

∀h ∈ C 0∞ ((\)T ).

(2.20)

0 \

Now since that C 0∞ ((\) T ) is dense in L 2 ((\) T ). We therefore must have

−∂t p 1 (x, t ) − ∇ · (β(x)∇ p 1 (x, t )) = y (x, t ) − yd (x, t ),

in (\) × (0, T ).

(2.21)

In particular, the integral over  \  × (0, T ) vanishes in the above equation (2.21). Next, we no longer require that h on ∂  T , and consider the set of all functions h ∈ C ((\) T ) with homogeneous jump condition on interface , note that 2, 1 p 1 ∈ H 0 ( T ), we have [ p 1 ] = 0 and it follows that



[ p 1 ] = 0,

∂ p1 and β ∂n



= 0, 

across  × (0, T ).

(2.22)

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.5 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

5

Finally, we no longer require h( T ) = 0, similarly we conclude that

p 1 (x, T ) = 0, Since p 1 ∈

2, 1 H 0 ( T ),

p 1 (x, t ) = 0,

x ∈ .

(2.23)

we obtain

on ∂  × (0, T ).

(2.24)

In summary, by putting p = p 1 our argument yields that the system as the adjoint equation as (2.13)-(2.16). In addition, if we no long require that [h] = 0 on the  × (0, T ), and h = 0 on the ∂  × (0, T ), respectively, we find that

p2 = β −

∂ p− 1 , ∂ n−

on  × (0, T ),

and p 3 = β +

∂ p+ 1 ∂n

on ∂  × (0, T ).

(2.25)

Secondly, we claim that the equation D p L ( y , u , p )h = 0, for all h leads to the state equation. Using the definition of the directional derivative of with respect to the adjoint variable, set equal to zero and evaluated at the optimal pair and note that (2.25), we easily get

T  D p L ( y , u , p )h = −

(∂t y − ∇ · (β∇ y ) − ( f + u )) hdxdt 0 \



T   T  ∂y ∂ h− + β + − g h dsdt − [ y ] β − − dsdt ∂n  ∂n 0 

T  −

(2.26)

0 

( y − ϕ )β +

∂ h+ ∂n

 dsdt −

0 ∂

( y (x, 0) − y 0 )h(0)dx. 

For all h ∈ C 0∞ ((\) T ) the expressions h(0), and h, ∂∂ nh vanish on  and ∂  T ,  T , respectively. We obtain

T  (∂t y − ∇ · (β∇ y ) − ( f + u )) hdxdt = 0,

∀h ∈ C 0∞ ((\)T ).

(2.27)

0 \

Now observe that C 0∞ ((\) T ) is dense in L 2 ((\) T ), we hence must have

∂t y (x, t ) − ∇ · (β(x)∇ y (x, t )) = f (x, t ) + u (x, t ), − We no longer require ∂∂ nh−

= 0 and h−

[ y] = 0, [β∂n y] = g (x, t ),

in (\) × (0, T ).

(2.28)

= 0 on  × (0, T ), respectively, we conclude that across  × (0, T ).

(2.29)

∂ h+

We no longer require that ∂ n = 0 on ∂  × (0, T ) and h(0) = 0 on , respectively, it follows that the boundary condition and initial condition hold, respectively.

y (x, t ) = ϕ (x, t ), y (x, 0) = y 0 ,

on ∂  × (0, T ),

(2.30)

x ∈ .

(2.31)

Thirdly, with respect to u, the known variational inequality follows. Evaluation of the variational inequality for D u L ( y , u , p ) and putting p = p 1 yields

T  D u L ( y , p , u )( w − u ) =

T 

α u ( w − u )dxdt + 0 

p ( w − u )dxdt ≥ 0,

∀ w ∈ U ad .

(2.32)

0 \

It follows that

(α u + p , w − u ) I ≥ 0, Since

∀ w ∈ U ad .

α > 0, then the optimal control has to obey the projection formula

1 u = P[ua ,ub ] − p , f or (x, t ) ∈  × (0, T ).

α

The assertion is proved.

2

(2.33)

(2.34)

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.6 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

6

Fig. 1. A typical interface element.

In order to derive the regularity of the optimal control problem (P) with a discontinuous coefficient β in the state equation, by the similar arguments in [13,2,4], we can prove the following lemma. Lemma 2.3. Assume that control u ∈ L 2 (0, T ; L 2 ()), f ∈ H 1 (0, T ; L 2 ()), g ∈ L 2 (0, T ; H 1/2 ()), ϕ ∈ L 2 (0, T ; H 1/2 (∂ )), y 0 ∈ ˜ 2 ()) ∩ H 1 (0, T ; H˜ 1 ()) ∩ H 1 () and the interface  ∈ C 2 . Then the problem (2.1)-(2.4) has a unique solution y ∈ L 2 (0, T ; H { y |∂  = ϕ } which satisfies for some constant C > 0

  y  H 2,1 (T ) ≤ C u L 2 (0,T ; L 2 ()) +  f L 2 (0,T ; L 2 ()) +  g L 2 (0,T ; H 1/2 ())  + ϕ L 2 (0,T ; H 1/2 (∂ )) + ∇ u 0 L 2 () .

(2.35)

˜ 2 ()) ∩ H 1 (0, T ; H˜ 1 ()) ∩ { y |∂  = 0}. Using the above lemma, we have the following regularity Let H 0 ( T ) = L 2 (0, T ; H result. 2, 1

2, 1

Theorem 2.4. The problem (P) has solution (u , y , p ) satisfying (u , y , p ) ∈ L 2 (0, T ; L 2 ()) × H 2,1 ( T ) × H 0 ( T ). Proof. It is obvious that u ∈ U ad ⊆ L 2 ( T ). Then by Lemma 2.3 and equation (2.9), we conclude that y ∈ H 2,1 ( T ). Hence, 2, 1 by using the fact that y − yd ∈ L 2 (0, T ; L 2 ()), we also have p ∈ H 0 ( T ). 2 3. The discretization and error estimates Since the analytic solution to problem (P) is rarely available, we seek numerical solutions. We use the variational discretization concept, see [11], where the state is replaced by a finite element approximation of the state equation (2.1)-(2.3). In order to get an accurate solution without using a body fitted mesh to resolve the interface, we use the immersed finite element method (IFEM) which is described in the next subsection. 3.1. The immersed finite element method

˜ 2 () because of the discontinuity in the coefficient, see The solution of interface problem (2.1)-(2.3) belongs to H Lemma 2.3. The standard finite element methods can not achieve the optimal convergence unless the mesh fits with the interface. The immersed finite element can achieve optimal convergence using a uniform mesh. Let Th be a uniform triangulation of  with mesh size h. If h is sufficiently small, then it is reasonable to assume the following: • The interface  will not intersect an edge of any element at more than two points unless the edge is part of ; • If  intersects the boundary of an element at two points, then these two points must be on different edges of this element. We call an element T an interface element if  passes through the interior of the element T ; otherwise we call T a non-interface element. Note that if  intersects ∂ T at two vertices of T , then the element is a non-interface element. The sets of all interface elements and non-interface elements are denoted by Thint and Thnon , respectively. On an element T ∈ Th , we consider the local finite element space ( T ,  T , T ), with

T = Span{1, x1 , x2 }, T = {ψT ( A i ) : i ∈ I , ∀ψT ∈ T },

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.7 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

7

where I = {1, 2, 3} and A i are vertices of T . It is well known that ( T ,  T , T ) has a set of shape functions ψi such that ψi ( A j ) = δi j , where δi j is the K ronecker delta function. For a non-interface element T ∈ Thnon , we use the standard piecewise linear polynomials as local basis functions and use S h ( T ) to denote the linear space spanned by the three nodal basis functions on T . For this space, we have the following well-known approximation property:

 v − I hT non v L 2 (T ) + h v − I hT non v  H 1 (T ) ≤ Ch2  v  H 2 (T ) , where I hT non : H 2 ( T ) −→ S h ( T ) is the interpolation operator. For interface element T ∈ Thint , we use piecewise linear polynomial separated by the interface  in the element. For a typical interface element  A 1 A 2 A 3 , its geometric configuration is given in Fig. 1. Assume that the interface  meets the triangle at points D and E, and the actual interface curve line  separates T into T + and T − . We let n( F ) be the normal of  at F and n D E be the unit normal vector on the line D E. Without loss of generality, we assume that β + ≥ β − and let

ρ=

β− β+

≤ 1. Furthermore, we partition I into two index sets: I + = {i : A i ∈ T + } and I − = {i : A i ∈ T − } according to the

locations of A i . The IFE functions introduced in [8] are constructed in the following piecewise polynomial form,



φ(x1 , x2 ) =

φ − (x1 , x2 ) ∈ T , if (x1 , x2 ) ∈ T − , φ + (x1 , x2 ) ∈ T , if (x1 , x2 ) ∈ T + ,

(3.1)

which satisfies

φ( A i ) = V i , i = 1, 2, 3, +

−

+

(3.2) −

φ ( D ) = φ ( D ), φ ( E ) = φ ( E ), ∂φ + ∂φ − β+ = β− , ∂ n( F ) ∂ n( F )

(3.3) (3.4)

where V i , i = 1, 2, 3 are given values, F is an arbitrary point on  ∩ T . We let  S h ( T ) denote the three dimensional linear space spanned by these piecewise linear functions. Without loss of generality, we assume that |I + | > |I − |. Then we can expand φ(x1 , x2 ) on T + . From equation (3.3), we have



φ( X ) =

+ if ( X ) ∈ T − , φ−( X ) = φ  ( X ) + c 0 L ( X ),  φ + ( X ) = i ∈I − c i ψi ( X ) + i ∈I + V i ψi ( X ), if ( X ) ∈ T + ,

where L ( X ) = n D E · ( X − D ) and X = (x1 , x2 ). By equation (3.4), we have

⎛

c0 = μ ⎝



c i ∇ψi · n( F ) +

i ∈I −

where

μ=



(3.5)

⎞

V i ∇ψi · n( F )⎠ ,

(3.6)

i ∈I +

1 − − ρ − 1. Using the above equation in (3.5) and setting φ ( A j ) = V j for j ∈ I leads to the following linear system

for c i , i ∈ I − :



   ψi ( A j ) + μ∇ψi · n( F ) L ( A j ) c i = V j − ψi ( A j ) + μ∇ψi · n( F ) L ( A j ) V i , j ∈ I − .

i ∈I −

i ∈I +

(3.7)

Due to ψi ( A j ) = δi j for i , j ∈ I − , the linear system (3.7) can be written in the following form,

(1 + μδ γ )c = b,

(3.8)



where c = c i , γ = ∇ψi · n( F ), δ = L ( A i ), and b = V i − μL( A i ) j ∈I + ∇ψ j · n( F ) V j , i ∈ I − . When h is sufficiently small, there exists one and only one IFE function φ(x1 , x2 ) in the form of (3.1) satisfying (3.2)-(3.4). Applying the Sherman-Morrison formula, the solution to (3.8) can be expressed explicitly as

c=b−μ Then, we have

S h (T ) =

γ bδ 1 + μγ δ



(3.9)

.

S h ( T ), if T ∈ Thnon ,  S h ( T ), if T ∈ Thint . J

Definition 3.1. (IFE space) The IFE space V h () is defined as the set of all piecewise linear functions that satisfy

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.8 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

8

• φ| T ∈ S h ( T ), • φ is continuous at all nodal points, • φ(xb ) = 0 if xb is a nodal point on ∂ . Following the idea from [9,8], for an interface element T ∈ Thint , we similarly obtain the following well-known approximation property:

 v − I hT int v L 2 (T ) + h v − I hT int v  H 1 (T ) ≤ Ch2  v  H 2 (T ) , J

where I h int : H 2 ( T ) −→ S h ( T ) is the interpolation operator. It has been proved in [8] that the IFE space V h () has an T

J

optimal approximation capability. We note that a function in the IFE space V h () is likely discontinuous across the common J

edges of two adjacent interface elements, see [8] for the detailed description of the phenomenon. In other words, V h () ∈ / H 01 (). Thus we define a larger space J

H h () := H 01 () + V h (), which is endowed with the broken H 1 norms as | v |21,h :=



T ∈Th

| v |2H 1 (T ) and  v 21,h :=



T ∈Th

 v 2H 1 (T ) . Due to the non-

conforming property of the IFE space, the following finite element approximation for (2.1)-(2.4) is called nonconforming immersed finite element method. J

Problem 3.2. (IFE approximation) For any u ∈ L 2 ( T ), find y h (t ) ∈ H 1 (0, T ; V h ()) such that J

(∂t yh , v h ) + a( yh , v h ) = (u , v h ) −  g , v h  ∀ v h ∈ V h () where

a(u h , v h ) :=



with yh (0) = y 0 ,

(3.10)

β∇ uh · ∇ v h dx ∀uh , v h ∈ H h ().

T ∈Th T

We list a useful lemma results in the following: Lemma 3.3. [1] (Discrete Poincaré inequality) There exists a constant C independent of h and the interface  such that

φL 2 () ≤ Cah (φ, φ),

J

∀φ ∈ V h ().

(3.11)

From [1], we similarly obtain the following theorem result, J

Theorem 3.4. (Error estimates) Let y ∈ H 2,1 ( T ) and y h ∈ H 1 (0, T ; V h ()) for all t ∈ [0, T ] be the solutions of (2.1)-(2.4) and (3.10) respectively. Then there exists a constant C > 0 such that

   y − yh (t )1,β ≤ Ch  y  H˜ 2 () +  yt  H˜ 2 () ,    y − yh (t ) L 2 () ≤ Ch2  y  H˜ 2 () +  yt  H˜ 2 () ,

(3.12) (3.13)

where  y 21,β = β + ∇ y 2L 2 (+ ) + β − ∇ y 2L 2 (− ) . 3.2. Semi-discretization of the optimal control problem and error estimates We are now ready to define the semi-discrete problem (P˜ h ) using the variational discretization concept introduced in [11]. Problem 3.5. (P˜ h ) Consider the problem of minimizing

J h ( yh , u ) =

1

T 2

 yh − yd  dt +

2

α

0

over all ( y h (t ), u ) ∈ H

1

T u 2 dt ,

2

(3.14)

0 J (0, T ; V h ()) ×

L (0, T ; L 2 ()) subject to 2

(∂t yh , v h ) + ah ( yh , v h ) = ( f + u , v h ) −  g , v h  ,

J

∀ v h ∈ V h (),

with yh (0) = y 0 ,

(3.15)

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.9 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

9

and the control constraints

ua ≤ u ≤ ub .

(3.16)

Similar to Theorem 2.2, we have the following theorem for (P˜ h ). Theorem 3.6. The problem (P˜ h ) has a unique solution u h ∈ L 2 ( T ) with associated state y h ∈ H 1 (0, T ; V h ()) and adjoint state J

J

p h ∈ H 1 (0, T ; V h ()) that satisfy the state equation J

(∂t yh , v h ) + ah ( yh , v h ) = ( f + uh , v h ) −  g , v h  ,

∀ v h ∈ V h (),

with y h (0) = y 0 ,

(3.17)

the adjoint equation J

−(∂t ph , v h ) + ah ( v h , ph ) = ( v h , yh − yd ), ∀ v h ∈ V h (),

with p h ( T ) = 0,

(3.18)

and the variational inequality

T (α uh + ph , w − uh )dt ≥ 0, ∀ w ∈ U ad .

(3.19)

0

Moreover, the variational inequality is equivalent to the projection equation



u h = P[ u a , u b ] −

1

α



ph .

(3.20)

From the theorem, we can see that the control is implicitly discretized by projecting a discrete adjoint state onto U ad . J To get error estimates between (P) and (P˜ h ), we introduce auxiliary functions y˜ ∈ H 1 (0, T ; V h ()) and p˜ ∈ J

H 1 (0, T ; V h ()) which are solutions of the following problems J

(∂t y˜ , v h ) + ah ( y˜ , v h ) = ( f + u , v h ) −  g , v h  , ∀ v h ∈ V h (), − (∂t p˜ , v h ) + ah ( v h , p˜ ) = ( v h , y − yd ), ∀ v h ∈

J V h ()

with y˜ (0) = y 0 ,

(3.21)

with p˜ ( T ) = 0.

(3.22)

Using Theorem 3.4, we have the following lemma 2, 1

Lemma 3.7. Let y ∈ H 2,1 ( T ) and p ∈ H 0 ( T ) be the solution of (2.6) and (2.10), for all t ∈ [0, T ], respectively. y˜ ∈ J (0, T ; V h ())

J (0, T ; V h ())

H and p˜ ∈ H be the solution of (3.21) and (3.22), for all t ∈ [0, T ], respectively. Then there exists a constant C > 0, independent of h, such that 1

1

 y (t ) − y˜ (t ) + h y (t ) − y˜ (t )1,β ≤ Ch2 ,

∀t ∈ [0, T ],

(3.23)

 p (t ) − p˜ (t ) + h p (t ) − p˜ (t )1,β ≤ Ch ,

∀t ∈ [0, T ].

(3.24)

2

Theorem 3.8. Let (u , y , p ) and (u h (t ), y h (t ), p h (t )) be the solutions of the problems (P) and (P˜ h ) respectively. Then there exists a constant C > 0, independent of h, such that

T

T 2

2

u − uh  dt +

α 0

 y − yh  dt ≤ 0

1

T

T  p − p˜  dt +

α 0

⎛

max  y (t ) − yh (t )2 ≤ C ⎝ max  y (t ) − y˜ (t )2 +

0≤t ≤ T

0

T

0≤t ≤ T

⎛

max  p (t ) − p h (t )2 ≤ C ⎝ max  p (t ) − p˜ (t )2 +

0≤t ≤ T

⎞

0

T

0≤t ≤ T

⎛

0

T

0≤t ≤ T

0

(3.25)

u − uh 2 dt ⎠ ,

max  y (t ) − yh (t )21,β ≤ C ⎝ max  y (t ) − y˜ (t )21,β +

0≤t ≤ T

 y − y˜ 2 dt ,

2

(3.26)

⎞ u − uh 2 dt +  y˜ (0) − yh (0)21,β ⎠ ,

(3.27)

⎞

 y − yh 2 dt ⎠ ,

(3.28)

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.10 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

10

⎛ max  p (t ) −

0≤t ≤ T

p h (t )21,β

≤ C ⎝ max  p (t ) − 0≤t ≤ T

⎞

T p˜ (t )21,β

 y − yh  dt +  p˜ ( T ) − 2

+

p h ( T )21,β ⎠ .

(3.29)

0

Proof. The error equations can be obtained from (3.17), (3.21) and (3.18), (3.22), respectively, J (∂t ( y˜ − yh ), v h ) + ah ( y˜ − yh , v h ) = (u − uh , v h ), ∀ v h ∈ V h (), with y˜ (0) − yh (0) = 0,

−(∂t ( p˜ − ph ), v h ) + ah ( v h , p˜ − ph ) = ( v h , y − yh ), ∀ v h ∈

J V h (),

with p˜ ( T ) − p h ( T ) = 0.

(3.30) (3.31)

We follow the idea of the proof in [12] to prove the first inequality. Taking w = u h in (2.11) and w h = u in (3.19) and summing up these two inequalities, we obtain

(α (u − uh ) + ( p − ph ), uh − u ) I ≥ 0.

(3.32)

It follows that

T

T 2

u − uh  dt ≤

α 0

(u − uh , ph − p )dt 0

T

T (u − uh , ph − p˜ )dt +

= 0

(u − uh , p˜ − p )dt

(3.33)

0

T (u − uh , ph − p˜ )dt +

≤

α

T 2

u − uh  dt +

2

0

T

1

 p − p˜ 2 dt .

2α

0

0

From error equation (3.30) with v h = p h − p˜ and error equation (3.31) with v h = y˜ − y h , we have

T

T (u − uh , ph − p˜ )dt =

0

T (∂t ( y˜ − yh ), ph − p˜ )dt +

0

ah ( y˜ − yh , p h − p˜ )dt 0

T

T ( y˜ − yh , ∂t ( ph − p˜ ))dt +

=− 0

ah ( y˜ − yh , p h − p˜ )dt 0

(3.34)

T ( y − yh , y˜ − yh )dt

=− 0

≤−

1

T  y − yh 2 dt +

2 0

1

T  y − y˜ 2 dt .

2 0

The estimate (3.25) then follows from (3.33) and (3.34). To prove the second inequality (3.26), we first use the triangle inequality

max  y (t ) − yh (t ) ≤ max  y (t ) − y˜ (t ) + max  y˜ − yh (t ),

0≤t ≤ T

0≤t ≤ T

(3.35)

0≤t ≤ T

for any t ∈ (0, T ], integral from 0 to t on the both sides of error equation (3.30) with v h = y˜ − y h , we have

1 2

1

t

 y˜ (t ) − yh (t )2 −  y˜ (0) − yh (0)2 +

t  y˜ − yh 21,β dτ =

2

0

(u − uh , y˜ − yh )dτ 0

≤

1

T u − uh 2 dτ +

2 0

1

(3.36)

T  y˜ − yh 2 dτ .

2 0

Noting that y˜ (0) − y (0) = 0, using Gronwall inequality, and combining (3.35), we can obtain the error estimate (3.26).

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.11 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

11

In order to get (3.27), for any t ∈ (0, T ], integral from 0 to t on the both sides of error equation (3.30) with v h =

∂t ( y˜ − yh ), we have

t

t

1

∂t ( y˜ − yh ) dτ +  y˜ (t ) − 2

2

yh (t )21,β

0

1

(u − uh , ∂t ( y˜ − yh ))dτ +  y˜ (0) − yh (0)21,β

=

2

0

≤

1

t

ε

u − uh 2 dτ +

2ε

1

∂t ( y˜ − yh )2 dτ +  y˜ (0) − yh (0)21,β .

2

0

(3.37)

t

2

0

Taking appropriate ε , and using the triangle inequality, error estimate (3.27) holds. To get (3.28), for any t ∈ [0, T ), integral from t to T on the both sides of error equation (3.31) with v h = p˜ − p h , we have

1 2

T

1

 p˜ (t ) − ph (t )2 −  p˜ ( T ) − ph ( T )2 +

T  p˜ − ph 21,β dτ =

2

t

( y − yh , p˜ − ph )dτ t

≤

1

T  y − yh 2 dτ +

2

1

 p˜ − ph 2 dτ .

2

0

(3.38)

T 0

Noting that p˜ ( T ) − p ( T ) = 0, using Gronwall inequality, and combining the triangle inequality, we can obtain the error estimate (3.28). Next, we will give the proof of (3.29). For any t ∈ [0, T ), integral from T to t on the both sides of error equation (3.31) with v h = ∂t ( p˜ − p h ), we have

T

T

1

∂t ( p˜ − ph ) dτ +  p˜ (t ) − 2

2

p h (t )21,β

t

1

( y − yh , ∂t ( p˜ − ph ))dτ +  p˜ ( T ) − ph ( T )21,β

=

2

t

≤

1

T  y − yh 2 dτ +

2ε

ε

1

∂t ( p˜ − ph )2 dτ +  p˜ ( T ) − ph ( T )21,β ,

2

t

(3.39)

T

2

t

taking appropriate ε , and using the triangle inequality, we obtain error estimate (3.29).

2

As a consequence of Theorem 3.8 and Lemma 3.7, we have the following results. Theorem 3.9. Let (u , y , p ) and (u h , y h , p h ) be the solutions of the optimal control problem (P) and (P˜ h ) respectively. Then there exists a generic constant C , independent of the mesh size h, such that

T u − uh 2 dt ≤ Ch4 ,

max  y (t ) − yh (t ) ≤ Ch2 ,

0≤t ≤ T

max  y (t ) − yh (t )1,β ≤ Ch,

0≤t ≤ T

(3.40)

0

max  p (t ) − p h (t ) ≤ Ch2 ,

0≤t ≤ T

max  p (t ) − p h (t )1,β ≤ Ch.

0≤t ≤ T

(3.41)

Remark 3.10. In the case of U ad = L 2 () (unconstrained problem), the projection equations (2.12) and (3.20) become u ∗ = − α1 p ∗ and uh∗ = − α1 ph∗ , respectively. Using the properties of p ∗ and ph∗ , we then have the regularity of the control u ∗ ∈

˜ 2 () and the following H 1 error estimate H

u (t ) − uh (t )1,β ≤ Ch. 3.3. Full-discretization of the optimal control problem and error estimates We are now ready to define the full discrete problem (Ph ) using the variational discretization concept introduced in [11]. Problem 3.11. (Ph ) Consider the problem of minimizing

1 N

J h ( ynh , u ) =

2



n =1 I

n

α N

 ynh − yd 2 dt +

2



n =1 I

u 2 dt , n

(3.42)

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.12 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

12 J

over all ( ynh , u ) ∈ L 2 (0, T ; V h ()) × L 2 (0, T ; L 2 ()) subject to

(∂ t ynh , v h ) + ah ( ynh , v h )

tn

1

=

tn

( f + u , v h )dt − t n −1

tn

1

tn

J

ynh − ynh−1 tn

∀ v h ∈ V h (),

(3.43)

t n −1

n = 1, 2, · · · , N , where on each time interval I n , ynh ∈ V h (), ∂ t ynh =

J

 g , v h  dt ,

with yh0 = h y 0 ,

and the control constraints

ua ≤ u ≤ ub .

(3.44)

Similar to Theorem 2.2, we have the following theorem for (Ph ). J

Theorem 3.12. The problem (Ph ) has a unique solution u h ∈ L 2 ( T ) with associated state ynh ∈ L 2 (0, T ; V h ()) and adjoint state ph ∈

J L 2 (0, T ; V h ())

that satisfy the state equation

tn

1

(∂ t ynh , v h ) + ah ( ynh , v h ) =

tn

( f + uh (t ), v h )dt − t n −1

tn

1

tn

J

 g , v h  dt ,

∀ v h ∈ V h (),

t n −1

n = 1, 2, · · · , N ,

(3.45)

with y h0 = h y 0 ,

the adjoint equation

−(∂ t pnh , v h ) + ah ( v h , pnh−1 ) = ( v h , ynh − yd ), ∀ v h ∈ V h (), J

n = 1, 2, · · · , N ,

with p hN = 0,

(3.46)

and the variational inequality N   n =1 I

(α uh (t ) + pnh−1 , w h − uh (t ))dt ≥ 0, ∀ w h ∈ U ad .

(3.47)

n

Moreover, the variational inequality is equivalent to the projection equation



unh = P[ua ,ub ] −

J

1

α



pnh−1 .

(3.48)

From the theorem, we can see that the control is implicitly discretized by projecting a discrete adjoint state onto U ad . J To get error estimates between (P) and (Ph ), we introduce the following auxiliary functions y˜ ∈ V h () and p˜ ∈

V h (),

n = 1, 2, · · · , N which are the full discrete solutions of the problems (3.21)

(∂ t y˜ n , v h ) + ah ( y˜ n , v h ) =

1

tn

tn ( f + u (t ), v h )dt − t n −1

1

tn

tn

J

 g , v h  dt , ∀ v h ∈ V h (), t n −1

n = 1, 2, · · · , N ,

−(∂ t p˜ , v h ) + ah ( v h , p˜ n

n −1

with y˜ 0 = h y 0 ,

n

) = ( v h , y (t ) − yd ), ∀ v h ∈ n = 1, 2, · · · , N ,

J V h (), N

with p˜ = 0.

(3.49)

(3.50)

Using Theorem 3.4, we have the following lemma J

2, 1 Lemma 3.13. Let y ∈ H 2,1 ( T ) and p ∈ H 0 ( T ) be the solution of (2.6) and (2.10), for all t ∈ [0, T ], respectively. y˜ n ∈ V h () J

and p˜ n ∈ V h () be the solution of (3.49) and (3.50), n = 1, 2, · · · , N, respectively. Then there exists a constant C > 0, independent of h and tn , such that

 y (t n ) − y˜ n  + h y (t n ) − y˜ n 1,β ≤ C (h2 + t ),

n = 1, 2, · · · , N ,

(3.51)

 p (t ) − p˜  + h p (t ) − p˜ 1,β ≤ C (h + t ),

n = 1, 2, · · · , N .

(3.52)

n

n

where t = max tn . 1≤n≤ N

n

n

2

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.13 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

13

Theorem 3.14. Let (u , y , p ) and (u h (t n ), ynh , pnh ) be the solutions of the problems (P) and (Ph ) respectively. Then there exists a constant C > 0, independent of h and tn , such that

α

N   n =1 I

2

u (t ) − uh (t ) dt +

N 

n

 y (t ) −

ynh 2 tn

≤

n =1

n

⎛

α n =1

N  

⎜

max  y (t n ) − ynh 2 ≤ C ⎝ max  y (t n ) − y˜ n 2 +

1≤n≤ N

N  1

1≤n≤ N

n =1 I

⎛ ⎜

max  y (t n ) − ynh 21,β ≤ C ⎝ max  y (t n ) − y˜ n 21,β +

1≤n≤ N

1≤n≤ N

n

max  p (t ) −

1≤t ≤ N

pnh 2

≤C

max  p (t ) − p˜  + n

n 2

1≤t ≤ N

max  p (t ) −

1≤t ≤ N

pnh 21,β

≤C

 y (t n ) − y˜ n 2 tn ,

In

⎞ (3.54)

n

N  

n

max  p (t ) −

1≤t ≤ N

p˜ n 21,β

+

(3.53)

n =1

⎞ ⎟ u (t ) − uh (t )2 dt +  y˜ 0 − yh0 21,β ⎠ ,

(3.55)

n

 n

 y (t ) −

ynh 2 tn

(3.56)

,

n =1

 n

N 

N 

⎟ u (t ) − uh (t )2 dt ⎠ ,

n =1 I



 p (t ) − p˜ n−1 2 dt +

N 

 n

 y (t ) −

+  p˜ −

ynh 2 tn

N

p hN 21,β

.

(3.57)

n =1

Proof. The fully discrete error equations can be obtained from (3.45), (3.49) and (3.46), (3.50), respectively,

(∂ t ( y˜ n − ynh ), v h ) + ah ( y˜ n − ynh , v h ) =

tn

1

tn

J

(u (t ) − uh (t ), v h ), ∀ v h ∈ V h (), t n −1

n = 1, 2, · · · , N , with y˜ 0 − yh0 = 0.

−(∂ t ( p˜ n − pnh ), v h ) + ah ( v h , p˜ n−1 − pnh−1 ) = ( v h , y (t n ) − ynh ), ∀ v h ∈

(3.58)

J V h ()

n = 1, 2, · · · , N , with p˜ N − p hN = 0.

(3.59)

Similarly with the proof of semi-discrete case, taking w = u h in (3.47) and w h = u in (3.19) and summing up these two inequalities, we obtain N    n =1 I



α (u (t ) − uh (t )) + ( p (t ) − pnh−1 ), uh − u dt ≥ 0.

(3.60)

n

It follows that

α

N   n =1 I

=

u (t ) − uh (t )2 dt ≤

≤

n =1 I

n

N   n =1 I

( pnh−1 − p (t ), u (t ) − uh (t ))dt

n

( pnh−1 − p˜ n−1 , u (t ) − uh (t ))dt +

n

N   n =1 I

n

N   n =1 I

N  

( pnh−1 − p˜ n−1 , u (t ) − uh (t ))dt +

( p˜ n−1 − p (t ), u (t ) − uh (t ))dt

n

N  1 

2α

n =1 I

n

 p˜ n−1 − p (t )2 dt +

α N

2

(3.61)



n =1 I

u (t ) − uh (t )2 dt . n

In order to estimate the first term in (3.61), from error equation (3.58) with v h = pnh−1 − p˜ n−1 and error equation (3.59) with v h = y˜ n − ynh , using discrete partial sum formula, we have

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.14 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

14

N   n =1 I

( pnh−1 − p˜ n−1 , u (t ) − uh (t ))dt =

N N   (∂ t ( y˜ n − ynh ), pnh−1 − p˜ n−1 )tn + ah ( y˜ n − ynh , pnh−1 − p˜ n−1 )tn n =1

n

n =1

N N   =− ( y˜ n − ynh , ∂ t ( pnh − p˜ n ))tn + ah ( y˜ n − ynh , pnh−1 − p˜ n−1 )tn n =1

n =1

(3.62)

N  =− ( y (t n ) − ynh , y˜ n − ynh )dt n =1 N

≤−

1 2

1 N

 y (t n ) − ynh 2 tn +

n =1

2

 y (t n ) − y˜ n 2 tn ,

n =1

then it follows from (3.61) and (3.63) that

α N

2



n =1 I

≤

n

1 2α

1 N

u (t ) − uh (t )2 dt + N   n =1 I

2

 y (t n ) − ynh 2 tn

n =1

 p˜ n−1 − p (t )2 dt +

n

1

(3.63)

N

2

 y (t n ) − y˜ n 2 tn .

n =1

The estimate (3.53) is obtained. To prove the second inequality (3.54), we first use the triangle inequality

max  y (t n ) − ynh  ≤ max  y (t n ) − y˜ n  + max  y˜ n − ynh .

1≤n≤ N

1≤n≤ N

(3.64)

1≤n≤ N

For any 1 ≤ n ≤ N, multiplying tn and summing over n on the both sides of error equation (3.58) with v h = y˜ n − ynh , we have

1 2

1

 y˜ n − ynh 2 −  y˜ 0 − yh0 2 + 2

n 

 y˜ n − ynh 21,β tn =

n =1

N   n =1 I

≤

(u (t ) − uh (t ), y˜ n − ynh )dτ

n

N  1

2

n =1 I

1 N

u (t ) − uh (t )2 dτ +

2

n



n =1 I

 y˜ n − ynh 2 tn . n

(3.65) Noting that y˜ 0 − = 0, using Gronwall inequality, and combining (3.64), we can obtain the error estimate (3.54). In order to get (3.55), for any 1 ≤ n ≤ N, multiplying tn and summing over n on the both sides of error equation (3.58) with v h = ∂ t ( y˜ n − ynh ), we have y h0

n 

1

∂ t ( y˜ n − ynh )2 tn +  y˜ n − ynh 21,β ≤ 2

n =1

≤

n  1 

2ε

n =1 I

n   n =1 I

1

(u (t ) − uh (t ), ∂t ( y˜ n − ynh ))dτ +  y˜ 0 − yh0 21,β 2

n

ε

(3.66)

n

2

u (t ) − uh (t ) dτ +

n

2

∂t ( y˜ − n

n =1

ynh )2 tn

1

+  y˜ − 0

2

y h0 21,β .

Taking appropriate ε , and using the triangle inequality, error estimate (3.55) holds. To get (3.56), for any 1 ≤ n ≤ N, multiplying tn and summing over n from n to N on the both sides of error equation (3.59) with v h = p˜ n−1 − pnh−1 , we have

1 2

1

 p˜ n − pnh (t )2 −  p˜ N − phN 2 + 2

N  n=n

 p˜ n − pnh 21,β tn =

N  ( y (t n ) − ynh , p˜ n−1 − pnh−1 )tn n=n

1 N

≤

2

n =1

1 N

 y (t n ) − ynh 2 tn +

2

 p˜ n−1 − pnh−1 2 tn .

n =1

(3.67)

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.15 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

15

Noting that p˜ N − p hN = 0, using Gronwall inequality, and combining the triangle inequality, we can obtain the error estimate (3.56). Next, we will give the proof of (3.57). for any 1 ≤ n ≤ N, multiplying tn and summing over n from n to N on the both sides of error equation (3.59) with v h = ∂ t ( p˜ n − pnh ), we have N  n=n

1

∂ t ( p˜ n − pnh )2 tn +  p˜ n − pnh 21,β ≤ 2

≤

N 1 

2ε

N  1 ( y (t n ) − ynh , ∂ t ( p˜ n − pnh ))tn +  p˜ N − phN 21,β

2

n=n

ε

(3.68)

N

n

 y (t ) −

ynh 2 tn

n =1

+

2 n=n

∂ t ( p˜ − n

pnh )2 tn

1

+  p˜ − N

2

Taking appropriate ε , and using the triangle inequality, we obtain error estimate (3.57).

p hN 21,β .

2

As a consequence of Theorem 3.14 and Lemma 3.13, we have the following results. Theorem 3.15. Let (u , y , p ) and (unh , ynh , pnh ) be the solutions of the optimal control problem (P) and (Ph ) respectively. Then there exists a generic constant C , independent of the mesh size h, such that

⎛

⎞1/2 N   ⎜ ⎟ u (t ) − uh (t )2 dt ⎠ ≤ C (h2 + t ), ⎝ n =1 I

(3.69)

n

max  y (t n ) − ynh  ≤ C (h2 + t ),

max  y (t n ) − ynh 1,β ≤ C (h + t ),

1≤n≤ N

max  p (t n ) − pnh  ≤ C (h2 + t ),

1≤n≤ N

(3.70)

1≤n≤ N

max  p (t n ) − pnh 1,β ≤ C (h + t ).

(3.71)

1≤n≤ N

Remark 3.16. In the case of U ad = L 2 () (unconstrained problem), the projection equations (2.12) and (3.48) become u ∗ = − α1 p ∗ and uh∗ = − α1 ph∗ , respectively. Using the properties of p ∗ and ph∗ , we then have the regularity of the control u ∗ ∈

˜ 2 () and the following H 1 error estimate H

⎛ ⎞1/2 N   ⎜ ⎟ u (t ) − uh (t )21,β dt ⎠ ≤ C (h + t ). ⎝ n =1 I

n

Remark 3.17. If we adopt the Crank-Nicolson time stepping method for solving this problem, following the idea of above numerical analysis, we can obtain the convergence rate O (h2 + t 2 ) with the norm  ·  and O (h + t 2 ) with the norm  · 1,β for the state, adjoint state and control, respectively. 4. Implementation details In this section, we give details about how to find the solution of the finite dimensional optimal control problem (Ph ). We consider in the sequel a general problem where the state equation of the problem (P) is replaced by

∂t y − ∇ · (β∇ y ) = f + u

in (\) × (0, T ),

y=ϕ

on ∂  × (0, T ),

[ y]×(0, T ) = 0, [β∂n y]×(0, T ) = g

(4.1)

y (0) = y 0 .

(4.2)

Note that the adjoint equation and the optimality system are the same. The theoretical results can be generalized to the 3

general one if f ∈ L 2 () and g ∈ H 2 (∂ ). In order to find the solution of the problem (Ph ), we then have to solve the coupled problem

(∂ t ynh , v h ) + ah ( ynh , v h )

=

1

tn

tn ( f + uh (t ), v h )dt − t n −1

1

tn

tn

pnh−1 )

=

( v h , ynh

− yd ), ∀ v h ∈

J V h (),

∀ v h ∈ V h (),

t n −1

n = 1, 2, · · · , N ,

−(∂ t pnh , v h ) + ah ( v h ,

J

 g , v h  dt ,

with yh0 = h y 0 ,

n = 1, 2, · · · , N ,



unh = P[ua ,ub ] −

1

α



pnh−1 ,

p hN

(4.3)

= 0,

(4.4)

n = 1, 2, · · · , N .

(4.5)

with

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.16 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

16 J

J

Note that the function ynh ∈ V h,ϕ () satisfies the discrete boundary condition. Let {φ1 , · · · , φm } be basis functions of V h (). For the boundary condition to be satisfied, we also define functions φm+1 , · · · , φm+l so that the boundary data. The function ynh then can be written as

ynh =

m 

m+l 

Y n ( j )φ j +

j =1

m+l

j =m+1 r (x j )φ j

r n (x j )φ j .

interpolates

(4.6)

j =m+1

Substituting (4.6) into (4.3) and choosing v h = φi , i = 1, · · · , m, we get the matrix-vector form

( M + tn A )Y n = F 1n−1 + MY n−1 + M tn U n ,

(4.7)

where

A (i , j ) = ah (φ j , φi ), Y n ( j ) = ynh (x j ), M (i , j ) = (φ j , φi ), U n ( j ) =



F 1n−1 (i ) =

 ( f , φi )dt −

In

 g , φi  dt −

j = 1, · · · , m ,

tn

r n (x j )ah (φ j , φi )tn +

j =m+1

In

i = 1, · · · , m ,

m+l 



1

u h (x j , t )dt , In m+l 

(rn−1 (x j ) − rn (x j ))(φ j , φi ),

(4.8)

j =m+1

n = 1, · · · , N .

Similarly, we obtain the matrix-vector form of (4.4),

( M + tn A ) P n−1 = M tn Y n + M P n + F 2n ,

(4.9)

where

P n−1 ( j ) = pnh−1 (x j ),

m+l 

F 2n (i ) =

r n (x j )(φ j , φi )tn − ( yd , φi )tn ,

j =m+1

i = 1, · · · , m ,

j = 1, · · · , m ,

(4.10)

n = 1, · · · , N .

The numerical integrations needed in the computation are approximated using the three-point Gaussian quadrature formula. For an interface triangle T ∈ Thint , it is divided into one sub-triangle and two sub-triangle with curved sides according to the approximated interface for the numerical integration. The parabolic interface optimal control problem has essential difference with the elliptic interface optimal control problem [25]. The main challenge for solving the parabolic interface optimal control problem results in the fact that the state y and the adjoint state p are marching in opposite directions and the interface problem. The numerical discretizations of our approach will describe the details in the following for the cases with non-constraints and constraints, respectively. In the case of U ad = L 2 () (unconstrained problem), the projection equation (4.5) becomes

unh = −

1

α

pnh−1 ,

(4.11)

The vector form is U n = − α1 P n−1 . Therefore we have the following large scale symmetric indefinite linear system of equations

⎡

⎤

Un n⎣ Y n ⎦ = F n, Sh P n −1

(4.12)

where the matrix S nh , (U n , Y n , P n−1 ) T and F n have the following block structured matrix, respectively.

⎛

S nh

A1 ⎜ E ⎜

⎜ =⎜ ⎜ ⎝ n

n

(U , Y , P

⎞

E A2

E

..

..

.

..

.

E n −1 T

⎟ ⎟ ⎟ ⎟, ⎟ E ⎠

.

A n −1 E

An

) = ( U , Y , P , U 2 , Y 2 , P 1 , · · · , U N , Y N , P N −1 ) T , 1

1

0

F n = (0, F 21 , F 10 + MY 0 , 0, F 22 , F 11 , · · · , 0, F 2N + M P N , F 1N −1 ) T , where

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.17 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

⎡ An = ⎣

−α M tn

⎡

0 E =⎣0 0

0 − M tn

17

⎤

0

− M tn M + A tn ⎤

− M tn M + A tn ⎦ 0

n = 1, 2, · · · , N ,

0 0 0 −M ⎦ . 0 0

In the case of U ad = {u ∈ L 2 () : ua ≤ u ≤ u b } (control-constrained problem), the problem (Ph ) leads to a single equation for the optimal control unh , that is,



G (u ) := unh − P[ua ,ub ] −

1

α



pnh−1 (u ) = 0,

n = 1, 2, · · · , N ,

(4.13)

where pnh−1 (u ) is the solution of (4.4) and (4.5) with unh = u. An iteration algorithm can be used to solve the nonlinear and non-smooth equation. Algorithm 1. Give an initial function u 0 ∈ L 2 ( T ). J 2. Get ynh ∈ V h () by solving (4.3). J

3. Get pnh ∈ V h () by solving (4.4). 4. Get un1 by the equation (4.5), where P[a,b] ( v (x)) := max{a, min{ v (x), b}}. 5. If |u 0 − u 1 | ≤ 1.0 × 10−6 , then we set unh = u 1 , else we set u 0 = u 1 and goto step 2. It has been shown in [12] that this algorithm is convergent if the parameter α is large enough. 5. Numerical experiments In this section we first present two test examples for which the exact solutions are known explicitly to verify our theory. In all examples, the computation domain  is a square [−1, 1] × [−1, 1], so that the uniform triangulation with 2N elements in each direction can be used. All the experiments are computed with double precision. The time step is chosen as t = O (h2 ). In all tables listed in this section, “Order” represents the convergence order which is evaluated by

Order =

1 log 2

log

u − u N  , u − u 2N 

for a specific norm  · . 5.1. An example without control constraints Example 1. In this example we set U ad = L 2 (). The interface is a circle centered at the origin with radius r0 = 0.5. We use the method in [23] to construct the analytical state, adjoint state and control. Let

⎧ sin(t − T )(x21 + x22 − r02 )(x1 − 1)(x1 + 1)(x2 − 1)(x2 + 1) ⎪ ⎪ ⎪ in − , ⎨ β− ϕ (x) = (5.1) ⎪ sin(t − T )(x21 + x22 − r02 )(x1 − 1)(x1 + 1)(x2 − 1)(x2 + 1) ⎪ + ⎪ ⎩ in  . β+ Then we have u ∗ = ϕ and p ∗ = −α u ∗ . It is easy to verify that p ∗ satisfies the interface jump conditions and the homoge-

neous Dirichlet boundary condition. We choose

⎧ cos(t − T )(x21 + x22 − r02 )(x1 − 1)(x1 + 1)(x2 − 1)(x2 + 1) ⎪ ⎪ ⎪ in − , − ⎨ β− ∗ y (x) = ⎪ ⎪ cos(t − T )(x21 + x22 − r02 )(x1 − 1)(x1 + 1)(x2 − 1)(x2 + 1) ⎪ ⎩− in + . β+

(5.2)

Hence, the functions f and g can be determined by (4.1), y ∗ and u ∗ . The desired state yd can be determined by (2.13), α = 0.1 in this example. An illustration of the geometry of the interface problem and a triangulation is depicted in Fig. 2. We list some numerical results obtained by the IFEM for three cases with different coefficient jump ratios: β − /β + = 1/1000 (large jump), β − /β + = 1000 (large jump), and β − /β + = 1/10 (moderate jump) in Table 1, Table 2, and Table 3, respectively. From these tables, one can find that the numerical results are in accordance with the theoretical findings. Fig. 3 depicts the development of the L ∞ (0, T ; L 2 ()) and L ∞ (0, T ; H 1 ) error for the control and state under uniform refinement of the mesh with β − /β + = 1/1000. From the figure, the expected order O (h2 ) in L ∞ (0, T ; L 2 ) norm for the control and state is observed, and the order O (h) in L ∞ (0, T ; H 1 ) norm is shown. y ∗ and p ∗ . We set

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.18 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

18

Fig. 2. The geometry of the domain and the interface. A triangulation with N = 5 is plotted.

Table 1 Grid refinement analysis for Example 1 with β − = 1, β + = 1000. (a) Errors in the L ∞ ( L 2 )-norm N

u ∗ − uh∗  L ∞ ( L 2 () )

16 32 64 128

3.3383E-03 7.2702E-04 1.8843E-04 4.0460E-05

Order

 y ∗ − yh∗  L ∞ ( L 2 )

2.19 1.95 2.21

5.0239E-03 1.1440E-03 2.8516E-04 6.4216E-05

Order

 p ∗ − ph∗  L ∞ ( L 2 )

Order

2.13 2.00 2.15

3.3383E-04 7.2702E-05 1.8843E-05 4.0460E-06

2.19 1.95 2.21

(b) Errors in the L ∞ ( H 1 )-norm N

|u ∗ − uh∗ | L ∞ ( H 1 )

16 32 64 128

8.6818E-02 3.9402E-02 1.8508E-02 9.0425E-03

Order

| y ∗ − yh∗ | L ∞ ( H 1 )

1.14 1.09 1.03

1.0310E-01 4.6824E-02 2.1996E-02 1.0746E-02

Order

 p ∗ − ph∗  L ∞ ( H 1 )

Order

1.14 1.09 1.03

8.6818E-03 3.9402E-03 1.8508E-03 9.0425E-04

1.14 1.09 1.03

Order

 p ∗ − ph∗  L ∞ ( L 2 )

Order

1.97 1.88 1.92

9.3133E-04 2.3502E-04 6.1025E-05 1.5247E-05

1.99 1.95 2.00

Order

 p ∗ − ph∗  L ∞ ( H 1 )

Order

0.99 1.00 1.00

2.7060E-02 1.3631E-02 6.8122E-03 3.4067E-03

0.99 1.00 1.00

Table 2 Grid refinement analysis for Example 1 with β − = 1000, β + = 1. (a) Errors in the L ∞ ( L 2 )-norm N

u ∗ − uh∗  L ∞ ( L 2 )

16 32 64 128

9.3133E-03 2.3502E-03 6.1025E-04 1.5247E-04

Order

 y ∗ − yh∗  L ∞ ( L 2 )

1.99 1.95 2.00

1.2500E-02 3.1938E-03 8.6686E-04 2.2951E-04

(b) Errors in the L ∞ ( H 1 )-norm N

|u ∗ − uh∗ | L ∞ ( H 1 )

16 32 64 128

2.7060E-01 1.3631E-01 6.8122E-02 3.4067E-02

Order

| y ∗ − yh∗ | L ∞ ( H 1 )

0.99 1.00 1.00

3.2125E-01 1.6195E-01 8.0952E-02 4.0485E-02

5.2. An example with control constraints Example 2. The interface is also a circle centered at the origin with radius r0 = 0.5. The space U ad = {u ∈ L 2 () : ua ≤ u ≤ u b }. In this example, we choose ua = −0.05, u b = 0.15 and α = 1.0. Then we set the optimal control and the associated state by

u ∗ (x) = P[−0.05,0.15] (ϕ (x)) = min {0.15, max{−0.05, ϕ (x)}} , p ∗ (x) = −αϕ (x).

(5.3)

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.19 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

19

Table 3 Grid refinement analysis for Example 1 with β − = 1, β + = 10. (a) Errors in the L ∞ ( L 2 )-norm N

u ∗ − uh∗  L ∞ ( L 2 )

16

2.4989E-03

32 64 128

6.5124E-04 1.6618E-04 4.5297E-05

Order

 y ∗ − yh∗  L ∞ ( L 2 )

Order

4.0977E-03 1.94 1.97 1.88

1.0637E-03 2.5466E-04 6.4239E-05

 p ∗ − ph∗  L ∞ ( L 2 )

Order

2.4989E-04 1.95 2.06 1.99

6.5124E-05 1.6618E-05 4.5297E-06

1.94 1.97 1.88

Order

 p ∗ − ph∗  L ∞ ( H 1 )

Order

0.99 0.94 0.91

7.8236E-03 3.9476E-03 2.0595E-03 1.0968E-03

0.99 0.94 0.91

(b) Errors in the L ∞ ( H 1 )-norm N

|u ∗ − uh∗ | L ∞ ( H 1 )

16 32 64 128

7.8236E-02 3.9476E-02 2.0595E-02 1.0968E-02

Order

| y ∗ − yh∗ | L ∞ ( H 1 )

0.99 0.94 0.91

9.2980E-02 4.6921E-02 2.4479E-02 1.3036E-02

Fig. 3. Errors of control, state, and adjoint state plotted against mesh size h for Example 1 with β − = 1 and β + = 1000. Left: L 2 -errors. Right: H 1 -errors.

We also take the optimal state as

⎧ cos(t − T )(x21 + x22 − r02 )(x1 − 1)(x1 + 1)(x2 − 1)(x2 + 1) ⎪ ⎪ ⎪ in − , ⎨− β− ∗ y = ⎪ cos(t − T )(x21 + x22 − r02 )(x1 − 1)(x1 + 1)(x2 − 1)(x2 + 1) ⎪ ⎪ ⎩− in + . β+

Then we can determine the functions f , g and yd accordingly. In this example, we also consider three cases: β − /β + = 1/1000 (large jump), β − /β + = 1000 (large jump), and β − /β + = 1/10 (moderate jump). Numerical results obtained by the IFEM are listed in Table 4, Table 5 and Table 6. These results are in agreement with the theoretical predictions in Theorem 3.15. Fig. 4 depicts the development of the L ∞ (0, T ; L 2 ) and L ∞ (0, T ; H 1 ) error for the control, state and adjoint state under uniform refinement of the mesh with β − /β + = 1/1000. From the figure, the expected order O (h2 ) in L ∞ (0, T ; L 2 ) norm for the control, state and adjoint state is observed, and the order O (h) in L ∞ (0, T ; H 1 ) norm is shown for the state and adjoint state. The exact and computed controls are shown in Fig. 5 and Fig. 6 for the cases of β − /β + = 1/1000 and β − /β + = 1000, respectively. In these figures the computed controls are obtained by the IFEM with N = 32. Examination of the figures shows that the approximate solutions nearly coincide with the exact solutions. 6. Conclusion In this paper, we developed an efficient numerical method for optimal control problems governed by parabolic interface problems. The rigorous numerical analysis of semi-discrete and fully discrete cases for the control problem of parabolic equation with interface were given. Due to low regularity of the state equation, the numerical analysis of the consideration control problem is challenging. On the other hand, the corresponding adjoint state is more regular. In contrast, optimal

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.20 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

20

Table 4 Grid refinement analysis for Example 2 with β − = 1, β + = 1000. (a) Errors in the L ∞ ( L 2 )-norm N

u ∗ − uh∗  L ∞ ( L 2 )

16 32 64 128

3.0620E-03 6.8886E-04 1.8117E-04 3.7723E-05

Order

 y ∗ − yh∗  L ∞ ( L 2 )

2.15 1.93 2.26

5.0537E-03 1.1518E-03 2.8688E-04 6.4670E-05

Order

 p ∗ − ph∗  L ∞ ( L 2 )

Order

2.14 2.00 2.15

4.0772E-03 9.2378E-04 2.3137E-04 5.1845E-05

2.14 2.00 2.16

(b) Errors in the L ∞ ( H 1 )-norm N

| y ∗ − yh∗ | L ∞ ( H 1 )

16 32 64 128

1.0310E-01 4.6824E-02 2.1996E-02 1.0746E-02

Order

 p ∗ − ph∗  L ∞ ( H 1 )

Order

1.14 1.09 1.03

8.6741E-02 3.9397E-02 1.8508E-02 9.0425E-03

1.14 1.09 1.03

Table 5 Grid refinement analysis for Example 2 with β − = 1000, β + = 1. (a) Errors in the L ∞ ( L 2 )-norm N

u ∗ − uh∗  L ∞ ( L 2 )

16

5.7044E-03

32 64 128

1.3575E-03 3.3981E-04 8.2110E-05

Order

 y ∗ − yh∗  L ∞ ( L 2 )

Order

1.2582E-02 2.07 2.00 2.05

3.2112E-03 8.1986E-04 2.1302E-04

 p ∗ − ph∗  L ∞ ( L 2 )

Order

1.0267E-02 1.97 1.97 1.94

2.6156E-03 6.6662E-04 1.7121E-04

1.97 1.97 1.96

(b) Errors in the L ∞ ( H 1 )-norm N

| y ∗ − yh∗ | L ∞ ( H 1 )

16 32 64 128

3.2125E-01 1.6196E-01 8.0950E-02 4.0485E-02

Order

 p ∗ − ph∗  L ∞ ( H 1 )

Order

0.99 1.00 1.00

2.7033E-01 1.3628E-01 6.8117E-02 3.4067E-02

0.99 1.00 1.00

Table 6 Grid refinement analysis for Example 2 with β − = 1, β + = 10. (a) Errors in the L ∞ ( L 2 )-norm N

u ∗ − uh∗  L ∞ ( L 2 )

16 32 64 128

2.4144E-03 6.4104E-04 1.6278E-04 4.3853E-05

Order

 y ∗ − yh∗  L ∞ ( L 2 )

1.91 1.98 1.89

4.1330E-03 1.0731E-03 2.5647E-04 6.4566E-05

Order

 p ∗ − ph∗  L ∞ ( L 2 )

Order

1.95 2.06 1.99

3.2951E-03 8.5571E-04 2.0639E-04 5.2584E-05

1.95 2.05 1.97

(b) Errors in the L ∞ ( H 1 )-norm N

| y ∗ − yh∗ | L ∞ ( H 1 )

16 32 64 128

9.2985E-02 4.6922E-02 2.4479E-02 1.3036E-02

Order

 p ∗ − ph∗  L ∞ ( H 1 )

Order

0.99 0.94 0.91

7.8221E-02 3.9479E-02 2.0597E-02 1.0969E-02

0.99 0.94 0.91

control problems with state constraints results in optimality systems with lower regularity of the adjoint state and more regular state. This case is under consideration. In particular, the method used a uniform mesh which is independent of the interface. Optimal convergence in space was achieved by modifying basis functions on interface elements according to interface jump conditions. We restricted this paper to the steady interface problem. In fact, an extension to the optimal control problem with the moving interface, the method developed in this paper is also easy to implement and is well suited for time-dependent moving interface problems. Acknowledgements Z. Zhang’s work and Q. Wang’s work are partially supported by the National Natural Science Foundation of China grant Nos. 11971241, 11701283. D. Liang’s work was partially supported by the Natural Sciences and Engineering Research Council

JID:APNUM AID:3644 /FLA

[m3G; v1.260; Prn:30/08/2019; 15:50] P.21 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

21

Fig. 4. Errors of control, state, and adjoint state plotted against mesh size h for Example 2 with β − = 1 and β + = 1000. Left: L ∞ ( L 2 )-errors. Right: L ∞ ( H 1 )-errors.

Fig. 5. The exact control and the computed control obtained by the IFEM with N = 32 for Example 2 with β − = 1 and β + = 1000.

Fig. 6. The exact control and the computed control obtained by the IFEM with N = 32 for Example 2 with β − = 1000 and β + = 1.

JID:APNUM AID:3644 /FLA

22

[m3G; v1.260; Prn:30/08/2019; 15:50] P.22 (1-22)

Z. Zhang et al. / Applied Numerical Mathematics ••• (••••) •••–•••

of Canada. Authors are very grateful to anonymous referees for their valuable comments and suggestions, which have greatly improved the paper. References [1] Champike Attanayake, Deepthika Senaratne, Convergence of an immersed finite element method for semilinear parabolic interface problems, Appl. Math. Sci. 5 (2011) 135–147. [2] Ivo Babuška, The finite element method for elliptic equations with discontinuous coefficients, Computing 5 (1970) 207–213. [3] Eduardo Casas, Florian Kruse, Karl Kunisch, Optimal control of semilinear parabolic equations by BV-functions, SIAM J. Control Optim. 55 (2017) 1752–1788. [4] Zhiming Chen, Jun Zou, Finite element methods and their convergence for elliptic and parabolic interface problems, Numer. Math. 79 (1998) 175–202. [5] Yanping Chen, Zuliang Lu, Error estimates of fully discrete mixed finite element methods for semilinear quadratic parabolic optimal control problem, Comput. Methods Appl. Mech. Eng. 199 (2010) 1415–1423. [6] So-Hsiang Chou, Y. Do Kwak, T. Wee Kye, Optimal convergence analysis of an immersed interface finite element method, Adv. Comput. Math. 33 (2010) 149–168. [7] Wei Gong, Michael Hinze, Zhaojie Zhou, A priori error analysis for finite element approximation of parabolic optimal control problems with pointwise control, SIAM J. Control Optim. 52 (2014) 97–119. [8] Ruchi Guo, Tao Lin, A group of immersed finite-element spaces for elliptic interface problems, IMA J. Numer. Anal. (2016). [9] Ruchi Guo, Tao Lin, Xu Zhang, Nonconforming immersed finite element spaces for elliptic interface problems, Comput. Math. Appl. 75 (2018) 2002–2016. [10] Xiaoming He, Tao Lin, Yanping Lin, Xu Zhang, Immersed finite element methods for parabolic equations with moving interface, Numer. Methods Partial Differ. Equ. 29 (2013) 619–646. [11] Michael Hinze, A variational discretization concept in control constrained optimization: the linear-quadratic case, Comput. Optim. Appl. 30 (2005) 45–61. [12] Michael Hinze, René Pinnau, Michael Ulbrich, Stefan Ulbrich, Optimization with PDE Constraints, 2008, p. 23. [13] Jianguo Huang, Jun Zou, Some new a priori estimates for second-order elliptic and parabolic interface problems, J. Differ. Equ. 184 (2002) 570–586. [14] Arbaz Khan, Chandra Shekhar Upadhyay, Marc Gerritsma, Spectral element method for parabolic interface problems, Comput. Methods Appl. Mech. Eng. 337 (2018) 66–94. [15] Zhilin Li, Tao Lin, Yanping Lin, Robert C. Rogers, An immersed finite element space and its approximation capability, Numer. Methods Partial Differ. Equ., Int. J. 20 (2004) 338–367. [16] Zhilin Li, Tao Lin, Xiaohui Wu, New Cartesian grid methods for interface problems using the finite element formulation, Numer. Math. 96 (2003) 61–98. [17] Wenbin Liu, Ningning Yan, A posteriori error estimates for optimal control problems governed by parabolic equations, Numer. Math. 93 (2003) 497–521. [18] Dominik Meidner, B. Vexler, A priori error analysis of the Petrov–Galerkin Crank–Nicolson scheme for parabolic optimal control problems, SIAM J. Control Optim. 49 (2011) 2183–2211. [19] Christian Meyer, Livia M. Susu, Optimal control of nonsmooth, semilinear parabolic equations, SIAM J. Control Optim. 55 (2017) 2206–2234. [20] Ira Neitzel, Boris Vexler, A priori error estimates for space–time finite element discretization of semilinear parabolic optimal control problems, Numer. Math. 120 (2012) 345–386. [21] Jhuma Sen Gupta, Rajen Kumar Sinha, G. Murali Mohan Reddy, Jain Jinank, New interpolation error estimates and a posteriori error analysis for linear parabolic interface problems, Numer. Methods Partial Differ. Equ. 33 (2017) 570–598. [22] Lunji Song, Shan Zhao, Symmetric interior penalty Galerkin approaches for two-dimensional parabolic interface problems with low regularity solutions, Am. J. Comput. Appl. Math. 330 (2018) 356–379. [23] Fredi Tröltzsch, Optimal Control of Partial Differential Equations: Theory, Methods, and Applications, vol. 112, American Mathematical Soc., 2010. [24] Y.J. Yuan, D. Liang, H.M. Zhu, Optimal control of groundwater pollution combined with source abatement costs and taxes, J. Comput. Sci. 20 (2017) 17–29. [25] Qian Zhang, Kazufumi Ito, Zhilin Li, Zhiyue Zhang, Immersed finite elements for optimal control problems of elliptic PDEs with interfaces, J. Comput. Phys. 298 (2015) 305–319.