Pointwise-in-time error estimates for an optimal control problem with subdiffusion constraint
Bangti Jin Buyang Li Zhi Zhou
math.NA


]Pointwise-in-time error estimates for an optimal control problem with subdiffusion constraint


Bangti Jin]Bangti Jin
Department of Computer Science, University College London, Gower Street, London, WC1E 2BT, UK.
b.jin@ucl.ac.uk


Buyang Li]Buyang Li
Department of Applied Mathematics,
The Hong Kong Polytechnic University, Kowloon, Hong Kong
buyang.li@ polyu.edu.hk

Zhi Zhou]Zhi Zhou
Department of Applied Mathematics,
The Hong Kong Polytechnic University, Kowloon, Hong Kong
zhizhou@ polyu.edu.hk




[
    [
    December 30, 2023
=====================





In this work, we present numerical analysis for a distributed optimal control problem, with box constraint on
the control, governed by a subdiffusion equation which involves a fractional
derivative of order α∈(0,1) in time. The fully discrete scheme is obtained by applying the conforming
linear Galerkin finite element method in space, L1 scheme/backward Euler convolution quadrature in time, and
the control variable by a variational type discretization. With a space mesh size h and time stepsize τ,
we establish the following order of convergence for the numerical solutions of the optimal control problem:
O(τ^min(1/2+α-ϵ,1)+h^2) in the discrete L^2(0,T;L^2(Ω)) norm and
O(τ^α-ϵ+ℓ_h^2h^2) in the discrete L^∞(0,T;L^2(Ω)) norm, with any
small ϵ>0 and ℓ_h=ln(2+1/h). The analysis relies essentially on the maximal
L^p-regularity and its discrete analogue for the subdiffusion problem. Numerical experiments are
provided to support the theoretical results.


optimal control, time-fractional diffusion, L1 scheme, convolution quadrature, pointwise-in-time
error estimate, maximal regularity.







§ INTRODUCTION

Let Ω⊂ℝ^d (d=1,2,3) be a convex polyhedral domain with a boundary ∂Ω.
Consider the distributed optimal control problem

    min_q ∈ J(u,q)=12  u - u_d  _L^2(0,T;L^2)^2 + γ2 q _L^2(0,T;L^2)^2,

subject to the following fractional-order partial differential equation

    u-Δ u= f+q,   0<t≤ T,   u(0)=0,

where T>0 is a fixed final time, γ>0 a fixed penalty parameter, Δ:H^1_0(Ω)∩
H^2(Ω)→ L^2(Ω) the Dirichlet Laplacian, f:(0,T)→ L^2(Ω) a given
source term, and u_d:(0,T)→ L^2(Ω) the target function.
The admissible set  for the control q is defined by

    ={ q∈ L^2(0,T;L^2):  a≤ q ≤ b   a.e. in  Ω×(0,T)},

with a,b∈ℝ and a< b. The notation u in (<ref>) denotes the left-sided Riemann-Liouville
fractional derivative in time t of order α∈(0,1), defined by <cit.>

    u(t)= 1/Γ(1-α)/ṭ∫_0^t(t-s)^-α u(s)ṣ.

Since u(0)=0, the Riemann-Liouville derivative u(t) coincides  with the usual Caputo derivative
<cit.>.

Further, when α=1, _0∂_t^α u(t) coincides with the first-order derivative u'(t),
and thus the model (<ref>) recovers the standard parabolic problem.

The fractional derivative u in the model (<ref>) is motivated by an ever-growing list of practical applications
related to subdiffusion processes, in which the mean square displacement grows sublinearly with time t, as opposed to linear
growth for normal diffusion. The list includes thermal diffusion in fractal media, protein transport in plasma membrane
and column experiments etc <cit.>.
The numerical analysis of the model (<ref>)
has received much attention. However, the design and analysis of numerical methods for related optimal
control problems only started to attract attention <cit.>.
The controllability of (<ref>) was discussed in <cit.> and <cit.>.
Ye and Xu <cit.> proposed space-time spectral type methods for optimal control problems
under a subdiffusion constraint, and derived error estimates by assuming sufficiently smooth state
and control variables. Antil et al <cit.> studied an optimal
control problem with space- and time-fractional models, and showed the convergence of the discrete
approximations via a compactness argument. However, no error estimate for the optimal control was
given for the time-fractional case. Zhou and Gong <cit.> proved the
well-posedness of problem (<ref>)–(<ref>) and derived L^2(0,T;L^2(Ω)) error estimates
for the spatially semidiscrete finite element method, and described a time discretization method without error
estimate. To the best of our knowledge, there is no error estimate for time discretizations of
(<ref>)–(<ref>). It is the main goal of this work to fill this gap.

This work is devoted to the error analysis of both time and space discretizations of
(<ref>)-(<ref>). The model (<ref>) is discretized by the continuous piecewise
linear Galerkin FEM in space and the L1 approximation <cit.> or backward Euler convolution quadrature
<cit.> in time, and the control q by a variational type discretization in
<cit.>. The analysis relies crucially on ℓ^p(L^2(Ω)) error estimates for the fully discrete
finite element solutions of the direct problem with nonsmooth source term. Such results are still unavailable in
the existing literature. We shall derive such estimates in Theorems <ref>
and <ref>, and use them to derive an O(τ^min(1/2+α-ϵ,1)+h^2) error
estimate in the discrete L^2(0,T;L^2(Ω)) norm for the numerical solutions of problem
(<ref>)–(<ref>), where h and τ denote the mesh size and time stepsize, respectively,
and ϵ>0 is small, cf. Theorems <ref> and <ref>. The O(τ^min(1/2+α-ϵ,1))
rate contrasts with the O(τ) rate for the parabolic counterpart
(see, e.g., <cit.>).
The lower rate for α≤1/2 is due to the limited smoothing property of
problem (<ref>), cf. Theorem <ref>. This also constitutes
the main technical challenge in the analysis. Based on the error estimate in the discrete L^2(0,T;L^2
(Ω)) norm, we further derive a pointwise-in-time error estimate O(τ^α-ϵ+ℓ_h^2h^2)
(with ℓ_h=log(2+1/h), cf. Theorems <ref> and <ref>). Our analysis
relies essentially on the maximal L^p-regularity of fractional evolution equations and its discrete
analogue <cit.>. Numerical experiments are provided to support the theoretical analysis.

The rest of the paper is organized as follows. In Section <ref>, we discuss the solution regularity
and numerical approximation for problem (<ref>). In Section <ref>, we prove error bounds on
fully discrete approximations to problem (<ref>)–(<ref>). Finally in Section <ref>,
we provide numerical experiments to support the theoretical results. Throughout, the notation c denotes a generic
constant which may differ at each occurrence, but it is always independent of the mesh size h and time stepsize τ.



§ REGULARITY THEORY AND NUMERICAL APPROXIMATION OF THE DIRECT PROBLEM

In this section, we recall preliminaries and present analysis for the direct problem

    u-Δ u= g,    0<t≤ T,
        
    u(0)=0 ,

and its adjoint problem

    z-Δ z= η,    0≤ t<T,
        
    z(T)=0,

where the fractional derivative z is defined in (<ref>) below.
In the case α∈(0,1/2], the initial condition should be understood properly: for a rough source  term g, the temporal
trace may not exist and the initial condition should be interpreted in a weak sense <cit.>.
Thus we refrain from the case of nonzero initial condition, and leave it to a future work.



 §.§ Sobolev spaces of functions vanishing at t=0

We shall use extensively Bochner-Sobolev spaces W^s,p(0,T;L^2(Ω)). For any s≥ 0 and
1≤ p< ∞, we denote by W^s,p(0,T;L^2(Ω)) the space of functions v:(0,T)→ L^2(Ω),
with the norm defined by interpolation. Equivalently, the space is equipped with the quotient norm

    v_W^s,p(0,T;L^2(Ω)):=
    inf_vv_W^s,p(ℝ;L^2(Ω)) ,

where the infimum is taken over all possible extensions v that extend v from (0,T) to ℝ.
For any 0<s< 1, one can define Sobolev–Slobodeckiǐ seminorm |·|_W^s,p(0,T;L^2(Ω)) by

    | v  |_W^s,p(0,T;L^2(Ω))^p := ∫_0^T∫_0^T v(t)-v(ξ)_L^2(Ω)^p/|t-ξ|^1+ps ṭξ̣,

and the full norm ·_W^s,p(0,T;L^2(Ω)) by

    v_W^s,p(0,T;L^2(Ω))^p = v_L^p(0,T;L^2(Ω))^p+|v|_W^s,p
    (0,T;L^2(Ω))^p .

For s>1, one can define similar seminorms and norms. Let

    C^∞_L(0,T;L^2(Ω)):={v=w|_(0,T): w∈ C^∞(ℝ;L^2(Ω)): (w)⊂ [0,∞)},

and denote by W_L^s,p(0,T;L^2(Ω)) the closure of C^∞_L(0,T;L^2(Ω)) in W^s,p(0,T;L^2(Ω)),
and by W^s,p_R(0,T;L^2(Ω)) the closure of C^∞_R(0,T;L^2(Ω)) in W^s,p(0,T;L^2(Ω)), with

    C^∞_R(0,T;L^2(Ω)):={v=w|_(0,T): w∈ C^∞(ℝ;L^2(Ω)): (w)⊂ (-∞,T]}.

By Sobolev embedding, for v∈ W_L^s,p(0,T;L^2(Ω)), there holds v^(j)(0)=0 for j=0,…,[s]-1 (with [s] being
the integral part of s>0), and also v^(j)=0 if (s-[s])p>1. For v∈ W_L^s,p(0,T;L^2(Ω)), the zero extension
of v to the left belongs to W^s,p(-∞,T;L^2(Ω)), and W_L^s,p(0,T;L^2(Ω))=W^s,p(0,T;L^2(Ω)), if
s<1/p. We abbreviate W^s,2_L(0,T;L^2(Ω)) as H^s_L(0,T;L^2(Ω)), and likewise H^s_R(0,T;L^2(Ω)) for
W^s,2_R(0,T;L^2(Ω)).

Similar to the left-sided fractional derivative u in (<ref>), the right-sided Riemann-Liouville fractional
derivative v(t) in (<ref>) is defined by

    v(t):= - 1/Γ(1-)/ṭ∫_t^T(s-t)^-v(s) ṣ .

Let any p∈(1,∞) and p'∈(1,∞) be conjugate to each other, i.e., 1/p+1/p'=1.
Since for u∈ W^α,p_L(0,T;L^2(Ω)),v∈ W^α,p'_R(0,T;L^2(Ω)), we have
_0∂_t^α u∈ L^p(0,T;L^2(Ω)), _t∂_T^α v∈ L^p'(0,T;L^2(Ω)).
Thus, there holds <cit.>:

    ∫_0^T (_0∂_t^α u(t)) v(t) ṭ = ∫_0^T u(t) (_t∂_T^α v(t)) ṭ ,
      ∀ u ∈ W^α,p_L(0,T;L^2(Ω)),  v ∈  W^α,p'_R(0,T;L^2(Ω)).




 §.§ Regularity of the direct problem

The next maximal L^p-regularity holds <cit.>, and an analogous result holds for
 (<ref>).
 
If u_0=0 and f∈ L^p(0,T;L^2(Ω)) with 1<p<∞, then (<ref>)
has a unique solution u∈ L^p(0,T; H_0^1(Ω)∩ H^2(Ω)) such that u∈ L^p(0,T;L^2(Ω)) such that

    u_L^p(0,T; H^2(Ω))
    + u_L^p(0,T;L^2(Ω))≤ cf_L^p(0,T;L^2(Ω)),

where the constant c is independent of f and T.


Now we give a regularity result. 

For g∈W^s,p( 0,T;L^2(Ω )),  s∈ [0,1/p) and p∈(1,∞), problem (<ref>) has a unique solution
u∈ W^α+s,p( 0,T;L^2(Ω )) ∩ W^s,p(0,T;H_0^1(Ω) ∩ H^2(Ω)), which satisfies

    u_W^α+s,p (0,T;L^2(Ω)) + u_W^s,p(0,T;H^2(Ω))≤ cg_W^s,p (0,T;L^2(Ω)).

Similarly, for η∈ W^s,p(0,T;L^2(Ω )),  s∈ [0,1/p) and p∈(1,∞), problem (<ref>) has
a unique solution z∈ W^α+s,p(0,T;L^2(Ω )) ∩ W^s,p(0,T;H_0^1(Ω) ∩ H^2(Ω)), which satisfies

    z_W^α+s,p (0,T;L^2(Ω)) + z_W^s,p(0,T;H^2(Ω))≤ cη_W^s,p (0,T;L^2(Ω)).



For g∈ W^s,p(0,T;L^2(Ω)),  s∈ [0,1/p) and p∈(1,∞),
extending g to be zero on Ω×[(-∞,0)∪(T,∞)] yields g∈  W^s,p(ℝ;L^2(Ω)) and

    g_W^s,p (ℝ;L^2(Ω))≤ cg_W^s,p (0,T;L^2(Ω)) .

For the zero extension to the left, we have the identity
_0∂_t^α g(t) = _-∞∂_t^α g(t) for t∈[0,T],
and there holds the relation _-∞∂_t^α g = (iξ)^αg(ξ)
<cit.>, where    denotes taking Fourier transform in
 t, and g the Fourier transform of g. Then, with ^∨ being the inverse Fourier
transform in ξ, u= [((iξ)^α-Δ)^-1g(ξ)]^∨ is a solution
of (<ref>) and

    (1+|ξ|^2)^α+s/2u(ξ)= (1+|ξ|^2)^α/2((iξ)^α-Δ)^-1(1+|ξ|^2)^s/2g(ξ) .


Let Δ be the Dirichlet Laplacian, with domain D(Δ) = H_0^1(Ω)∩ H^2(Ω).
Then the self-adjoint operator is invertible from
L^2(Ω) to D(Δ), and generates a bounded analytic semigroup <cit.>. Thus the operator

    (1+|ξ|^2)^α/2((iξ)^α-Δ)^-1

is bounded from L^2(Ω) to D(Δ) in a small neighborhood 𝒩 of ξ=0. Further, in 𝒩, the operator

    ξ/ξ̣(1+|ξ|^2)^α/2((iξ)^α-Δ)^-1
    =   α |ξ|^2/1+|ξ|^2 (1+|ξ|^2)^α/2((iξ)^α-Δ)^-1
       +α (1+|ξ|^2)^α/2((iξ)^α-Δ)^-1
     (iξ)^α((iξ)^α-Δ)^-1

is also bounded. If ξ is away from 0, then

    (1+|ξ|^2)^α/2((iξ)^α-Δ)^-1= (iξ)^-α(1+|ξ|^2)^α/2(iξ)^α((iξ)^α-Δ)^-1 ,

and in the latter case the inequalities

    |(iξ)^-α(1+|ξ|^2)^α/2|≤ c     (iξ)^α((iξ)^α-Δ)^-1≤ c

imply the boundedness of (<ref>) and (<ref>).

Since boundedness of operators is equivalent to R-boundedness of operators in L^2(Ω) (see <cit.>
for the concept of R-boundedness), the boundedness of (<ref>) and (<ref>) implies that
(<ref>) is an operator-valued Fourier multiplier <cit.>, and thus

    u_W^α+s,p(ℝ;L^2(Ω))   ≤[(1+|ξ|^2)^α+s/2u(ξ)]^∨_L^p(ℝ;L^2(Ω))
       =[(1+|ξ|^2)^α/2((iξ)^α-Δ)^-1(1+|ξ|^2)^s/2g(ξ)]^∨_L^p(ℝ;L^2(Ω))
       ≤ c[(1+|ξ|^2)^s/2g(ξ)]^∨_L^p(ℝ;L^2(Ω))≤ cg_W^s,p(ℝ;L^2(Ω)) .

This and  (<ref>) implies the desired bound on u_W^α+s,p(0,T;L^2(Ω)).
The estimate

    u_W^s,p(0,T;H^2(Ω))≤ cg_W^s,p(0,T;L^2(Ω))

follows similarly by replacing (1+|ξ|^2)^α/2((iξ)^α-Δ)^-1 with
Δ((iξ)^α-Δ)^-1 in the proof.



Below we only use the cases “p=2, s=min(1/2-ϵ,α-ϵ)” and “p>max(1/α,1/(1-α)),
s=1/p-ϵ” of Theorem <ref>. Both cases satisfy the conditions of Theorem <ref>.
A similar assertion holds for the more general case g∈ W_L^s,p(0,T;L^2(Ω)), s>0 and 1<p<∞:

    u_W^α+s,p (0,T;L^2(Ω)) + u_W^s,p(0,T;H^2(Ω))≤ cg_W_L^s,p (0,T;L^2(Ω)).

In fact, for g∈ W_L^s,p(0,T;L^2(Ω)), the zero extension of g to t≤ 0 belongs to W^s,p(-∞,T;L^2(Ω)),
which can further be boundedly extended to a function in W^s,p(ℝ;L^2(Ω)). Then the argument in
Theorem <ref> gives the desired assertion. This also indicates a certain compatibility condition for regularity pickup.




 §.§ Numerical scheme

Now we describe numerical treatment of the forward problem (<ref>), which forms the basis for the
fully discrete scheme of problem (<ref>)–(<ref>) in Section
<ref>. 
We denote by 𝒯_h a shape-regular and quasi-uniform triangulation of the
domain Ω into d-dimensional simplexes, and let

    X_h= {v_h∈ H_0^1(Ω): v_h|_K , ∀  K ∈𝒯_h}

be the finite element space consisting of continuous piecewise linear functions.
The L^2(Ω)-orthogonal projection P_h:L^2(Ω)→ X_h is defined by
(P_h ,χ_h)  =(,χ_h), for all ∈ L^2(Ω), χ_h∈ X_h,
where (·,·) denotes the L^2 inner product. Then the spatially semidiscrete Galerkin FEM for
 problem (<ref>) is to find u_h(t)∈ X_h such that u_h(0)=0 and

    ( u_h(t),χ_h) + (∇ u_h(t),∇χ_h) = (g(t),χ_h),  ∀χ_h∈ X_h,   ∀  t∈(0,T],

By introducing the discrete Laplacian Δ_h: X_h→ X_h, defined by
-(Δ_h_h,χ_h)=(∇_h,∇χ_h), for all _h, χ_h∈ X_h,
problem (<ref>) can be written as

    u_h(t) -Δ_h u_h(t) =P_hg(t),   ∀  t∈(0,T] ,  u_h(0)=0.

Similar to Theorem <ref>,  for s<1/p there holds

    u_h_W^α+s,p(0,T;L^2(Ω))+ Δ_h u_h_W^s,p(0,T;L^2(Ω))≤ c g_W^s,p(0,T;L^2(Ω)) ,

where the constant c is independent of h (following the proof of Theorem <ref>). Lemma <ref>
and Remark <ref> remain valid for the semidiscrete solution u_h, e.g.,

    Δ_hu_h_L^p(0,T;L^2(Ω))
    + u_h_L^p(0,T;L^2(Ω))≤ cf_L^p(0,T;L^2(Ω)).

These assertions will be used extensively below without explicitly referencing.

To discretize (<ref>) in time, we uniformly partition [0,T]
with grid points t_n=nτ, n,=0,1,2,…,N and a time stepsize τ=T/N≤1, and approximate
φ(t_n) by (with φ^j=φ(t_j)):

    φ^n  = τ^-α∑_j=0^nβ_n-jφ^j,

where β_j are suitable weights. We consider
two methods: L1 scheme <cit.> and backward Euler convolution quadrature (BE-CQ)
<cit.>, for which β_j are respectively given by

    β_0=1,   β_j=(j+1)^1-α-2j^1-α+(j-1)^1-α,  j=1,2,…,N,
             β_0 = 1,   β_j = -β_j-1(α-j+1)/j,  j = 1,2,…,N.

Both schemes extend the classical backward Euler scheme to the fractional case.
Then we discretize problem (<ref>) by: with g_h^n=P_hg(t_n), find U_h^n∈ X_h such that

    U_h^n -Δ_h  U_h^n  = g_h^n,   n=1,2,…,N,    U_h^0 = 0.

By <cit.> and <cit.>, we have the following error bound.


For g∈ W^1,p(0,T;L^2(Ω)), 1≤ p≤∞, let u_h and U_h^n be the semidiscrete solution and
 fully discrete solution, respectively, in (<ref>) and (<ref>). Then there holds

    U_h^n - u_h(t_n)   _L^2≤ cτ t_n^α-1g(0)_L^2(Ω)+cτ∫_0^t_n (t_n+1-s)^-1g'(s) _L^2ṣ.




Lemma <ref> slightly refines the estimates in <cit.>,
but can be proved in the same way using the following estimates in the proof of <cit.>:

    t_n^α-1≤ c (t+τ)^α-1 and  ∫_t^t_n s^α-1ṣ≤ c (t+τ)^α-1 ,



For any Banach space X, we define

    (U_h^n)_n=0^N_ℓ^p(X):=
    {   (∑_n=0^NτU_h^n_X^p)^1/p          1≤ p<∞,
       max_0≤ n≤ NU_h^n_X          p=∞ .
    .

Then the maximal ℓ^p-regularity estimate holds for (<ref>) <cit.>.

The solutions (U_h^n)_n=1^N of (<ref>) satisfy the following estimate:

    ( U_h^n)_n=1^N_ℓ^p(L^2(Ω)) + (Δ_h U_h^n)_n=1^N_ℓ^p(L^2(Ω))≤ c_p(g_h^n)_n=1^N_ℓ^p(L^2(Ω)),
      ∀  1<p<∞.





 §.§ Error estimates

Now we present ℓ^p(L^2) error estimates for g∈ W^s,p(0,T;L^2(Ω)),
0≤ s≤ 1, 1≤ p≤∞. Error analysis for
such g is unavailable in the literature. First, we give an interpolation error estimate.
This result seems standard, but we are unable to find a proof, and thus
include a proof in Appendix <ref>.


For v∈ W^s,p(0,T;L^2(Ω)), 1<p<∞ and s∈(1/p,1], let v̅^n = τ^-1∫_t_n-1^t_n v(t) ṭ,
there holds

    (v(t_n) - v̅^n)_n=1^N  _ℓ^p(L^2)≤ c τ^s v _W^s,p(0,T;L^2(Ω)) .




Our first result is an error estimate for g∈ W_L^s,p (0,T;L^2) (i.e., compatible
source). Since g may not be smooth enough in time for pointwise evaluation,
we define the averages g̅_h^n=τ^-1∫_t_n-1^t_nP_h g(s)ṣ, and
consider a variant of the scheme (<ref>) for problem (<ref>):
find U_h^n ∈ X_h such that

    U_h^n -Δ_h U_h^n = g̅_h^n ,   n = 1,…,N,   U_h^0=0.



For g∈ W^s,p_L(0,T; L^2), 1  < p<∞ and s∈[0,1], let u_h and U_h^n be the solutions of problems (<ref>)
and (<ref>), respectively, and u̅_h^n:=τ^-1∫_t_n-1^t_nu_h(s)ṣ. Then there holds

    (U_h^n - u̅_h^n)_n=1^N  _ℓ^p(L^2)≤ cτ^ s g _W^s,p_L(0,T;L^2).
 



By Hölder's inequality and (<ref>) (with s=0), we have

    τ∑_n=1^Nu̅_h^n_L^2^p     = τ∑_n=1^N τ^-1∫_t_n-1^t_nu_h(s)ṣ_L^2^p
      ≤τ^1-p∑_n=1^N(∫_t_n-1^t_nu_h(s)_L^2(Ω)ṣ)^p
       ≤∫_0^Tu_h(s)_L^2^pṣ≤ cg_L^p(0,T;L^2)^p .

Similarly, by applying Lemma <ref> to (<ref>) and the L^2(Ω) stability of P_h, we have

    (U_h^n)_n=1^N_ℓ^p(L^2(Ω))≤ c   (g̅_h^n)_n=1^N  _ℓ^p(L^2(Ω))≤ c  g _L^p(0,T;L^2) .

This and the triangle inequality show the assertion for s=0.

Next we consider g∈ W_L^1,p(0,T;L^2), and resort to (<ref>). Since g(0)=0, by Lemma
<ref>,

    U_h^n - u_h(t_n)   _L^2≤ cτ∫_0^t_n (t_n+τ-s)^-1 g'(s) _L^2ṣ.

This directly implies

    ( U_h^n -u_h(t_n))_n=1^N  _ℓ^∞(L^2)   ≤  c τ g _W_L^1,∞(0,T;L^2).

Further, let ψ(s)=τ∑_n=1^N(t_n+τ-s)^α-1χ_[0,t_n], where χ_S denotes the
characteristic function of a set S. Then clearly, we have

    sup_s∈ [0,T]ψ(s)=ψ(T)=τ∑_n=1^N(nτ)^α-1≤∫_0^Ts^α-1ṣ = α^-1T^α≤ c_T .

Therefore,

    (U_h^n -u_h(t_n))_n=1^N  _ℓ^1(L^2)   ≤ cτ^2 ∑_n=1^N ∫_0^t_n (t_n+τ-s)^-1 g'(s) _L^2 ṣ
       = cτ∫_0^Tψ(s)g'(s)_L^2(Ω)ṣ≤ c_Tτ g _W^1,1(0,T;L^2).

Then (<ref>), (<ref>) and Riesz-Thorin
interpolation theorem <cit.>, yield for any 1<p<∞

    (U_h^n -u_h(t_n))_n=1^N  _ℓ^p(L^2)≤  c τ g _W_L^1,p(0,T;L^2).

Since U_h^n-U_h^n satisfies the discrete scheme, cf. (<ref>)
and (<ref>), Lemmas <ref> and <ref> imply

    (U_h^n - U_h^n)_n=1^N_ℓ^p(L^2(Ω))≤  c   (g̅_h^n - P_hg(t_n))_n=1^N_ℓ^p(L^2(Ω))≤ cτg_W^1,p(0,T;L^2(Ω)).

Further, by Lemma <ref> and Remark <ref>, we have

    (u_h(t_n) - u̅_h^n)_n=1^N_ℓ^p(L^2)≤ cτu_h_W^1,p(0,T;L^2)≤ cτg_W_L^1,p(0,T;L^2) .

The last three estimates show the assertion for s=1.
The case 0<s<1 follows by interpolation.


In Theorem <ref>, we compare the numerical solution U_h^n to (<ref>) with the time-averaged
solution u̅_h^n, instead of u_h(t_n). This is due to possible insufficient temporal regularity of u_h:
it is unclear how to define u_h(t_n) for t_n∈(0,T] for g∈ W_L^s,p(0,T; L^2) with s+α<1/p. For s∈(1/p,1],
W_L^s,p(0,T; L^2) ↪ C([0,T] ;L^2) and so Theorem <ref> requires the condition
g(0)=0. Such a compatibility condition at t=0 is not necessarily satisfied by (<ref>)–(<ref>). Hence, we state an error estimate
below for a smooth but incompatible source g∈ W^s,p(0,T; L^2).


For g∈ W^s,p(0,T; L^2) with p∈(1,∞) and s∈(1/p,1), let u_h and U_h^n be the
solutions of (<ref>) and (<ref>), respectively. Then there holds

    (U_h^n -  u_h(t_n))_n=1^N  _ℓ^p(L^2)≤ c τ^min(1/p+α,s) g _W^s,p(0,T; L^2) .

Moreover, if p>1/α is so large that α∈(0,1/p'), then

    (U_h^n -  u_h(t_n))_n=1^N  _ℓ^∞(L^2)≤ c τ^α g _W^1/p+α,p(0,T; L^2).



For g(x,t)≡ g(x), which belongs to W^1,p(0,T;L^2(Ω)), by Lemma
<ref> we have

    (U_h^n - u_h(t_n))_n=1^N _ℓ^p(L^2)^p
          =  τ∑_n=1^N   U_h^n - u_h(t_n)  _L^2^p
      ≤  c τ^p+1  g _L^2^p ∑_n=1^N  t_n^p(α-1)
       ≤ c τ^pα+1  g _L^2^p  + c τ^p  g _L^2^p ∫_τ^T  t^p(α-1) ṭ.

Now for s<1, there holds

    ∫_τ^Tt^p(α-1) ṭ≤{   c τ^pα+1         α∈ (0,1/p') 
       c τ^p         α∈ (1/p' ,1)
       c τ^p (1+ln (T/τ))          α = 1/p'
                       .
         ≤{   c τ^p(1/p+α)         α∈ (0,1/p') 
       c τ^ps         α∈ (1/p' ,1)
       c τ^ps         α = 1/p'
                       .

which together with the preceding estimate implies

    (U_h^n - u_h(t_n))_n=1^N _ℓ^p(L^2(Ω))≤ cτ^min(1/p+α,s) g _W^s,p(0,T; L^2(Ω)) .

For g∈ W^s,p(0,T; L^2), by Sobolev embedding, g(0) exists,
and in the splitting g(t)=g(0) + (g(t)-g(0)), there holds
g-g(0)_W_L^s,p(0,T;L^2(Ω))≤ cg_W^s,p(0,T;L^2(Ω)).
Let v_h be the semidiscrete solution for the source g(t)-g(0), and v̅_h^n=∫_t_n-1^t_nv_h(t)ṭ.
Since g(t)-g(0)∈ W_L^s,p(0,T; L^2), by Theorem <ref>, the corresponding fully discrete solution
V_h^n by (<ref>) satisfies

    (V_h^n- v̅_h^n)_n=1^N_ℓ^p(L^2)≤ cτ^sg-g(0)_W_L^s,p(0,T; L^2).

Further, by Lemma <ref> and (<ref>), we have

    (v_h(t_n)-v̅_h^n)_n=1^N_ℓ^p(L^2)≤ c τ^sv_h_W^s,p(0,T;L^2(Ω))   ≤ cτ^sg-g(0)_W_L^s,p(0,T; L^2).

Similarly, for the fully discrete solution V_h^n for the source g(t)-g(0) by (<ref>),
from Lemmas <ref> and <ref>, we deduce

    (V_h^n-V_h^n)_n=1^N_ℓ^p(L^2(Ω))   ≤ c (P_hg(t_n)-g̅_h^n)_n=1^N_ℓ^p(L^2(Ω))
       ≤ cτ^sg_W^s,p(0,T;L^2(Ω)).

These estimates together with the triangle inequality give (<ref>).

















Finally, (<ref>) follows by the inverse inequality in time and (<ref>) with s=1/p+α, i.e.,

    (U_h^n -  u_h(t_n))_n=1^N  _ℓ^∞(L^2)   ≤ c τ^-1/p (U_h^n -  u_h(t_n))_n=1^N  _ℓ^p(L^2)
       ≤ cτ^αg_W^1/p+α,p(0,T; L^2) .

This completes the proof of the theorem.



The estimate (<ref>) is not sharp since,
according to the proof, the restriction s<1 is only needed for α = 1/p'. Nonetheless,
it is sufficient for the error analysis in Section <ref>.



§ THE OPTIMAL CONTROL PROBLEM AND ITS NUMERICAL APPROXIMATION

In this section, we develop a numerical scheme for problem (<ref>)–(<ref>),
and derive error bounds for the spatial and temporal discretizations.



 §.§ The continuous problem

The first-order optimality condition of  (<ref>)-(<ref>) was given in <cit.>.

Problem (<ref>)-(<ref>) admits a unique solution q∈.
There exist a state u∈ L^2(0,T; D(Δ))∩ H_L^(0,T; L^2) and an adjoint
z∈ L^2(0,T; D(Δ))∩ H_R^(0,T; L^2) such that (u,z,q) satisfies

    u-Δ u   = f+q,    0<t≤ T,     u(0)=0 ,
    
       _t∂_T^α z-Δ z   = u-u_d,     0≤ t<T,     z(T)=0,
    
       (γ q + z,   v-q)_L^2(0,T;L^2(Ω))≥ 0,    ∀ v∈.

where (·,·)_L^2(0,T;L^2(Ω)) denotes the L^2(0,T;L^2(Ω)) inner product.


Let P_ be the nonlinear pointwise projection operator defined by

    P_ (q) = max(a,min(q,b)).

It is bounded on W^s,p(0,T; L^2(Ω)) for 0≤ s≤ 1 and 1≤ p≤∞


    P_u_W^s,p(0,T;L^2(Ω))≤ cu_W^s,p(0,T; L^2(Ω)).

This estimate holds trivially for s=0 and s=1 (see
<cit.>), and the case 0<s<1 follows by interpolation.
Then (<ref>) is equivalent to the complementarity condition 

    q= P_(-γ^-1 z).


Now we give higher regularity of the triple (u,z,q).

For any s∈(0,1/2), let u_d∈ H^s(0,T;L^2(Ω)) and f∈ H^s(0,T;L^2(Ω)).
Then the solution (u,z,q) of problem (<ref>)–(<ref>) satisfies the following estimate

    q_H^min(1,α+s) (0,T;L^2(Ω)) +u_H^α+s(0,T;L^2(Ω))+
      z_H^α+s(0,T;L^2(Ω))≤ c.




Let r=min(1,α+s). By (<ref>) and (<ref>), we have

    q_H^r (0,T;L^2(Ω))≤ cz_H^r (0,T;L^2(Ω))≤ cz _H^α+s (0,T;L^2(Ω)).

Applying Theorem <ref> to (<ref>) yields

    z_H^α+s (0,T;L^2(Ω))   ≤ c(u_H^s(0,T;L^2(Ω))+u_d_H^s(0,T;L^2(Ω)))
      ≤ c u_H^s (0,T;L^2(Ω))+c.

Similarly, applying Theorem <ref> to (<ref>) gives

    u_H^α+s (0,T;L^2(Ω))   ≤ c(f_H^s (0,T;L^2(Ω))+q_H^s(0,T;L^2(Ω)))
       ≤ c+cq_H^s(0,T;L^2(Ω)) .

The last three estimates together imply

    q_H^r (0,T;L^2(Ω))   ≤ c+cq_H^s(0,T;L^2(Ω))≤ c+c_ϵ'q_L^2(0,T;L^2(Ω))  +ϵ'q_H^r (0,T;L^2(Ω)) ,

where the last step is due to the interpolation inequality <cit.>

    q_H^s(0,T;L^2(Ω))≤ c_ϵ'q_L^2(0,T;L^2(Ω))  +ϵ'q_H^r (0,T;L^2(Ω)).

By choosing a small ϵ'>0 and the pointwise boundedness of q,
cf. (<ref>), we obtain

    q_H^r (0,T;L^2(Ω))≤  c+c_ϵ'q_L^2(0,T;L^2(Ω))≤ c.

This shows the bound on q.  (<ref>) and (<ref>) give
the bound on u, and that of z follows similarly.


Next, we give an improved stability estimate on q.


Let p>1/α be sufficiently large so that α∈(0,1/p'), u_d∈
W^α,p(0,T;L^2(Ω)) and f∈ L^p(0,T;L^2(Ω)).
Then the optimal control q satisfies:

    q_W^1/p+α-ϵ,p(0,T;L^2(Ω))≤ c,

where the constant c depends on u_d_W^α,p(0,T;L^2(Ω)) and f_L^p(0,T;L^2(Ω)).


The condition α∈(0,1/p') implies r:=1/p+α-ϵ<1. Thus
(<ref>) and Theorem <ref> (with s=r-α) imply

    q_W^r,p(0,T;L^2(Ω))≤    c z_W^r,p(0,T;L^2(Ω))≤  cu-u_d_W^r-α,p(0,T;L^2(Ω)) ,

Since p>1/α, r-α=1/p-ϵ<α and thus Theorem <ref> (with s=0)
and (<ref>) give

    u-u_d_W^r-α,p(0,T;L^2(Ω))≤  cu-u_d_W^α,p(0,T;L^2(Ω))≤  cf+q_L^p(0,T;L^2(Ω))+c≤ c,



The last two estimates together imply the desired result.




 §.§ Spatially semidiscrete scheme

Now we give a spatially semidiscrete scheme for problem (<ref>)–(<ref>): find q_h∈ such that

    min_q_h∈ J(u_h, q_h) = 12 u_h - u_d  _L^2(0,T;L^2)^2 + γ2 q_h _L^2(0,T;L^2)^2,

subject to the semidiscrete problem

    u_h-Δ_h u_h= P_h(f+q_h),    0<t≤ T,     u_h(0)=0.

Similar to Theorem <ref>, problem (<ref>)-(<ref>) admits a unique solution
q_h∈. The first-order optimality system reads:

    u_h-Δ_h u_h   = P_h(f+q_h),   0<t≤ T,     u_h(0)=0,
    
       _t∂_T^α z_h-Δ_h z_h   = P_h(u_h-u_d),   0≤ t<T,     z_h(T)=0, 
    
         (γ q_h + z_h,   v-q_h)_L^2(0,T;L^2(Ω))≥ 0,   ∀ v∈.

The variational inequality (<ref>) is equivalent to

    q_h=P_(-γ^-1 z_h).


For the approximation (<ref>)–(<ref>), see <cit.> for an error estimate.

For f, u_d∈ L^2(0,T;L^2(Ω)), let (u,z,q) and (u_h,z_h,q_h) be the solutions of
problems (<ref>)–(<ref>) and (<ref>)–(<ref>), respectively.
Then there hold

    u- u_h  _L^2(0,T;L^2(Ω)) +  z-z_h  _L^2(0,T;L^2(Ω)) +  q-q_h_L^2(0,T;L^2(Ω))   ≤ c  h^2,
    ∇( u- u_h) _L^2(0,T;L^2(Ω)) + ∇( z-z_h  )_L^2(0,T;L^2(Ω))   ≤ ch.



Next, we present the regularity of the semidiscrete solution (u_h,z_h,q_h). The proof is similar
to the continuous case in Lemmas <ref> and <ref> and hence omitted.

Let s∈(0,1/2), u_d∈ H^s(0,T;L^2(Ω)) and f∈ H^s(0,T;L^2(Ω)).
Then the solution (u_h,z_h,q_h) of problem (<ref>)–(<ref>) satisfies the following estimate:

    q_h _H^min(1,α+s) (0,T;L^2(Ω)) +u_h_H^α+s(0,T;L^2(Ω))+
      z_h_H^α+s(0,T;L^2(Ω))≤ c.

Further, for u_d∈ W^α,p(0,T;L^2(Ω)), f∈ L^p(0,T;L^2(Ω)), with p>1/α and α∈(0,1/p'), there holds

    q_h_W^1/p+α-ϵ,p(0,T;L^2(Ω))≤ c.



Last, we derive a pointwise-in-time error estimate.

For f,  u_d∈ H^1(0,T;L^2(Ω)), let (u,z,q) and (u_h,z_h,q_h) be the solutions
of problems (<ref>)–(<ref>) and (<ref>)–(<ref>), respectively.
Then there holds

    u- u_h  _L^∞(0,T;L^2(Ω)) +  z-z_h  _L^∞(0,T;L^2(Ω)) +  q-q_h_L^∞(0,T;L^2(Ω))≤ c ℓ_h^2 h^2

with ℓ_h=log(2+1/h), where the constant c depends on f_H^1(0,T;L^2(Ω)) and u_d_H^1(0,T;L^2(Ω)).


We employ the splitting
u - u_h = (u - u_h(q)) + (u_h(q) - u_h):=ϱ+ϑ,
where u_h(q)∈ X_h solves

    ∂_t^α u_h(q) -Δ_h u_h(q) = P_h(f + q),   0<t≤ T,  with   u_h(q)(0)=0.

Then u_h(q) is the semidiscrete solution of (<ref>) with g=f+q, and ϱ
is the FEM error for the direct problem. By <cit.> and Lemma <ref>, we have

    ϱ_L^∞(0,T;L^2)≤ c ℓ_h^2  h^2f+q_L^∞(0,T;L^2(Ω))≤ cℓ_h^2h^2.

Since ϑ satisfies ∂_t^αϑ -Δ_h ϑ= P_h (q-q_h), for 0<t≤ T with ϑ(0)=0,
(<ref>), L^2(Ω)-stability of P_h, the conditions (<ref>)
and (<ref>), and the pointwise contractivity of P_ imply

    _0∂_t^αϑ_L^p(0,T;L^2(Ω))≤ P_h (q -  q_h)  _L^p(0,T;L^2(Ω))
    ≤    c   q - q_h  _L^p(0,T;L^2(Ω))≤ c  z-z_h _L^p(0,T;L^2(Ω))  .

Next, it follows from (<ref>), (<ref>) and the identity P_hΔ=Δ_hR_h
(with R_h:H^1(Ω)→ X_h being Ritz projection) that w_h:=P_hz-z_h satisfies w_h(T)=0

    _t∂_T^α w_h-Δ_h w_h= P_hu-u_h
    -Δ_h(P_hz-R_hz),   0≤ t<T,

and thus

    _t∂_T^αΔ_h^-1w_h-Δ_h Δ_h^-1w_h   = Δ_h^-1(P_hu-u_h) - (P_hz-R_hz) .

The maximal L^p regularity (<ref>) and triangle inequality imply

    w_h_L^p(0,T;L^2(Ω))   ≤ cΔ_h^-1(P_hu-u_h) - (P_hz-R_hz)_L^p(0,T;L^2(Ω))
       ≤ cP_hu-u_h_L^p(0,T;L^2(Ω))
         +cP_hz-R_hz_L^p(0,T;L^2(Ω)).

The L^2(Ω)-stability of P_h and triangle inequality yield

    P_hu-u_h_L^p(0,T;L^2(Ω))≤ cu-u_h_L^p(0,T;L^2(Ω))≤ c(ϑ_L^p(0,T;L^2(Ω))+ϱ_L^p(0,T;L^2(Ω))),

and Theorem <ref> (with s=0) and lemma <ref> give

    P_hz-R_hz_L^p(0,T;L^2(Ω))≤ cz_L^p(0,T;H^2(Ω)) h^2 ≤ cu-u_d_L^p(0,T;L^2(Ω)) h^2≤ ch^2.












The last three estimates and (<ref>) yield

    w_h_L^p(0,T;L^2(Ω))≤ cϑ_L^p(0,T;L^2(Ω))+cℓ_h^2h^2.

Thus repeating the preceding argument yields

    z-z_h_L^p(0,T;L^2(Ω))   ≤z-P_hz_L^p(0,T;L^2(Ω))+w_h_L^p(0,T;L^2(Ω))
       ≤ cϑ_L^p(0,T;L^2(Ω))+cℓ_h^2h^2.

Substituting it into (<ref>) and by Sobolev embedding W^α,p
(0,T;L^2(Ω))↪ L^p_α(0,T;L^2(Ω)), with the critical exponent p_α=p/(1-pα)
if pα<1, and p_α=∞ if pα>1:

    ϑ_L^p_α(0,T;L^2(Ω))≤ c_0∂_t^αϑ_L^p(0,T;L^2(Ω))≤ cϑ_L^p(0,T;L^2(Ω))+cℓ_h^2h^2 .

A finite number of repeated applications of this inequality yields

    ϑ_L^∞(0,T;L^2(Ω))≤ c ϑ_L^2(0,T;L^2(Ω))+ cℓ_h^2 h^2  ≤ c ℓ_h^2 h^2,

where we have used the fact that, by maximal L^p regularity (<ref>) and Theorem <ref>,

    ϑ_L^2(0,T;L^2(Ω))=u_h(q)-u_h_L^2(0,T;L^2(Ω))≤ cq-q_h_L^2(0,T;L^2(Ω))≤ ch^2.

This gives the desired bound on u-u_h_L^∞(0,T;L^2). The bounds on z - z_h_L^∞(0,T;L^2)
and  q  -q_h_L^∞(0,T;L^2) follow similarly by the contraction property of P_U_ad.




 §.§ Fully discrete scheme

Now we turn to the fully discrete approximation of (<ref>)–(<ref>), with L1 scheme
or BE-CQ time stepping. First, we define a discrete admissible set

    U_ad^τ = {Q⃗_h=(Q_h^n-1)_n=1^N: a≤ Q^n-1≤ b,   n=1,2,...,N },

and consider the following fully discrete problem:

    min_Q⃗_h ∈ U_ad^ττ/2∑_n=1^N ( U_h^n - u_d^n  _L^2^2 + γ Q_h^n-1_L^2^2 ),

subject to the fully discrete problem

    U_h^n-Δ_h U_h^n = f_h^n+P_h Q_h^n-1,     n=1,2,...,N,   with U_h^0 =0,

with u_d^n=u_d(t_n) and f_h^n=P_hf(t_n).
Let ∂̅_τ^αφ^n be the L1/BE-CQ approximation of _t∂_T^αφ (t_n):

    ∂̅_τ^αφ^N-n
      =τ^-α∑_j=0^nβ_n-jφ^N-j .

Then the fully discrete problem is to find (U_h^n, Z_h^n, Q_h^n) such that

    U_h^n-Δ_h U_h^n    = f_h^n+P_h Q_h^n-1,         n=1,2,...,N,   with  U_h^0 =0 ,
    
    ∂̅_τ^α Z_h^n-1-Δ_h Z_h^n-1   = U_h^n-P_h u_d^n ,       n=1,2,...,N,  with  Z_h^N =0 ,
    
    
       (γ Q_h^n-1  + Z_h^n-1  , v-Q_h^n-1 )    ≥ 0,       ∀  v∈ L^2    a≤ v≤ b.

Similar to (<ref>), (<ref>) can be rewritten as

    Q_h^n-1=P_U_ad(-γ^-1Z_h^n-1),   n=1,2,…,N .


To simplify the notation, we define a discrete L^2(0,T;L^2(Ω)) inner product [·,·]_τ by

    [ v, w]_τ = τ∑_n=1^N (v_n,w_n)   ∀  v=(v_n)_n=1^N,
     w=(w_n)_n=1^N∈ L^2(Ω)^N ,

and denote by ·_τ the induced norm.
Let v=( v^n_h)_n=1^N ∈ L^2(Ω)^N and ∂̅_τ^α w=(∂̅_τ^α w^n-1_h)_n=1^N∈ L^2(Ω)^N.
Then the discrete integration by parts formula holds <cit.>

    [ v,  w]_τ
    =[ v,∂̅_τ^α w]_τ  ∀  v , w∈ L^2(Ω)^N  .

Thus, ∂̅_τ^α is the adjoint to  with respect to [·,·]_τ. Let

    U_h=(U_h^n)_n=1^N,        Z_h=(Z_h^n-1)_n=1^N,        Q_h=(Q_h^n-1)_n=1^N ,
        u_h=(u_h(t_n))_n=1^N,        z_h=(z_h(t_n-1))_n=1^N,        q_h=(q_h(t_n-1))_n=1^N .


Next we introduce two auxiliary problems. Let 𝐔_h(q_h)=(U_h^n(q_h)) )_n=1^N∈ X_h^N solve

    U_h^n(q_h)  -Δ_hU_h^n(q_h)= f_h^n + q_h(t_n-1) ,   n=1,…,N,  U_h^0(q_h)=0 .

By Lemma <ref>, the pointwise evaluation q_h(t_n) does make sense, and thus problem
(<ref>) is well defined. For any v_h=(v_h^n)_n=1^N, let Z_h( v_h)=
(Z_h^n-1( v_h))_n=1^N∈ X_h^N solve

    ∂̅_τ^α Z_h^n-1( v_h) -Δ_h Z_h^n-1( v_h)
    = v_h^n - P_h u_d^n ,   n=1,2,…,N,  with  Z_h^N( v_h) =0 .


The rest of this part is devoted to error analysis. First, we bound q_h- Q_h_τ.

For 𝐐_h, 𝐪_h, Z_h and Z_h( U_h(q_h)) defined as above, there holds

    γ Q_h- q_h_τ ^2
      ≤  [ q_h- Q_h , Z_h( U_h(q_h))-  z_h]_τ.



It follows from  (<ref>) and (<ref>), similarly from (<ref>) and (<ref>), that

    ( -Δ_h) ( U_h(q_h)- U_h) =   q_h-  Q_h    
       (∂̅_τ^α -Δ_h) ( Z_h( v_h)- Z_h) =  v_h -  U_h.

Together with (<ref>), these identities imply

    [ q_h- Q_h ,  Z_h- Z_h( U_h(q_h))]_τ   =[( - Δ_h)( U_h (q_h)- U_h)  ,  Z_h- Z_h( U_h(q_h))]_τ
       = [ U_h (q_h)- U_h  , (∂̅_τ^α - Δ_h)( Z_h- Z_h( U_h(q_h)))]_τ
       = -  U_h (q_h)- U_h _τ^2≤0.

Next, since (<ref>) holds pointwise in time, i.e.,
q_h(t_n-1) = P_U_ad(-γ^-1z_h(t_n-1)), we have

    (q_h(t_n-1)+γ^-1z_h(t_n-1),χ-q_h(t_n-1))≥ 0,  ∀ χ∈ L^2(Ω)        a≤χ≤ b.

Upon setting v=q_h(t_n-1) in (<ref>) and χ=Q_h^n-1 in (<ref>), we deduce

    γ Q_h- q_h_τ ^2    
      =γ[ Q_h - q_h,  Q_h]_τ  - γ[ Q_h - q_h,  q_h]_τ
       ≤ [ q_h- Q_h ,  Z_h]_τ - [ q_h- Q_h ,  z_h]_τ
       = [ q_h- Q_h ,  Z_h- Z_h( U_h(q_h))]_τ + [ q_h- Q_h , Z_h( U_h(q_h))-  z_h]_τ.

Now invoking (<ref>) completes the proof of the lemma.


The next result gives an error estimate for the approximate state 𝐔_h(q_h).

Let f,u_d∈ H^1(0,T;L^2(Ω)). For any ϵ∈(0,min(1/2,α)), there holds

    U_h(q_h) -  u_h_τ≤ cτ^1/2+min(1/2,α-ϵ) .




By the triangle inequality, we have

    U_h(q_h) -  u_h_τ≤ U_h(q_h)- U_h(q_h)_τ +  U_h(q_h)- u_h_τ,

where 𝐔_h(q_h)=(U_h^n(q_h))_n=1^N is the solution to

    U_h^n(q_h)-Δ_h^nU_h^n(q_h)  = f_h^n + q̃_h^n,    n =1,2,…,N  U_h^0(q_h)=0 ,

with q̃_h^n=P_hq_h(t_n) (and q̃_h=(q̃_h^n)_n=1^N). That is, U_h^n(q_h)
is the fully discrete solution of problem (<ref>) with
g=f+q_h. By Lemmas <ref> and <ref>, we have

    U_h(q_h)- U_h(q_h)_τ≤ c q_h-q̃_h_τ≤ cτ^min(1/2+α-ϵ,1)q_h_H^1/2+α-ϵ(0,T;L^2(Ω)).

Further, Theorem <ref> (with s=min(1,1/2+α-ϵ)∈(1/2,1)) implies

    U_h(q_h) -  u_h_τ   ≤ c P_hf+q_h _ H^s(0,T;L^2(Ω))τ^s
       ≤ c(P_hf_H^s(0,T;L^2(Ω))+q_h _H^s(0,T;L^2(Ω))) τ^s .

The last two estimates and Lemma <ref> (with s=1/2-ϵ) yield the desired assertion.


Now we can give an ℓ^2(L^2(Ω)) error estimate for the approximation (U_h^n,Z_h^n,Q_h^n).

For f∈  H^1 (0,T;L^2) and u_d∈  H^1(0,T;L^2(Ω)),
(u_h,z_h,q_h) and (U_h^n,Z_h^n,Q_h^n) be the solutions of problems (<ref>)-(<ref>)
and (<ref>)-(<ref>), respectively. Then there holds for any small ϵ>0

    u_h-  U_h  _τ +  z_h-  Z_h  _τ +  q_h-  Q_h  _τ≤ c τ ^1/2+min(1/2,α-ϵ) ,

where the constant c depends on f_H^1(0,T;L^2(Ω))
and u_d_H^1(0,T;L^2(Ω)).


By Lemma <ref> and the triangle inequality, we deduce

    Q_h- q_h_τ   ≤ c  Z_h( U_h(q_h))-  Z_h( 𝐮_h)_τ+ c  Z_h(𝐮_h)-  z_h_τ.

It suffices to bound the two terms on the right hand side.
Lemmas <ref> and <ref> imply

    Z_h( U_h(q_h))-  Z_h( 𝐮_h)_τ≤ c U_h(q_h) -  u_h _τ≤ cτ^r.

with r=1/2+min(1/2,α-ϵ).
Further, since 𝐙_h(𝐮_h) is a fully discrete approximation to z_h(u_h), by
Theorem <ref> (with s=r) and Lemma <ref>, we have

    Z_h( u_h)-  z_h_τ≤ c u_h-P_hu_d _H^r (0,T;L^2(Ω))τ^r ≤ cτ ^r.

Thus, we obtain the estimate Q_h- q_h_τ≤ c τ^r.
Next, by Lemmas <ref> and <ref>, we deduce

    U_h -  u_h _τ   ≤ U_h - U_h ( q_h)_τ+ U_h ( q_h) -  u_h _τ
       ≤ c Q_h -  q_h _τ+ U_h ( q_h) -  u_h _τ≤ c τ^r.

Similarly, Z_h -  z_h _τ can be bounded by

    Z_h -  z_h _τ   ≤ Z_h - Z_h ( u_h)_τ+ Z_h ( u_h) -  z_h _τ
       ≤ c U_h -  u_h _τ
         + Z_h ( u_h) -  z_h _τ≤ cτ^r.

This completes the proof of Theorem <ref>.


Last, we give a pointwise-in-time error estimate for the approximation (U_h^n,Q_h^n,Z_h^n).

For f,u_d∈ W^1,p (0,T;L^2)∩ H^1(0,T;L^2(Ω)), p>1/α with α∈(0,1/p'),
let (u_h,z_h,q_h) and (U_h^n,Z_h^n,Q_h^n) be the solutions of problems (<ref>)-(<ref>)
and (<ref>)-(<ref>), respectively. Then there holds for any small ϵ>0

    max_1≤ n≤ N(u^n_h- U^n_h _L^2(Ω) + z^n-1_h- Z^n-1_h _L^2(Ω)
    + q^n-1_h- Q^n-1_h _L^2(Ω))  ≤ c τ^α-ϵ ,

where the constant c depends on f_W^1,min(p,2)(0,T;L^2(Ω)) and u_d_W^1,min(p,2)(0,T;L^2(Ω)).


It follows from (<ref>) and (<ref>) that U_h^0-U_h^0(q_h) =0 and

    (U_h^n-U_h^n(q_h)) -Δ_h(U_h^n-U_h^n(q_h)) =  Q_h^n-1-q_h(t_n-1),   n=1,…,N .

By Lemma <ref> and the inverse inequality (in time), we obtain for any 1/α< p_1<∞

    (U_h^n-U_h^n(q_h))_n=1^N_ℓ^p_1(L^2(Ω))   ≤ c (Q_h^n-1-q_h(t_n-1))_n=1^N_ℓ^p_1(L^2(Ω))
       ≤   cτ^min(0,1/p_1-1/2)(Q_h^n-1-q_h(t_n-1))_n=1^N_ℓ^2(L^2(Ω)).

This and Theorem <ref> imply

    (U_h^n-U_h^n(q_h))_n=1^N_ℓ^p_1(L^2(Ω))≤  cτ^min(1/p_1,1/2) +min(1/2,α-ϵ).

By choosing p_1>1/α sufficiently close to 1/α and
discrete embedding <cit.>,

    (U_h^n-U_h^n(q_h))_n=1^N_ℓ^∞(L^2(Ω))   ≤ c (U_h^n-U_h^n(q_h))_n=1^N_ℓ^p_1(L^2(Ω))
       ≤ cτ^min(1/p_1,1/2)  +min(1/2,α-ϵ)≤ cτ^α-ϵ  ,

where the last inequality follows from the inequality min(1/p_1,1/2)  +min(1/2,α-ϵ)≥α-ϵ, due to the choice
of p_1. Further, by the definition of U_h^n(q_h) in  (<ref>), choosing p_2>1/α
sufficiently large so that α∈(0,1/p_2') and applying (<ref>) and Lemma <ref>, we get

    u_h(t_n)-U_h^n(q_h)_L^2(Ω)≤      cτ^α-ϵf+q_h_W^1/p_2+α-ϵ,p_2(0,T;L^2(Ω))
    ≤       cτ^α-ϵ(q_h_W^1/p_2+α-ϵ,p_2(0,T;L^2(Ω)) +c) ≤ cτ^α-ϵ.

Last, by choosing p_3>1/α sufficiently close to 1/α, Lemmas <ref>, <ref>,
and <ref>, and discrete embedding <cit.>, we obtain

    (U_h^n(q_h)-U_h^n(q_h))_n=1^N_ℓ^∞(L^2(Ω))   ≤ c (U_h^n(q_h)-U_h^n(q_h))_n=1^N_ℓ^p_3(L^2(Ω))
       ≤ c(q_h(t_n-1)-q_h(t_n))_n=1^N_ℓ^p_3(L^2(Ω))≤ cτ^α-ϵ.

The last three estimates yield the desired bound on u_h(t_n)-U_h^n_L^2(Ω).
The bound on z_h(t_n-1)-Z_h^n-1_L^2(Ω) follows similarly, and that
on q_h(t_n-1)-Q_h^n-1_L^2(Ω) by the contraction property of P_.




§ NUMERICAL RESULTS AND DISCUSSIONS

Now we present numerical experiments to illustrate the theoretical findings. We perform experiments on the unit
interval Ω=(0,1). The domain Ω is divided into M equally spaced subintervals with a mesh size
h = 1/M. To discretize the fractional derivatives _0∂_t^α u and _t∂_T^α z, we fix the
time stepsize τ = T/N. We present numerical results only for the fully discrete scheme by the Galerkin FEM
in space and the L1 scheme in time, since  BE-CQ gives nearly identical results.

We consider the following two examples to illustrate the analysis.


  (a) f≡0 and u_d(x,t)=e^t x(1-x).

  (b) f=(1+cos(t))χ_(1/2,1)(x) and u_d(x,t)=5e^tx(1-x).

Throughout, unless otherwise specified, the penalty parameter γ is set to γ=1, and the lower and upper bounds a and b
in the admissible set  to a=0 and b=0.05. The final time T is fixed at T=0.1. The conditions from Theorems
<ref>, <ref> and <ref> are satisfied for both examples, and thus the error estimates therein hold.

In Tables <ref> and <ref>, we present the spatial error e_h(u) in the L^∞
(0,T;L^2(Ω))-norm for the semidiscrete solution u_h, defined by

    e_h(u) = max_1≤ n≤ Nu_h(t_n)-u(t_n)_L^2(Ω),

and similarly for the approximations z_h and q_h. The numbers in the bracket denote the theoretical rates. Since
the exact solution to problem (<ref>) is unavailable, we
compute reference solutions on a finer mesh, i.e., the continuous solution u(t_n) with a fixed time step τ=T/1000 and
mesh size h=1/1280. The empirical rate for the spatial error e_h is of order O(h^2), which is consistent with the theoretical result in
Theorem <ref>. For case (a), the box constraint is inactive, and thus the errors for the
control q and adjoint z are identical (since γ=1).





Next, to examine the convergence in time, we compute the ℓ^2(L^2(Ω)) and ℓ^∞
(L^2(Ω)) temporal errors e_τ,2(u) and e_τ,∞(u) for the fully discrete solutions U_h^n,
respectively, defined by

    e_τ,2(u) =   (U_h^n-u_h(t_n))_n=1^N_ℓ^2(L^2(Ω))    
      e_τ,∞(u) = max_1≤ n≤ NU_h^n-u_h(t_n)_L^2(Ω),

and similarly for the approximations Z_h^n and Q_h^n.
The reference semidiscrete solution u_h is computed with h=1/50 and τ=1/(64×10^4).
Numerical experiments show that the empirical rate for the temporal discretization error is of
order O(τ^min(1/2+α,1)) and
O(τ^α) in the ℓ^2(L^2(Ω)) and ℓ^∞(L^2(Ω))-norms, respectively,
cf. Tables <ref>–<ref> and <ref>–<ref>, for
cases (a) and (b). These results agree well with the theoretical predictions from Theorems
 <ref> and <ref>, and thus fully support the
error analysis in Section <ref>. In Fig. <ref>,
we plot the optimal control q, the state u and the adjoint z. One clearly observes the weak
solution singularity at t=0 for the state u and at t=T for the adjoint z. The latter is
especially pronounced for case (b). The weak solution singularity is due to the incompatibility of the
source term with the zero initial/terminal data.












§ CONCLUSIONS


In this work, we have developed a complete numerical analysis of a fully discrete scheme for a distributed
optimal control problem governed by a subdiffusion equation, with box constraint on the control
variable, and derived nearly sharp pointwise-in-time error estimates for both space and time discretizations. These estimates agree well with
the empirical rates observed in the numerical experiments. The theoretical and numerical results show  the adverse
influence of the fractional derivatives on the convergence rate when the fractional order α is small.














§ PROOF OF LEMMA <REF>


By Sobolev embedding, W^s,p(0,1;L^2(Ω))↪ C([0,1];L^2(Ω)) for s∈ (1/p,1], and thus
we can define an interpolation operator Π by Π v(t̂) = v(1), for t̂∈(0,1), for any v∈ W^s,p(0,1;L^2(Ω)).
The operator E =I-Π is bounded from W^s,p(0,1;L^2(Ω)) to L^2(0,1;L^2(Ω)):

    Ev_L^2(0,1;L^2(Ω)) =  (I-Π)v_L^2(0,1;L^2(Ω))≤ c  v _W^s,p(0,1;L^2(Ω)).

By the fractional Poincaré inequality (cf. <cit.>), we have

    E v_L^p(0,1;L^2(Ω))   =   inf_p∈ℝ E  (v - p) _L^p(0,1;L^2(Ω))≤   c inf_p∈ℝ v- p  _W^s,p(0,1;L^2(Ω))
       ≤  c | v|_W^s,p(0,1;L^2(Ω)) ,

where the seminorm |·|_W^s,p(0,T;L^2(Ω)) is defined in (<ref>).
By Hölder's inequality, we obtain

    (v(t_n) - v̅^n)_n=1^N  _ℓ^p(L^2)  ^p
         = τ∑_n=1^N v(t_n)-τ^-1∫_t_n-1^t_nv(t)ṭ_L^2(Ω)^p
    
       =   τ^1-p∑_n=1^N ∫_t_n-1^t_n(v(t_n)-v(t))ṭ_L^2(Ω)^p
        ≤∑_n=1^N ∫_t_n-1^t_nv(t_n) - v(t)_L^2(Ω)^p ṭ.

Let v_n(t̂)=v(t_n-1+τt̂), for t̂∈[0,1], n=1,…,N.
Then v_n(t̂) ∈ W^s,p(0,1;L^2(Ω)) and by (<ref>), we have

    (v(t_n) - v̅^n)_n=1^N  _ℓ^p(L^2)  ^p
           ≤τ∑_n=1^N ∫_0^1Πv_n  -  v_n_L^2(Ω)^p( t̂) ̣̂t
        ≤ c τ∑_n=1^N  | v_n|_W^s,p(0,1;L^2(Ω))^p 
       ≤ c τ∑_n=1^N   ∫_0^1 ∫_0^1 v_n(t̂)-v_n(ξ̂)_L^2(Ω)^p/|t̂-ξ̂|^1+ps ̣̂t ̣̂ξ
       = c τ^ps∑_n=1^N   ∫_t_n-1^t_n∫_t_n-1^t_nv(t)-v(ξ)_L^2(Ω)^p/|t-ξ|^1+ps ṭξ̣
    
           ≤ c τ^ps∫_0^T  ∫_0^T v(t)-v(ξ)_L^2(Ω)^p/|t-ξ|^1+ps ṭξ̣=  c τ^ps  |v|_W^s,p(0,T;L^2(Ω))^p ,

which implies the desired assertion.


abbrv


