Stability and numerical analysis of backward problem for subdiffusion with time-dependent coefficients

Let $\Omega\subset\mathbb{R}^d$ ($d = 1,2,3$) be a convex polyhedral domain with boundary $\partial\Omega$. We are interested in the fractional evolution model with time-dependent coefficient:

$$\begin{align} \begin{aligned} \partial_t^\alpha u (x,t)+ \nabla \cdot(a(x,t)\nabla u) &= f(x,t),\ &&\text{in}\ \Omega\times(0,T],\\ u(x,t) &= 0,\ &&\text{on}\ \partial\Omega, \\ u(x,0)& = u_0(x),\ &&\text{in}\ \Omega, \\ \end{aligned} \end{align} \tag{ 1.1 }$$

where T > 0 is a fixed final time, $f \in L^\infty(0,T;L^2(\Omega))$ and $u_0\in L^2(\Omega)$ are given source term and initial data, respectively. $a(x,t)\in\mathbb{R}^{d\times d}$ is a symmetric matrix-valued diffusion coefficient such that for constants $c_0\geqslant 1$ and $c_1\gt0$

$$\begin{align} c_0^{-1}|\xi|^2\leqslant a(x,t)\xi\cdot\xi \leqslant c_0 |\xi|^2,\quad\quad\quad\quad~~ \forall\, \xi\in\mathbb{R}^d, \,\, \forall\, (x,t)\in \Omega\times\mathbb{R}^+, \end{align} \tag{ 1.2 }$$

$$\begin{align} |\partial_t a(x,t)|+|\nabla_xa(x,t)|+|\nabla_x\partial_t a(x,t)|\leqslant c_1,\quad \forall\, (x,t)\in \Omega\times\mathbb{R}^+. \end{align} \tag{ 1.3 }$$

Here · and $|\cdot|$ denote the standard Euclidean inner product and norm, respectively, and $\mathbb{R}^+ = [0,\infty)$. In (1.1), ${\partial_t^\alpha} u(t)$ denotes the Caputo fractional derivative in time t of order $\alpha\in(0,1)$ [16, p 70]

$$\begin{align*} {\partial_t^\alpha} u(t)= \frac{1}{\Gamma(1-\alpha)}\int_0^t(t-s)^{-\alpha} \tfrac{\partial}{\partial s} u(s){\,\mathrm{d}}s, \end{align*}$$

where $\Gamma(z) = \int_0^\infty s^{z-1}e^{-s}{\,\mathrm{d}}s$(for $\Re(z)\gt0$) denotes the Euler's Gamma function. As the order $\alpha\to 1^-$, the fractional derivative converges to the standard first derivative u^ʹ for sufficiently smooth u. Due to its numerous applications in physics, engineering, biology, and finance, there has been a surge in interest in fractional/nonlocal models in recent years. In particular, the time-fractional diffusion equation models the mean squared displacement of particles grows only sublinearly with time, compared to classical diffusion process growing linearly in time, so it often is called subdiffusion model. The model is extensively used to describe subdiffusion process in nature, such as highly heterogeneous aquifers in media [1, 5] and fractal geometry [31]. For interested readers, see a long list of applications of fractional models discovered from biology and physics in [28, 29].

Inverse problems for subdiffusion models have been extensively studied and there has already been a vast literature; we recommend review papers [14, 20, 21] and reference therein. In this paper, we focus on backward problem for the subdiffusion model (1.1): to recover the initial data $u_0(x)$ with $x\in\Omega$ from terminal observation

$$\begin{equation*} u(x,T) = g(x),~~\text{for all }x\in \Omega. \end{equation*}$$

In recent years, the backward subdiffusion problem has received a lot of attention. All existing works are for the case that the coefficients are independent of time, a.e. $a(x,t) \equiv a(x)$. The uniqueness and some stability estimate can be found in the pioneer work [32]. We also refer [22, 35, 37–39] for different regularization methods, and [41] for error analysis of fully discrete schemes. The analysis heavily relies on the asymptotic behaviors of Mittag–Leffler functions, or equivalently the smoothing properties of solution operators. Unfortunately, this strategy is not directly applicable for subdiffusion models with time dependent coefficients. Note that the backward problem for the parabolic equation with time-dependent coefficient has been intensively studied. However, the studies of fractional diffusion model (1.1) is much more challenging since many useful mathematical tools, including product rule and chain rule, are not directly applicable.

There have been some existing works for the direct problem of time-fractional model (1.1). For time-dependent elliptic operators or nonlinear problems, energy arguments [36] or perturbation arguments [17] can be used to show existence and uniqueness of the solution. However, more refined stability estimates, needed for numerical analysis of nonsmooth problem data, often have to be derived separately. Mustapha [30] analyzed the spatially semidiscrete Galerkin finite element method (FEM) approximation of problem (1.1) using a novel energy argument, and established optimal-order convergence rates for both smooth and nonsmooth initial data. See also [25–27] for time-fractional advection diffusion equation. In [12], a perturbation argument of freezing the diffusion coefficients was proposed to analyze the partial differential equation (PDE) (1.1) and its numerical treatment. The argument was then modified and adapted to the error analysis of high-order discretization scheme in [13]. However, the analysis for the uniqueness and stability of backward problem is still missing in the literature. We also refer interested readers to [4, 6] for the inverse source problem with time-dependent coefficients, where the uniqueness was proved using some nonstandard energy argument.

The first contribution of the paper is to develop a conditional stability in Sobolev spaces. Under some assumptions on diffusion coefficients and some smallness/largeness conditions on terminal time T, there holds the Lipschitz stability (theorems 2.2 and 2.4)

$$\begin{equation} \| u_0 \|_{L^2(\Omega)} \leqslant c(1+T^{\,\alpha}) \| u(T) \|_{H^2(\Omega)}. \end{equation} \tag{ 1.4 }$$

The proof heavily relies on several a priori estimates of the direct problem, the smoothing properties of solution operators with frozen diffusion coefficients, and a perturbation argument. The stability estimate (1.4) plays a key role in the analysis of regularization scheme and completely discrete approximation.

In practice, the observational data often involves random noise. In this work, we denote the empirical observation by g_δ and assume that it is noisy with a level δ > 0 in the sense that

$$\begin{equation} \|g_\delta-g\|_{L^2(\Omega)} = \delta. \end{equation} \tag{ 1.5 }$$

According to the stability estimate (1.4), we know the backward problem is mildly ill-posed, and it is necessary to apply some regularization in case of noisy data. In this work, we propose a quasi-boundary value method for regularization and analyze convergence of the regularized solution. Specifically, under the assumptions in the stability estimate (1.4), if $\| u_0 \|_{\dot H^{q}(\Omega)}\leqslant c$ with $q\in (0,2]$, we prove that (theorem 3.1)

$$\begin{equation*} \| {u_\gamma^\delta}(0) - u_0 \|_{L^2(\Omega)} \leqslant c \,\Big(\delta \gamma^{-1} + \gamma^{\frac q2}\Big). \end{equation*}$$

Here γ denotes the regularization parameter and ${u_\gamma^\delta}(0)$ denotes the reconstruction via the regularization scheme. Then the optimal convergence rate (in noise level) is $O(\delta^{\frac{q}{q+2}})$ provided that $\gamma \sim \delta^{\frac2{q+2}}$. Moreover, for $u_0\in L^2(\Omega)$, there holds

$$\begin{equation*} \| {u_\gamma^\delta}(0) - u_0\|_{L^2(\Omega)} \rightarrow 0\quad \text{as}~~\delta,~\gamma\rightarrow0 ~~\text{and}~~\frac\delta\gamma\rightarrow0. \end{equation*}$$

The next contribution of this paper is to develop a fully discrete scheme with thorough error analysis. To numerically recover the initial data, we discretize the proposed regularization scheme by using piecewise linear FEM in space with mesh size h, and backward Euler convolution quadrature scheme (CQ-BE) in time with time step size τ. Then the numerical discretization introduces additional discretization error. We establish an a priori error bound for the fully discrete scheme. This estimate provides a useful guideline to choose suitable discretization parameters h and τ and regularization parameter γ according to the a priori known noise level δ. In particular, let $U_{h,\gamma}^{0,\delta}$ be the numerical reconstruction of initial condition. Under those assumptions in the stability estimate (1.4), we show that (theorem 4.1) there holds

$$\begin{equation*} \|U_{h,\gamma}^{0,\delta}-u_0\|_{L^2(\Omega)}\leqslant c\Big(\gamma^\frac q2+ \delta\gamma^{-1} + h^2\gamma^{-1}+\tau |\log \tau| (h^2\gamma^{-1}+1) \Big), \end{equation*}$$

if $\| u_0 \|_{\dot H^{q}(\Omega)}\leqslant c$ with $q\in (0,2]$. Then with the choice $\gamma \sim \delta^{\frac{2}{q+2}}$, $h\sim \delta^\frac 12$ and $\tau |\log \tau| \sim \delta^{\frac{q}{q+2}}$, we obtain the optimal approximation error of order $O(\delta^{\frac{q}{q+2}})$. Moreover, for $u_0\in L^2(\Omega)$, there holds

$$\begin{equation*} \| U_{h,\gamma}^{0,\delta} - u_0\|_{L^2(\Omega)} \rightarrow 0\quad \text{as}~~\delta,\gamma,h,\tau\rightarrow0, ~~\frac\delta\gamma\rightarrow0~~\text{and}~~\frac{h^2}{\gamma} \rightarrow0. \end{equation*}$$

In this proof, we combine several useful tools, including the error analysis of the direct problem with respect to the problem data [12, 13], the smoothing properties of discrete solution operators with frozen diffusion coefficients [41], and the analysis of the conditional stability estimate (1.4) in the discrete sense.

The rest of the paper is organized as follows. In section 2 we provide some preliminary results about solution regularity, smoothing properties of solution operators and derive conditional stability of the inverse problem. In section 3 we discuss the regularization scheme by quasi-boundary value method. In section 4 we propose and analyze a fully discrete scheme for solving the backward problem. Finally, in section 5 we present some numerical examples to illustrate and complete the theoretical analysis.

Here we introduce some notations used throughout the paper. Under conditions (1.2) and (1.3), the abstract time-dependent elliptic operator $A(t): H_0^1(\Omega)\cap H^2(\Omega)\rightarrow L^2(\Omega)$ is defined by

$$\begin{align*} A(t)\phi=- \nabla\cdot(a(x,t)\nabla \phi), \end{align*}$$

with the domain $\text{Dom}(A(t)) = H^1_0(\Omega)\cap H^2(\Omega)$ for all $t\in[0,T]$. By the complex interpolation method [34], this implies

$$\begin{align*} \text{Dom}(A(t)^\gamma)=\dot H^{2\gamma}(\Omega)=(L^2(\Omega),H^1_0(\Omega)\cap H^2(\Omega))_{[\gamma]} ,\quad\forall\, t\in[0,T], \,\,\,\forall\,\gamma\in[0,1]. \end{align*}$$

Equivalently, it can be defined via spectral decomposition of the operator A(t) [33, chapter 3]. Let $\{(\lambda_j,\varphi_j)\}_{j = 1}^n$ be the eigenpairs of $A(t_*)$ for a fixed $t_* \in [0,T]$ with multiplicity counted and $\{\varphi_j\}_{j = 1}^\infty$ be an orthonormal basis in $L^2(\Omega)$. Then the Hilbert space $\dot{H}^{\gamma} (\Omega)$ can be equivalently defined as

$$\begin{equation*} \dot H^{\gamma}(\Omega) = \Big\{v\in L^2(\Omega): \sum_{j=1}^\infty \lambda_j^{\gamma}(v,\varphi_j)^2\lt\infty\Big\}. \end{equation*}$$

For $\gamma\in[0,2]$ we also denote by $\dot H^{-\gamma}(\Omega)$ the dual space of $\dot H^{\gamma}(\Omega)$. Then the norm of $\dot H^{-\gamma}(\Omega)$ satisfies

$$\begin{equation*} \|v\|_{\dot H^{-\gamma}(\Omega)} \sim \|A(t)^{-\frac{\gamma}{2}}v\|_{L^2(\Omega)}\quad\forall\, v\in\dot H^{-\gamma}(\Omega),\,\,\forall\, t\in[0,T]. \end{equation*}$$

First we recall basic properties of the subdiffusion model with a time-independent diffusion coefficient, i.e. $a(x,t_*)$ for some $t_*\geqslant0$. Accordingly, consider the problem

$$\begin{equation} {\partial_t^\alpha} u(t) + A(t_*)u(t) = f(t)\ \,\,\, \forall t\in(0,T], \quad \mbox{with }u(0)=u_0. \end{equation} \tag{ 2.1 }$$

By means of Laplace transform, the solution u(t) can be represented by [11, section 4]

$$\begin{align} u(t)= F(t;t_*)u_0 + \int_0^t E(t-s;t_*) f(s) {\,\mathrm{d}}s , \end{align} \tag{ 2.2 }$$

where the solution operators $F(t;t_*)$ and $E(t;t_*)$ are defined by

$$\begin{align} F(t;t_*) &= \frac{1}{2\pi\mathrm{i}} \int_{\Gamma_{\theta,\kappa}} e^{zt} z^{\alpha-1} (z^\alpha + A(t_*))^{-1} \,{\,\mathrm{d}}z,~~\text{and} \nonumber\\ E(t;t_*) &= \frac{1}{2\pi\mathrm{i}} \int_{\Gamma_{\theta,\kappa}} e^{zt} (z^\alpha +A(t_*))^{-1}\,{\,\mathrm{d}}z, \end{align} \tag{ 2.3 }$$

with integration over a contour $\Gamma_{\theta,\kappa}\subset\mathbb{C}$ (oriented with an increasing imaginary part):

$$\begin{equation*} \Gamma_{\theta,\kappa}=\left\{z\in \mathbb{C}: |z|=\kappa, |\arg z|\leqslant \theta\right\}\cup \{z\in \mathbb{C}: z=\rho e^{\pm\mathrm{i}\theta}, \rho\geqslant \kappa\} . \end{equation*}$$

Throughout, we fix $\theta \in(\frac{\pi}{2},\pi)$ so that $z^{\alpha} \in \Sigma_{\alpha\theta}\subset \Sigma_{\theta}: = \{0\neq z\in\mathbb{C}: \textrm{arg}(z)\leqslant\theta\},$ for all $z\in\Sigma_{\theta}$.

The next lemma gives smoothing properties and asymptotics of $F(t;t_*)$ and $E(t;t_*)$. The proof follows from the resolvent estimate [2, example 3.7.5 and theorem 3.7.11]:

$$\begin{equation} \| (z+A)^{-1} \|\leqslant c_\phi (|z|^{-1},\lambda^{-1}) \quad \forall z \in \Sigma_{\phi}, \,\,\,\forall\,\phi\in(0,\pi) , \end{equation} \tag{ 2.4 }$$

where $\|\cdot\|$ denotes the operator norm from $L^2(\Omega)$ to $L^2(\Omega)$, and λ denotes the smallest eigenvalue of $-\Delta$ with homogeneous Dirichlet boundary condition. The proof of (a) and (b) were given in [8, theorems 6.4 and 3.2], and (c) were proved by Sakamoto and Yamamoto in [32, theorem 4.1].

Lemma 2.1. Let $F(t;t_*)$ and $E(t;t_*)$ be the solution operators defined in (2.3) for any $t_*\geqslant0$ Then they satisfy the following properties for all t > 0

• (a)  
  $\|A(t_*) F (t;t_*)v\|_{L^2(\Omega)} + t^{1-(2-k)\alpha} \| A(t_*)^{k} E (t;t_*)v \|_{L^2(\Omega)} \leqslant c t^{-\alpha} \|v\|_{L^2(\Omega)}$ with $k = 1,2$;
• (b)  
  $\|F(t;t_*)v\|_{L^2(\Omega)}+ t^{1-\alpha}\|E(t;t_*)v\|_{L^2(\Omega)} \leqslant c \min(1, t^{-\alpha}) \|v\|_{L^2(\Omega)}$;
• (c)  
  $\| F(t; t_*)^{-1} v \|_{L^2(\Omega)} \leqslant c (1+ t^\alpha) \| v \|_{\dot H^{2}(\Omega)}$ for all $v\in \dot H^{2}(\Omega)$.

The constants in all above estimates are uniform in t, but they are only dependent of $t_*$ and T.

Next, we turn to the subdiffusion with a time-dependent coefficient. The overall proof strategy is to employ a perturbation argument [12], and then to properly resolve the singularity. Specifically, for any fixed $t_*\in(0,T]$, we rewrite problem (1.1) into

$$\begin{align} \left\{\begin{aligned}{\partial_t^\alpha} u(t) + A(t_*)u(t) &= (A(t_*)-A(t))u(t) + f(t) ,\quad \forall t\in(0,T],\\ u(0)&=u_0. \end{aligned} \right. \end{align} \tag{ 2.5 }$$

By (2.2), the solution u(t) of (2.5) is given by

$$\begin{align} u(t)&= F(t;t_*)u_0+ \int_0^t E(t-s;t_*)(f(s) + (A(t_*)-A(s))u(s)){\,\mathrm{d}}s. \end{align} \tag{ 2.6 }$$

The following perturbation estimate will be used extensively. See similar results in [12, corollary 3.1].

Lemma 2.2. Under conditions (1.2) and (1.3), there holds that

$$\begin{equation*} \|(A(t)-A(s))v\|_{\dot H^{p}(\Omega)}\leqslant c\min(1,|t-s|)\|v\|_{\dot H^{p+2}(\Omega)},~~p\in [-2,0]. \end{equation*}$$

Proof. The condition (1.3) implies the case that p = 0. The case $p = -2$ is has been proved in [12, corollary 3.1]. Then the intermediate case follows from the interpolation [7, section 2.5].

Next, we state a few regularity results. The proof of these results can be found in, e.g. [3, 12, 32].

Theorem 2.1. Let u(t) be the solution to (1.1). Then the following statements hold.

• (a)  
  If $u_0 \in \dot H^{q}(\Omega)$ with $s\in[0,2]$ and f = 0, then there holds

  $$\begin{equation*} \| \partial_t^{(m)} u(t) \|_{\dot H^{p}(\Omega)} \leqslant c t^{\frac{(s-p)\alpha}{2}-m} \| u_0 \|_{\dot H^{q}(\Omega)} \end{equation*}$$

  with $0\leqslant p-q\leqslant 2$ and $m = 0,1$. The constant c in the estimate depends on T and α.
• (b)  
  If $u_0 = 0$ and $f\in L^p(0,T;L^2(\Omega))$ with $1\lt p\lt\infty$, then there holds

  $$\begin{equation*} \| u\|_{L^p(0,T;\dot H^2(\Omega))} +\|{\partial_t^\alpha} u\|_{L^p(0,T;L^2(\Omega))} \leqslant c\|f\|_{L^p(0,T;L^2(\Omega))}. \end{equation*}$$

  Moreover, if $f\in L^p(0,T;L^2(\Omega))$ with $1/\alpha\lt p\lt\infty$, then u(t) is the solution to problem (1.1) such that $u\in C([0,T];L^2(\Omega))$. The constant c in the estimate depends on T and α.

The next lemma provides an a priori estimate similar to theorem 2.1 (a). Note that the generic constant in the new estimate is independent of T.

Lemma 2.3. Suppose that $u_0 \in L^2(\Omega)$ and f = 0. Let u(t) be the solution to the subdiffusion problem (1.1). Under conditions (1.2) and (1.3), there holds

$$\begin{align} \| u(t) \|_{L^2(\Omega)} &\leqslant c \min(1,t^{-\alpha}) \| u_0 \|_{L^2(\Omega)}~~\text{and} \nonumber\\ \| u(t) \|_{\dot H^{2}(\Omega)} &\leqslant c \,e^{ct} t^{-\alpha} \| u_0 \|_{L^2(\Omega)} \quad \text{for all}~~t\gt0. \end{align} \tag{ 2.7 }$$

Meanwhile, for any $\epsilon\in(0,1/\alpha-1)$ and t > 0, there holds that

$$\begin{equation} \| u(t)\|_{\dot H^{2}(\Omega)} \leqslant c t ^{-(1-\epsilon)\alpha}\|u_0\|_{L^2(\Omega)}. \end{equation} \tag{ 2.8 }$$

All the positive constants c in above estimates are independent of t and T.

Proof. We define an operator $\underline{A} = -c_0 \Delta$. Then by condition (1.2), the operator $A(t) - \underline{A}$ is selfadjoint and positive semidefinite for all $t\geqslant 0$. Then we rewrite the equation (1.1) as

$$\begin{align*} \partial_t^\alpha u(t) + \underline{A} u(t) = (\underline{A}-A(t)) u(t) \quad \text{for all}~~t\in(0,\infty). \end{align*}$$

Taking inner product with u(t) on the above equation and integrating by parts, we obtain

$$\begin{align*} (\partial_t^\alpha u(t), u(t)) + c_0 \| \nabla u(t) \|_{L^2(\Omega)}^2 = \big((c_0-a(\cdot,t)) \nabla u(t), \nabla u(t)\big)\leqslant 0 \quad \text{for all}~~t\in(0,\infty). \end{align*}$$

Using the facts that $(\partial_t^\alpha u(t), u(t)) \geqslant \| u(t) \|_{L^2(\Omega)} \partial_t^\alpha \| u(t) \|_{L^2(\Omega)}$ [8, lemma 6.1(iii)] and Poincaré inequality we arrive at

$$\begin{align*} \partial_t^\alpha \| u(t)\|_{L^2(\Omega)} + c \| u(t) \|_{L^2(\Omega)} \leqslant 0 \quad \text{for all}~~t\in(0,\infty), \end{align*}$$

for some constant c uniform in t. Then the comparison principle for fractional ordinary differential equation (ODEs) [18, theorem 2.3] leads to

$$\begin{align*} \| u(t)\|_{L^2(\Omega)} \leqslant E_{\alpha,1}(-c t^\alpha) \| u_0 \| \leqslant \frac{c}{1+ct^\alpha} \| u_0 \|_{L^2(\Omega)}. \end{align*}$$

This immediately leads to the desired claim (2.7).

Next, we apply the relation (2.6), lemmas 2.1 and 2.2 (with p = 2) to obtain for any $t_*\in(0,T]$

$$\begin{align*} \| u(t_*)\|_{\dot H^{2}(\Omega)} &\leqslant \| F(t_*;t_*) u_0\|_{\dot H^{2}(\Omega)} + c \int_0^{t_*} \| A(t_*) E(t_*-s;t_*) \| \, \| (A(t_*) -A(s))u(s) \|_{L^2(\Omega)}\,{\,\mathrm{d}}s \\ &\leqslant c t_*^{-\alpha} \| u_0 \|_{L^2(\Omega)} + c \int_0^{t_*} \| u(s) \|_{\dot H^{2}(\Omega)}\,{\,\mathrm{d}}s. \end{align*}$$

Then the Gronwall's inequality implies for any t > 0

$$\begin{align*}\| u(t) \|_{\dot H^{2}(\Omega)} \leqslant c \,e^{ct} t^{-\alpha} \| u_0 \|_{L^2(\Omega)}.\end{align*}$$

Meanwhile, lemma 2.2 leads to the estimate for $\beta = (1+\epsilon)\alpha$ with $\epsilon\in(0,1/\alpha-1)$

$$\begin{align*} \| u(t_*)\|_{\dot H^{2}(\Omega)} &\leqslant \| F(t_*;t_*) u_0\|_{\dot H^{2}(\Omega)} + \int_0^{t_*} \| A(t_*)^2 E(t_*-s;t_*) \| \, \| I - A(t_*)^{-1} A(s) \| \, \| u(s) \|_{L^2(\Omega)}\,{\,\mathrm{d}}s \\ & \leqslant c t_*^{-\alpha} \| u_0 \|_{L^2(\Omega)} + \int_{0}^{t_*}(t_*-s)^{-1+\epsilon\alpha} s^{-\alpha} \,{\,\mathrm{d}}s \leqslant c_\epsilon t_*^{-(1-\epsilon)\alpha}, \end{align*}$$

for any $t_*\gt0$. This completes the proof of (2.8).

Using the superposition principle, we consider the homogeneous source condition, i.e. $f\equiv 0$, without loss of generality. Then the corresponding backward subdiffusion problem reads: find u(0) such that

$$\begin{align} \partial_t^\alpha u + A(t) u = 0 \,\,\,\forall\, t\in(0,T] \quad \text{with} \quad u(T)&= g \qquad\mbox{in}\,\,\,\Omega. \end{align} \tag{ 2.9 }$$

The next theorem provides a stability estimate for the backward problem of (2.9) when T is sufficiently small.

Theorem 2.2. Suppose that $u_0 \in L^2(\Omega)$ and f = 0. Let u(t) be the solution to (1.1). Under conditions (1.2) and (1.3), there exists a positive constant T₀ such that for any $T\leqslant T_0$ there holds

$$\begin{align*} \| u_0 \|_{L^2(\Omega)} \leqslant c(1+T^{\,\alpha}) \| u(T) \|_{H^2(\Omega)}, \end{align*}$$

where the constant c depends on T₀ and T.

Proof. We rearrange the terms in relation (2.6) with $t_* = T$ to obtain

$$\begin{equation} u_0 = F(T;T)^{-1} \Big[u(T) - \int_0^T E(T-s;T) (A(T) - A(s)) u(s) {\,\mathrm{d}}s\Big]. \end{equation} \tag{ 2.10 }$$

Taking $L^2(\Omega)$ norm on both sides of the above relation, we apply lemma 2.1 (c) to obtain

$$\begin{align*}\begin{array}{*{20}{c}} {\parallel {u_0}{\parallel _{{L^2}({\textrm{\Omega }})}}}&{C(1 + {T^\alpha })\left( {\parallel u(T){\parallel _{{{\dot H}^2}({\textrm{\Omega }})}} + \mathop \smallint \limits_0^T \parallel A(T)E(T - s;T)\parallel } \right.}\\ {}&{\left. {\parallel (A(T) - A(s))u(s){\parallel _{{L^2}({\textrm{\Omega }})}}{\textrm{d}}s} \right).} \end{array}\end{align*}$$

According to lemma 2.2 with p = 0 and 2.1 (a) we arrive at

$$\begin{align*} \| u_0 \|_{L^2(\Omega)} & \leqslant c(1+T^{\,\alpha}) \Big(\| u(T) \|_{\dot H^{2}(\Omega)} + \int_0^{T} \| u(s) \|_{\dot H^{2}(\Omega)}\,{\,\mathrm{d}}s \Big). \end{align*}$$

Then this together with the estimate (2.7) implies

$$\begin{align*} \| u_0 \|_{L^2(\Omega)} &\leqslant c (1+T^{\,\alpha}) \Big(\| u(T) \|_{\dot H^{2}(\Omega)} + \int_0^{T} e^{cs} s^{-\alpha} \,{\,\mathrm{d}}s \| u_0 \|_{L^2(\Omega)}\Big)\\ &\leqslant c (1+T^{\,\alpha}) \Big(\| u(T) \|_{\dot H^{2}(\Omega)} + ce^{cT}T^{1-\alpha} \| u_0 \|_{L^2(\Omega)}\Big). \end{align*}$$

Let be the constant that

$$\begin{equation} c (1+T_0^\alpha)e^{cT_0}T_0^{1-\alpha} \lt \frac12,~~T_0\lt1. \end{equation} \tag{ 2.11 }$$

Then for any $T\leqslant T_0$

$$\begin{align*} \| u_0 \|_{L^2(\Omega)} \leqslant c (1+T^{\,\alpha}) \| u(T) \|_{\dot H^{2}(\Omega)} . \end{align*}$$

This completes the proof of the lemma.

Next, we derive a stability estimate for a large T. To this end, we need the following assumption.

Assumption 2.3. There exists constants $c_2\gt0$ and κ > 0 such that

$$\begin{align*} | \partial_t \nabla_x a(x,t)| + | \partial_t a(x,t)|\leqslant c_2 t^{-\kappa}\quad \forall\, (x,t)\in \Omega\times(0,\infty). \end{align*}$$

Under the condition, we have the following perturbation estimate. The proof is similar to that of lemma 2.2. The proof is provided in appendix A for completeness.

Lemma 2.4. Under conditions (1.2), (1.3) and assumption 2.3, there holds for all $t,s\geqslant1$

$$\begin{equation*} \|(A(t)-A(s))v\|_{\dot H^{p}(\Omega)} \leqslant c \min\big(1,\min(t,s)^{-\kappa}|t-s|\big)\|v\|_{\dot H^{p+2}(\Omega)},\qquad \forall~~p\in[-2,0]. \end{equation*}$$

The next theorem provides a stability result in case of sufficiently large T.

Theorem 2.4. Suppose that $u_0 \in L^2(\Omega)$ and f = 0. Let conditions (1.2), (1.3) and assumption 2.3 be valid. Let u(t) be the solution to the subdiffusion problem (1.1). Then there exists positive $T_1\gt1$ such that for any $T\geqslant T_1$ there holds

$$\begin{align*} \| u_0 \|_{L^2(\Omega)} \leqslant c (1+T^{\,\alpha}) \| u(T) \|_{H^2(\Omega)}, \end{align*}$$

where the constant c depends on T₁ and T.

Proof. Using (2.10) and taking L² norm on both sides, we apply again lemma 2.1 (c) to obtain

$$\begin{align*} \| u_0 \|_{L^2(\Omega)} & \leqslant c (1+T^{\,\alpha}) \Big(\| u(T) \|_{\dot H^{2}(\Omega)} + \Big\|\int_0^{T} A(T) E(T-s;T)\big( A(T) -A(s) \big) u(s) \,{\,\mathrm{d}}s \Big\|_{L^2(\Omega)}\Big). \end{align*}$$

Applying lemma 2.4 with $p = -2$, we have for sufficiently small ε > 0

$$\begin{align*}\|(I-A(t)^{-1}A(s))v\|_{L^2(\Omega)} \leqslant c \min(t,s)^{-\kappa\alpha(1+\epsilon)}|t-s|^{\alpha(1+\epsilon)}\|v\|_{L^2(\Omega)}. \end{align*}$$

This together with lemma 2.1 (a) and the a priori estimate (2.8), imply

$$\begin{align*} \| A(T) E(T-s;T)\big( A(T) \!-\!A(s) \big) u(s) \|_{L^2(\Omega)} &\!\leqslant\! \| A(T)^2 E(T-s;T)\| \, \| (I-A(T)^{-1}A(s)) u(s) \|_{L^2(\Omega)}\\ &\leqslant c T^{-\kappa\alpha(1+\epsilon)} (T-s)^{-1+\epsilon\alpha} \, \| u(s) \|_{L^2(\Omega)} \\ &\leqslant c T^{-\kappa\alpha(1+\epsilon)} (T-s)^{-1+\epsilon\alpha} \, s^{-(1-\epsilon)\alpha}\| u(0) \|_{L^2(\Omega)} \end{align*}$$

for all $s\in[T/2,T]$. Then we arrive at

$$\begin{align*} \Big\|\int_{{T}/2}^{T} A(T) E(T-s;T)\big( A(T) -A(s) \big) u(s) \,{\,\mathrm{d}}s \Big\|_{L^2(\Omega)} \leqslant c T^{-\alpha-\kappa\alpha(1+\epsilon)+2\epsilon\alpha} \| u_0 \|_{L^2(\Omega)} . \end{align*}$$

Meanwhile, we apply lemmas 2.1 and 2.2 again to derive

$$\begin{align*}\| A(T)^2 E(T-s;T)\| \, \| I - A(T)^{-1}A(s) \| \leqslant c (T-s)^{-1-\alpha}\quad \text{for all}~~s\in(0,T/2].\end{align*}$$

This together with the estimate (2.7) leads to

$$\begin{align*} &\int_0^{T/2} \| A(T)^2 E(T-s;T)\| \, \| I - A(T)^{-1}A(s) \| \, \| u(s) \|_{L^2(\Omega)}\,{\,\mathrm{d}}s \\ \leqslant c &\int_0^{T/2} (T-s)^{-1-\alpha} s^{-(1-\epsilon)\alpha} \,{\,\mathrm{d}}s \| u(0) \|_{L^2(\Omega)} \leqslant c T^{-(2-\epsilon)\alpha}\| u(0) \|_{\dot H^{2}(\Omega)}. \end{align*}$$

To sum up, we arrive at the estimate

$$\begin{align*} \| u_0 \|_{L^2(\Omega)} & \leqslant c (1+T^{\,\alpha}) \| u(T) \|_{\dot H^{2}(\Omega)} + c (1+T^{\,\alpha}) (T^{-\kappa\alpha(1+\epsilon)-\alpha+2\epsilon\alpha} + T^{-(2-\epsilon)\alpha})\| u_0 \|_{L^2(\Omega)}. \end{align*}$$

Then choosing a sufficiently small ε, there exists $T_1\gt1$ sufficiently large such that

$$\begin{equation} c (1+T_1^\alpha) (T_1^{-\kappa\alpha(1+\epsilon)-\alpha+2\epsilon\alpha} + T_1^{-(2-\epsilon)\alpha}) = \frac 12 \end{equation} \tag{ 2.12 }$$

and hence for any $T\geqslant T_1$, there holds the desired stability estimate.

In sections 3 and 4, we shall discuss respectively the regularization and a fully discrete scheme with rigorous numerical analysis. The stability estimate in theorems 2.2 and 2.4 provides a key tool in the coming numerical analysis. Therefore, from now on, we suppose the following assumption are valid.

Assumption 2.5. Suppose conditions (1.2) and (1.3) and one of the following conditions are valid.

• (a)  
  $T\leqslant T_0$, where T₀ be a sufficiently small constant;
• (b)  
  Assumption 2.3 holds and $T\geqslant T_1$ where T₁ be a sufficiently large constant.

In practice, the observational data often suffers from noise, i.e. (1.5). In this section, we study a simple regularization scheme by using the quasi boundary value method. Let ${u_\gamma^\delta}(t) \in \dot H^{1}(\Omega)$ be the regularizing solution such that

$$\begin{equation} \begin{aligned} \partial_t^\alpha {u_\gamma^\delta}(t) + A(t) {u_\gamma^\delta}(t) &= 0~~\,\,\,\forall\, t\in(0,T] \quad \text{with} \quad \gamma{u_\gamma^\delta}(0) + {u_\gamma^\delta}(T)&= g^\delta. \end{aligned} \end{equation} \tag{ 3.1 }$$

where γ denotes a positive regularization parameter. To derive an error estimate for ${u_\gamma^\delta}(0) - u(0)$, we introduce an auxiliary function ${u_\gamma}(t) \in \dot H^{1}(\Omega)$ satisfying

$$\begin{align} \partial_t^\alpha {u_\gamma}(t) + A(t) {u_\gamma}(t) = 0 \,\,\,\forall\, t\in(0,T] \quad \text{with} \quad \gamma{u_\gamma}(0) + {u_\gamma}(T)&= g. \end{align} \tag{ 3.2 }$$

Then using the solution representation

$$\begin{align*} {u_\gamma}(T) = F(T;T){u_\gamma}(0) + \int_0^T E(T-s;T) (A(T) - A(s)) {u_\gamma}(s) \,{\,\mathrm{d}}s \end{align*}$$

we have the relation

$$\begin{equation*} \gamma {u_\gamma}(0) + F(T;T){u_\gamma}(0) + \int_0^T E(T-s;T) (A(T) - A(s)) {u_\gamma}(s) \,{\,\mathrm{d}}s = g. \end{equation*}$$

Therefore, we derive

$$\begin{align} \begin{aligned} {u_\gamma}(0) &=(\gamma I + F(T;T))^{-1} \Big[ g - \int_0^T E(T-s;T) (A(T) - A(s)) {u_\gamma}(s) \,{\,\mathrm{d}}s \Big]. \end{aligned} \end{align} \tag{ 3.3 }$$

Similarly, we have

$$\begin{align} \begin{aligned} {u_\gamma^\delta}(0) &=(\gamma I + F(T;T))^{-1} \Big[ g^\delta - \int_0^T E(T-s;T) (A(T) - A(s)) {u_\gamma^\delta}(s) \,{\,\mathrm{d}}s \Big]. \end{aligned} \end{align} \tag{ 3.4 }$$

We begin with the following lemma on solution operator with fixed-time operator A(T). These estimates have been proved in [42, lemma 3.3] by means of spectral decomposition.

Lemma 3.1. Let $0\leqslant p\leqslant q\leqslant 2+p$. Then there holds the estimates for any $\gamma\in (0,1]$

$$\begin{align*} & \|(\gamma I+F(T;T))^{-1}\|_{\dot H^{p}(\Omega)}\leqslant c (1+T^{\,\alpha})^{\frac{q-p}{2}} \gamma^{-(1+\frac{p-q}{2})}\|v\|_{\dot H^{q}(\Omega)} \quad\text{and}\\ & \quad \|F(T;T)(\gamma I+F(T;T))^{-1}\|\leqslant c. \end{align*}$$

All the constants are independent of p, q, T and γ.

Also we need the following regularity of the regularized solution.

Lemma 3.2. Let ${u_\gamma}(t)$ be the solution to (3.2). Suppose conditions (1.2) and (1.3) and one of the following conditions are valid.

• (a)  
  $T\leqslant T_0$, where T₀ be a sufficiently small constant;
• (b)  
  Assumption 2.3 holds and $T\geqslant T_1$ where T₁ be a sufficiently large constant.

Then there holds for any $p\in [0,2]$,

$$\begin{equation*} \|{u_\gamma}(0)\|_{\dot H^{p}(\Omega)} \leqslant c \gamma^{-\frac p2}\|u_0\|_{L^2(\Omega)} \end{equation*}$$

where the constant c depends on T₀ and T₁.

Proof. By means of the representation (3.3), theorem 2.1 and lemma 3.1,

$$\begin{align*} \|{u_\gamma}(0)\|_{L^2(\Omega)} &= \Big\|(\gamma I+F(T;T))^{-1}\Big(g - \int_0^T E(T-s;T)(A(T)-A(s)){u_\gamma}(s){\,\mathrm{d}}s\Big)\Big\|_{L^2(\Omega)}\\ &\leqslant c(1+T^{\,\alpha})\|g\|_{\dot H^{2}(\Omega)}+\|F(T;T)^{-1} \int_0^TE(T-s;T)(A(T)-A(s)){u_\gamma}(s){\,\mathrm{d}}s\|_{L^2(\Omega)}\\ &\leqslant c_T\|u_0\|_{L^2(\Omega)}+\|F(T;T)^{-1} \int_0^TE(T-s;T)(A(T)-A(s)){u_\gamma}(s){\,\mathrm{d}}s\|_{L^2(\Omega)}. \end{align*}$$

Then the desired result with p = 0 follows immediately from the proof of theorems 2.2 and 2.4.

Next, we turn to the case that p = 2. Similarly, we apply the representation (3.3) and lemma 3.1 again to obtain

$$\begin{align*} \| {u_\gamma}(0)\|_{\dot H^{2}(\Omega)} &\leqslant c \Big\|A(T)(\gamma I+F(T;T))^{-1}\Big(g - \int_0^T E(T-s;T)(A(T)-A(s)){u_\gamma}(s){\,\mathrm{d}}s\Big)\Big\|_{L^2(\Omega)}\\ &\leqslant c\gamma^{-1} \|g \|_{\dot H^{2}(\Omega)} +\int_0^T \|F(T;T)^{-1} A(T) E(T-s;T)(A(T)-A(s)){u_\gamma}(s)\|_{L^2(\Omega)} {\,\mathrm{d}}s \\ &\leqslant c_T \gamma^{-1} \| u_0 \|_{L^2(\Omega)} +c(1+T^{\,\alpha})\int_0^T \| A(T)^2 E(T-s;T)(A(T)-A(s)){u_\gamma}(s)\|_{L^2(\Omega)} {\,\mathrm{d}}s. \end{align*}$$

Using lemma 2.3 and Poincare inequality, we have

$$\begin{align} \| {u_\gamma}(t)\|_{\dot H^{2}(\Omega)} \leqslant ce^{ct} t^{-\alpha} \| {u_\gamma}(0) \|_{\dot H^{2}(\Omega)} \quad \text{and}\quad \| {u_\gamma}(t)\|_{\dot H^{2}(\Omega)} \leqslant c t ^{-(1-\epsilon)\alpha}\|{u_\gamma}(0)\|_{\dot H^{2}(\Omega)}, \end{align} \tag{ 3.5 }$$

with any small parameter ε > 0 and t > 0, and all the positive constants c in above estimates are independent of t and T. Next, we repeat the argument in theorems 2.2 and 2.4. Now lemmas 2.1 and 2.2 (with p = 0) imply that

$$\begin{align*} \|{u_\gamma}(0)\|_{\dot H^{2}(\Omega)}&\leqslant c\gamma^{-1} T^{-\alpha}\| u_0 \|_{L^2(\Omega)} +c(1+T^{\,\alpha})\int_0^T \| A(T)^2 E(T-s;T)(A(T)-A(s)){u_\gamma}(s)\|_{L^2(\Omega)} {\,\mathrm{d}}s \\ &\leqslant c\gamma^{-1} T^{-\alpha}\| u_0 \|_{L^2(\Omega)} + c(1+T^{\,\alpha}) \int_0^T (T-s)^{-\alpha}\| {u_\gamma}(s)\|_{\dot H^{2}(\Omega)} {\,\mathrm{d}}s\\ &\leqslant c\gamma^{-1} T^{-\alpha}\| u_0 \|_{L^2(\Omega)} + c(1+T^{\,\alpha}) \int_0^T (T-s)^{-\alpha} e^{cs}s^{-\alpha} \| {u_\gamma}(0)\|_{\dot H^{2}(\Omega)} {\,\mathrm{d}}s. \end{align*}$$

We combine this and (3.5) to arrive at

$$\begin{align*} \|{u_\gamma}(0)\|_{\dot H^{2}(\Omega)} \leqslant c\gamma^{-1} T^{-\alpha}\| u_0 \|_{L^2(\Omega)}+ c(1+T^{\,\alpha})T^{1-2\alpha}e^{cT}\|{u_\gamma}(0)\|_{\dot H^{2}(\Omega)}. \end{align*}$$

Then by choosing small T₀ such that $c(1+T_0^\alpha)T_0^{1-2\alpha}e^{cT_0}\lt\frac 12$, we arrive at

$$\begin{equation*} \|{u_\gamma}(0)\|_{\dot H^{2}(\Omega)}\leqslant c\gamma^{-1} T^{-\alpha }\|u_0\|_{L^2(\Omega)}\qquad \text{for all}~~T\in(0,T_0). \end{equation*}$$

Next we consider the case that T is sufficiently large, and we let assumption 2.3 be valid. Then we apply lemma 2.4 with p = 0 to arrive at

$$\begin{align*}\|(A(t) - A(s))v\|_{L^2(\Omega)} \leqslant c \min(t,s)^{-\kappa\alpha(1+\epsilon)}|t-s|^{\alpha(1+\epsilon)}\big)\|v\|_{\dot H^{2}(\Omega)} \end{align*}$$

for sufficiently small ε. This together with lemma 2.1 and the estimate (3.5) lead to

$$\begin{align*} \| A(T)^2 E(T-s;T)\big( A(T) -A(s) \big) {u_\gamma}(s)\|_{L^2(\Omega)} &\leqslant\| A(T)^2 E(T-s;T)\| \, \|\big( A(T) -A(s) \big) {u_\gamma}(s) \|_{L^2(\Omega)}\\ &\leqslant c T^{-\kappa\alpha(1+\epsilon)} (T-s)^{-1+\epsilon\alpha}\,\| {u_\gamma}(s)\|_{\dot H^{2}(\Omega)} \\ &\leqslant c T^{-\kappa\alpha(1+\epsilon)} (T-s)^{-1+\epsilon\alpha}\,s^{-(1-\epsilon)\alpha}\| {u_\gamma}(0) \|_{\dot H^{2}(\Omega)} \end{align*}$$

for all $s\in[T/2,T]$. Then we arrive at

$$\begin{align*} \int_{{T}/2}^{T} \big\| A(T)^2 E(T-s;T)\big( A(T) -A(s) \big) {u_\gamma}(s)\big\|_{L^2(\Omega)} \,{\,\mathrm{d}}s \leqslant c T^{-\alpha-\kappa\alpha(1+\epsilon)+2\epsilon\alpha} \| {u_\gamma}(0) \|_{\dot H^{2}(\Omega)} . \end{align*}$$

Meanwhile,we apply lemmas 2.1 and 2.2 again to derive

$$\begin{align*} \|A(T)^2 E(T-s;T)(A(t)-A(s))v\|_{L^2(\Omega)}\leqslant c(T-s)^{-1-\alpha}\|v\|_{\dot H^{2}(\Omega)}\quad\text{for all}~~s\in(0,T/2].\end{align*}$$

This together with the estimate (3.5) leads to

$$\begin{align*} &\quad \int_0^{T/2} \big\| A(T)^2 E(T-s;T)\big( A(T) -A(s) \big) {u_\gamma}(s)\big\|_{L^2(\Omega)} \,{\,\mathrm{d}}s \\ &\leqslant c \int_0^{T/2} (T-s)^{-1-\alpha} s^{-(1-\epsilon)\alpha} \,{\,\mathrm{d}}s \| {u_\gamma}(0) \|_{\dot H^{2}(\Omega)} \leqslant c T^{-(2-\epsilon)\alpha}\| {u_\gamma(0)} \|_{\dot H^{2}(\Omega)}. \end{align*}$$

To sum up, we arrive at the estimate

$$\begin{align*} \| {u_\gamma(0)} \|_{L^2(\Omega)} & \leqslant c\gamma^{-1} T^{-\alpha}\| u_0 \|_{L^2(\Omega)} + c (1+T^{\,\alpha}) (T^{-\kappa\alpha(1+\epsilon)-\alpha+2\epsilon\alpha} + T^{-(2-\epsilon)\alpha})\| u_0 \|_{L^2(\Omega)}. \end{align*}$$

Then choosing a sufficiently small ε, there exists $T_1\gt1$ sufficiently large such that

$$\begin{equation} c (1+T_1^\alpha) \left(T_1^{-\kappa\alpha(1+\epsilon)-\alpha+2\epsilon\alpha} + T_1^{-(2-\epsilon)\alpha}\right) = \frac 12 \end{equation} \tag{ 3.6 }$$

and hence for any $T\geqslant T_1$, there holds the desired stability estimate for p = 2.

The following lemma is about the estimate of the regularization with the backward solution.

Lemma 3.3. Let u and u_γ be the solutions to the backward problem (2.9) and regularized problem (3.2), respectively. Suppose assumption 2.5 is valid. Then if $u_0\in \dot H^{q}(\Omega)$ with $q\in (0,2]$ there holds

$$\begin{equation*} \| {u_\gamma}(0) - u(0) \|_{L^2(\Omega)} \leqslant c \gamma^{\frac q2} \| u_0 \|_{\dot H^{q}(\Omega)}, \end{equation*}$$

where the constant c depends on T₀ and T₁. Moreover, for $u_0\in L^2(\Omega)$, there holds

$$\begin{equation*} \lim_{\gamma\rightarrow0^+}\|{u}_\gamma(0) - u(0)\|_{L^2(\Omega)} =0 . \end{equation*}$$

Proof. We let $e: = {u_\gamma} - u$, it would satisfy

$$\begin{equation*} {\partial_t^\alpha} e+A(t) e = 0,~~\gamma e(0) + e(T) = -\gamma u_0, \end{equation*}$$

which further implies

$$\begin{equation} e(0) =(\gamma I + F(T;T))^{-1} \Big[-\gamma u(0) - \int_0^T E(T-s;T) (A(T) - A(s)) e(s) \,{\,\mathrm{d}}s \Big], \end{equation} \tag{ 3.7 }$$

lemma 3.1 implies its estimate that

$$\begin{align*} \| e(0)\|_{L^2(\Omega)}&\leqslant c\gamma^\frac{q}{2}\|u_0\|_{\dot H^{q}(\Omega)} + \|(\gamma I+F(T;T))^{-1}\int_0^T E(T-s;T)(A(T)-A(s)) e(s){\,\mathrm{d}}s\|_{L^2(\Omega)}\\ &\leqslant c\gamma^\frac{q}{2}\|u_0\|_{\dot H^{q}(\Omega)} + \|F(T;T)^{-1}\int_0^T E(T-s;T)(A(T)-A(s)) e(s){\,\mathrm{d}}s\|_{L^2(\Omega)}. \end{align*}$$

Then the desired result follows immediately from the proof of theorems 2.2 and 2.4.

Next, we consider the case that $u_0 \in L^2(\Omega)$. For an arbitrary $\tilde u_0 \in \dot H^{2}(\Omega)$, let $\tilde u(t)$ and $\tilde u_\gamma(t)$ be the functions respectively satisfying

$$\begin{align*} \partial_t^\alpha \tilde u(t) + A(t) \tilde u(t) = 0 \,\,\,\forall\, t\in(0,T] \quad \text{with} \quad \tilde u (0) = \tilde u_0, \end{align*}$$

and

$$\begin{align*} \partial_t^\alpha \tilde u_\gamma(t) + A(t) \tilde u_\gamma(t) = 0 \,\,\,\forall\, t\in(0,T] \quad \text{with} \quad \gamma\tilde u_\gamma(0) + \tilde u_\gamma(T)&= \tilde u(T). \end{align*}$$

We have proved that

$$\begin{equation*} \| \tilde u_\gamma(0) - \tilde u(0) \|_{L^2(\Omega)} \leqslant c \gamma \| \tilde u_0 \|_{\dot H^{2}(\Omega)}. \end{equation*}$$

Meanwhile, using the argument in theorems 2.2 and 2.4, we have

$$\begin{align*} \| \tilde u_\gamma(0) - u_\gamma(0) \|_{L^2(\Omega)} \leqslant c \| u_0 - \tilde u_0 \|_{L^2(\Omega)} \leqslant c\epsilon. \end{align*}$$

As a result, we apply triangle inequality to obtain

$$\begin{align*} \| u_\gamma(0) - u_0 \|_{L^2(\Omega)} &\leqslant \| u_0 - \tilde u_0 \|_{L^2(\Omega)} + \| u_\gamma(0) - \tilde u_\gamma(0) \|_{L^2(\Omega)} + \| \tilde u_\gamma(0) - \tilde u_0 \|_{L^2(\Omega)} \\ &\leqslant c \| u_0 - \tilde u_0 \|_{L^2(\Omega)} + c \gamma \| \tilde u_0 \|_{\dot H^{2}(\Omega)}. \end{align*}$$

Let ε be an arbitrarily small number. Using the density of $\dot H^{2}(\Omega)$ in $L^2(\Omega)$, we choose $\tilde u_0$ such that $c \| u_0 - \tilde u_0 \|_{L^2(\Omega)} \leqslant \frac{\epsilon}{2}$. Moreover, let γ₀ be the constant that $c \gamma_0 \| \tilde u_0 \|_{\dot H^{2}(\Omega)} \lt \frac{\epsilon}{2}$. Therefore, for all $\gamma\leqslant \gamma_0$, we have $\|{u}_\gamma(0) - u(0)\|_{L^2(\Omega)} \leqslant \epsilon$. Then the proof is complete.

Then we are ready to state our main theorem to show the error for the regularizing solution ${u_\gamma^\delta}(0)$.

Theorem 3.1. Let u and ${u_\gamma^\delta}$ be the solutions to the backward problem (2.9) and regularized problem (3.1), respectively. Suppose assumption 2.5 is valid. Then if $\| u_0 \|_{\dot H^{q}(\Omega)}\leqslant c$ with $q\in (0,2]$ there holds

$$\begin{equation*} \| {u_\gamma^\delta}(0) - u(0) \|_{L^2(\Omega)} \leqslant c \,\Big(\delta \gamma^{-1} + \gamma^{\frac q2}\Big). \end{equation*}$$

Moreover, for $u_0\in L^2(\Omega)$, there holds

$$\begin{equation*} \| {u_\gamma^\delta}(0) - u(0)\|_{L^2(\Omega)} \rightarrow 0\quad \text{as}~~\delta,~\gamma\rightarrow0 ~~\text{and}~~\frac\delta\gamma\rightarrow0. \end{equation*}$$

Proof. To show the error estimate, we consider the splitting

$$\begin{align*} {u_\gamma^\delta}(t) - u(t) = ({u_\gamma^\delta}(t) - {u_\gamma}(t)) + ({u_\gamma}(t) - u(t)) = \vartheta(t) + \varrho(t).\end{align*}$$

Using the solution representation (3.3) and (3.4), we have

$$\begin{align*} \| \vartheta(0) \|_{L^2(\Omega)}&\leqslant \|(\gamma I + F(T;T))^{-1} (g - g^\delta)\|_{L^2(\Omega)} \\ &\quad +\|(\gamma I + F(T;T))^{-1}\int_0^T E(T-s;T) (A(T) - A(s)) \theta(s) \,{\,\mathrm{d}}s \|_{L^2(\Omega)}. \end{align*}$$

Using lemma (3.1), we derive

$$\begin{align*} \| \vartheta(0) \|_{L^2(\Omega)}&\leqslant c \gamma^{-1} \| g - g^\delta \|_{L^2(\Omega)} +\| F(T;T)^{-1}\int_0^T E(T-s;T) (A(T) - A(s)) \theta(s) \,{\,\mathrm{d}}s \|_{L^2(\Omega)}\\ &\leqslant c \gamma^{-1} \delta +\| F(T;T)^{-1}\int_0^T E(T-s;T) (A(T) - A(s)) \theta(s) \,{\,\mathrm{d}}s \|_{L^2(\Omega)}. \end{align*}$$

Applying the argument in theorems 2.2 and 2.4, we conclude that $\| \vartheta(0) \|_{L^2(\Omega)} \leqslant c \gamma^{-1} \delta$. This estimate and lemma 3.3 lead to the desired result.

In this section, we shall propose and analyze a completely discrete scheme for solving the backward problem. To begin with, we study the semidiscrete scheme using the FEMs. The semidiscrete solution plays an important role in the analysis of completely discrete scheme.

4.1. Semidiscrete scheme for solving the problem

To begin with, we study the semidiscrete scheme using the FEMs. Let ${\{\mathcal{T}_h\}}_{0\lt h\lt1}$ be a family of shape regular and quasi-uniform partitions of the domain Ω into d-simplexes, called finite elements, with h denoting the maximum diameter of the elements. We consider the finite element space X_h defined by

$$\begin{equation} X_h =\left\{\chi\in C(\bar\Omega)\cap H_0^1: \ \chi|_{K}\in P_1(K), \,\,\,\,\forall K \in \mathcal{T}_h\right\} \end{equation} \tag{ 4.1 }$$

where $P_1(K)$ denotes the space of linear polynomials on K. Then we define the $L^2(\Omega)$ projection $P_h:L^2(\Omega)\to X_h$, by

$$\begin{align*} (P_h \psi,\chi) & =(\psi,\chi) \quad\forall \chi\in X_h,\psi\in L^2(\Omega). \end{align*}$$

Then P_h satisfies the following approximation properties [33, chapter 1]

$$\begin{align} \|P_hv-v\|_{{L^2(\Omega)}}+h\|\nabla(P_hv-v)\|_{{L^2(\Omega)}}\leqslant ch^q\|v\|_{H^q(\Omega)},\ \forall v\in \dot H^{q}(\Omega),\ q=1,2. \end{align} \tag{ 4.2 }$$

The semidiscrete standard Galerkin FEM of problem (1.1) reads: find $u_h\in X_h$ such that

$$\begin{align} &({\partial_t^\alpha} u_h(t),\chi)+(a(\cdot,t)\nabla u_h(t),\nabla \chi)=(f(\cdot,t),\chi),~~\forall \chi\in X_h,~~t\in(0,T], \nonumber\\ & \quad \text{with }u_h(0) = P_hu_0. \end{align} \tag{ 4.3 }$$

We also need a time-dependent discrete elliptic operator $A_h(t):X_h\to X_h$ by

$$\begin{equation*} (A_h(t)v_h,\chi) = (a(\cdot,t)\nabla v_h,\nabla\chi),~~\forall v_h,\chi\in X_h. \end{equation*}$$

With conditions (1.2) and (1.3), $A_h(t)$ is bounded and invertible on X_h , and problem (4.3) can be written as

$$\begin{equation} {\partial_t^\alpha} u_h+A_hu_h = P_hf,~~\forall t\in(0,T],~~u_h(0) = P_hu_0. \end{equation} \tag{ 4.4 }$$

Besides, we have the following perturbation result, which has been proved in [12, Remark 3.1].

Lemma 4.1. Under condition (1.2) and (1.3), there holds

$$\begin{align*} \|(I-A_h(t)^{-1}A_h(s)) v_h\|_{L^2(\Omega)}\leqslant c\min(1,|t-s|)\| v_h\|_{L^2(\Omega)}. \end{align*}$$

Next, we introduce a time-dependent Ritz projection operator $R_h(t):H_0^1(\Omega)\to X_h$:

$$\begin{equation} (a(\cdot,t)\nabla R_h(t)\varphi,\nabla \chi)=(a(\cdot,t)\nabla\varphi,\nabla \chi),~~\forall \varphi\in H_0^1(\Omega),~~\chi\in X_h. \end{equation} \tag{ 4.5 }$$

It is well-known that the Ritz projection satisfies the following approximation property [24, p 99]:

$$\begin{equation} \|R_h(t)v-v\|_{L^2(\Omega)}+h\|\nabla(R_h(t)v-v)\|_{L^2(\Omega)}\leqslant ch^q\|v\|_{H^q(\Omega)},~~\forall v \in \dot H^{q}(\Omega),~~q=1,2. \end{equation} \tag{ 4.6 }$$

Next, with assumption 2.3, we have an updated version of the discrete perturbation estimate.

Lemma 4.2. With conditions (1.2), (1.3) and assumption 2.3, we have for all $v_h \in X_h$

$$\begin{align*} \|(I-A_h(t)^{-1}A_h(s)) v_h\|_{L^2(\Omega)}\leqslant c\min(1,\min(t,s)^{-\kappa}|t-s|)\|v_h\|_{L^2(\Omega)},~~\forall t,s\gt1. \end{align*}$$

Proof. Let $w_h = A_h(t)^{-1}A_h(s)v_h$. Then we have $A_h(t)w_h = A_h(s)v_h$ and hence

$$\begin{equation*} (a(\cdot,t)\nabla w_h,\nabla\chi_h) = (a(\cdot,s)\nabla v_h,\nabla \chi_h),~~\forall \chi_h\in X_h. \end{equation*}$$

This further implies the relation

$$\begin{equation*} (a(\cdot,t)\nabla (v_h-w_h),\nabla \chi_h) = ((a(\cdot,t)-a(\cdot,s))\nabla v_h,\nabla \chi_h),~~\forall \chi\in X_h. \end{equation*}$$

Let φ be the weak solution to the following elliptic problem:

$$\begin{equation*} (a(\cdot,t)\nabla \phi,\nabla \chi) = ((a(\cdot,t)-a(\cdot,s))\nabla v_h,\nabla \chi),~~\forall \chi\in \dot H^{1}(\Omega). \end{equation*}$$

Then Lax–Milgram lemma and assumption 2.3 implies the following a priori estimate

$$\begin{equation*} \|\phi\|_{\dot H^{1}(\Omega)}\leqslant c\|(a(\cdot,t)-a(\cdot,s))\nabla v_h\|_{L^2(\Omega)}\leqslant c\min(1,\min(t,s)^{-\kappa}|t-s|)\|v_h\|_{\dot H^{1}(\Omega)}. \end{equation*}$$

Using the fact that $w_h - v_h = R_h(t) \phi$, the approximation property (4.6), and the inverse inequality, we derive

$$\begin{align*} \|w_h-v_h-\phi\|_{L^2(\Omega)}&\leqslant ch\|\varphi\|_{\dot H^{1}(\Omega)}\leqslant ch\min(1,\min(t,s)^{-\kappa}|t-s|)\|v_h\|_{\dot H^{1}(\Omega)} \\ &\leqslant c\min(1,\min(t,s)^{-\kappa}|t-s|)\|v_h\|_{L^2(\Omega)}. \end{align*}$$

According triangle inequality we have

$$\begin{align*} \|w_h-v_h\|_{L^2(\Omega)}\leqslant c\min(1,\min(t,s)^{-\kappa}|t-s|)\|v_h\|_{L^2(\Omega)}+\|\phi\|_{L^2(\Omega)}. \end{align*}$$

Next, we apply the duality argument to derive a bound for $\|\phi\|_{L^2(\Omega)}$. Let $\xi \in \dot H^{2}(\Omega)$ be the function such that $A(t) \xi = \phi$. Then

$$\begin{align*} \| \phi \|_{L^2(\Omega)}^2 &= |(a(\cdot,t)\nabla\phi,\nabla\xi)| = |((a(\cdot,t)-a(\cdot,s))\nabla v_h,\nabla\xi)|\\ &\leqslant |(v_h,(a(\cdot,t)-a(\cdot,s)) \Delta\xi)| + |(v_h,\nabla (a(\cdot,t)-a(\cdot,s)) \cdot\nabla \xi)|\\ &\leqslant c\min(1,\min(t,s)^{-\kappa}|t-s|) \|v_h\|_{L^2(\Omega)}\|\xi\|_{\dot H^{2}(\Omega)}\\ &\leqslant c\min(1,\min(t,s)^{-\kappa}|t-s|) \|v_h\|_{L^2(\Omega)}\|\phi\|_{L^2(\Omega)}. \end{align*}$$

This completes the proof of the lemma.

Next we derive some semidiscrete solution representation analogue to (2.6), that is given any $t_*\in(0,T]$,

$$\begin{equation} u_h(t) = F_h(t;t_*)u_h(0)+\int_0^t E_h(t-s;t_*)(P_h\,f(s)+(A_h(t_*)-A_h(s))u_h(s))ds \end{equation} \tag{ 4.7 }$$

where the solution operators $F_h(t;t_*)$ and $E_h(t;t_*)$ can be written as

$$\begin{align} F_h(t;t_*) &= \frac{1}{2\pi\mathrm{i}} \int_{\Gamma_{\theta,\kappa}} e^{zt} z^{\alpha-1} (z^\alpha + A_h(t_*))^{-1} \,{\,\mathrm{d}}z,~~\text{and} \nonumber\\ E_h(t;t_*) &= \frac{1}{2\pi\mathrm{i}} \int_{\Gamma_{\theta,\kappa}} e^{zt} (z^\alpha +A_h(t_*))^{-1}\,{\,\mathrm{d}}z. \end{align} \tag{ 4.8 }$$

For any fixed $t_*$, the discrete operators $F_h(t;t_*)$ and $E_h(t;t_*)$ satisfy the following smoothing property, whose proof is identical to that of lemma 2.1.

Lemma 4.3. Let $F_h(t;t_*)$ and $E_h(t;t_*)$ be the discrete solution operators defined in (4.8) for any $t_*\in[0,T]$. Then they satisfy the following properties for all t > 0 and $v_h \in X_h$

• (a)  
  $\|A_h(t_*) F_h (t;t_*)v_h\|_{L^2(\Omega)} + t^{1-(2-k)\alpha} \| A_h(t_*)^k E_h (t;t_*)v_h \|_{L^2(\Omega)} \leqslant c t^{-\alpha} \|v_h\|_{L^2(\Omega)}$, with $k = 1,2$;
• (b)  
  $\|F_h(t;t_*)v_h\|_{L^2(\Omega)}+ t^{1-\alpha}\|E_h(t;t_*)v_h\|_{L^2(\Omega)} \leqslant c \min(1, t^{-\alpha}) \|v_h\|_{L^2(\Omega)}$;
• (c)  
  $\|F_h(t; t_*)^{-1} v_h \|_{L^2(\Omega)} \leqslant c (1+ t^\alpha) \| A_h(t_*) v_h \|_{L^2(\Omega)}$.

The constants in all above estimates are uniform in t, but they are only dependent of $t_*$ and T.

Analogue to lemma 3.1, we have the following result.

Lemma 4.4. Let $F_h(t;t_*)$ be the discrete solution operator defined in (4.8). For all $0\lt t\leqslant T$, $t_*\in(0,T]$ and $v_h\in X_h$, we have

$$\begin{align*} &\| (\gamma I+F_h(T;T))^{-1}v_h \| \leqslant c \gamma^{-1}\| v_h \|_{L^2(\Omega)}~~\text{and}\\ & \|F_h(T;T)(\gamma I+F_h(T;T))^{-1}v_h\|_{L^2(\Omega)}\leqslant c\|v_h\|_{L^2(\Omega)}, \end{align*}$$

where the constant c is independent of t, γ and h.

The following Lemma provides an error estimate for the semidiscrete error of the direct problem, see [12, theorem 3.2] for a detailed proof.

Lemma 4.5. Let u and u_h be the solutions to (1.1) and (4.3) respectively. If $u_0\in L^2$ and $f\equiv 0$, then there holds that

$$\begin{equation*} \|(u_h-u)(t)\|_{L^2(\Omega)}\leqslant ch^2 t^{-\alpha}\|u_0\|_{L^2(\Omega)}~~\text{for all}~~t\in(0,T], \end{equation*}$$

where the constant c is independent of t and h.

After proposing many results about solving direct problem, we shall propose a semidiscrete scheme for solving the backward problem.

We apply the regularized semidiscrete scheme: find $u_{\gamma,h}(t)\in X_h$ such that

$$\begin{equation} {\partial_t^\alpha} u_{\gamma,h}(t)+A_h(t)u_{\gamma,h}(t) = 0,~~0\lt t\leqslant T,\qquad \gamma u_{\gamma,h}(t) + u_{\gamma,h}(T) = P_hg. \end{equation} \tag{ 4.9 }$$

Then analogue to (3.3) we have

$$\begin{equation} u_{\gamma,h}(0)=(\gamma I+F_h(T;T))^{-1} \left[P_hg-\int_0^TE_h(T-s;T)(A_h(T)-A_h(s))u_{\gamma,h}(s)ds\right]. \end{equation} \tag{ 4.10 }$$

Next we shall derive a preliminary estimate for the proof of the semidiscrete error $u_{\gamma,h}-{u_\gamma}$.

Lemma 4.6. Let ${u_\gamma}(t)$ be the solution to the backward regularized problem (4.9). Then fix any $t_*\in (0,T]$ there holds that

$$\begin{equation*} \left\|\int_0^{t_*} E_h(t_*-s;t_*)A_h(s)(R_h(s)-P_h){u_\gamma}(s)ds\right\|_{L^2(\Omega)}\leqslant ch^2 {\max\{t_*^{-\alpha},t_*^{1-\alpha}\}}\|{u_\gamma}(0)\|_{L^2(\Omega)}. \end{equation*}$$

The constant c is independent of t and $t_*$.

Proof. Let ϕ_h be the solution to the following semidiscrete problem

$$\begin{equation*} {\partial_t^\alpha} \varphi_h(t)+A_h(t)\varphi_h(t) =0,\qquad \varphi_h(0) = P_h {u_\gamma}(0). \end{equation*}$$

Lemma 4.5 implies that

$$\begin{equation*} \|(\varphi_h-{u_\gamma})(t)\|_{L^2(\Omega)}\leqslant ch^2 t^{-\alpha} \|{u_\gamma}(0)\|_{L^2(\Omega)}. \end{equation*}$$

Then we consider the splitting

$$\begin{equation*} (\varphi_h-{u_\gamma})(t)=(\varphi_h-P_h{u_\gamma})(t)+(P_h{u_\gamma}-{u_\gamma})(t):=\zeta_h(t)+\rho(t). \end{equation*}$$

The approximation property (4.2) and the regularity estimate in theorem 2.1 give that

$$\begin{equation*} \|\rho(t)\|_{L^2(\Omega)}\leqslant ch^2\|{u_\gamma}(t)\|_{{\dot H^{2}(\Omega)}}\leqslant ch^2 t^{-\alpha} \|{u_\gamma}(0)\|_{L^2(\Omega)}. \end{equation*}$$

Then by triangle's inequality, we obtain

$$\begin{equation} \|\zeta_h(t)\|_{L^2(\Omega)}\leqslant \| \rho(t) \|_{L^2(\Omega)} + \| (\varphi_h-{u_\gamma})(t) \|_{L^2(\Omega)} \leqslant ch^2 t^{-\alpha} \|{u_\gamma}(0)\|_{L^2(\Omega)}. \end{equation} \tag{ 4.11 }$$

Meanwhile notice that

$$\begin{equation*} {\partial_t^\alpha} \zeta_h(t)+A_h(t)\zeta_h(t)=A_h(t)(R_h(t)-P_h){u_\gamma}(t),~~T\geqslant t\gt0,\qquad \zeta(0) = 0. \end{equation*}$$

Then for any $t_* \in (0,T]$, $\zeta_h(t_*)$ could be written as

$$\begin{align*} \zeta_h(t_*) &= \int_0^{t_*} E_h(t_*-s;t_*)A_h(s)(R_h(s)-P_h){u_\gamma}(s)ds\\ & \quad +\int_0^{t_*} E_h(t_*-s;t_*)(A_h(t_*)-A_h(s))\zeta_h(s)ds. \end{align*}$$

We apply lemmas 4.2 and 4.3, and the estimate (4.11) to derive

$$\begin{equation*} \begin{aligned} &\quad\|\int_0^t E_h(t-s;t_*)A_h(s)(R_h(s)-P_h){u_\gamma}(s)ds\|_{L^2(\Omega)} \\ &\leqslant \|\zeta_h(t)\|_{L^2(\Omega)} + c\int_0^{t_*} \|\zeta_h(s)\|_{L^2(\Omega)} \,{\,\mathrm{d}}s \leqslant c (t_*^{-\alpha}+t_*^{1-\alpha})h^2\|{u_\gamma}(0)\|_{L^2(\Omega)} \\ & \leqslant {c h^2 \max\{t_*^{-\alpha},t_*^{1-\alpha}\}\|{u_\gamma}(0)\|_{L^2(\Omega)}}. \end{aligned} \end{equation*}$$

Next, we state a key lemma providing an estimate for the discretization error $u_{\gamma,h}-{u_\gamma}$.

Lemma 4.7. Let ${u_\gamma}(t)$, $u_{\gamma,h}(t)$ be the solutions to problem (3.2) and (4.9) respectively. Suppose assumption 2.5 is valid. Then there holds

$$\begin{align*} \| u_{\gamma,h}(0) - {u_\gamma}(0) \| \leqslant c h^2\gamma^{-1}\|u_0\|_{L^2(\Omega)}, \end{align*}$$

where the constant c is independent on γ, h and t.

Proof. We use the splitting

$$\begin{equation*} (u_{\gamma,h}-{u_\gamma})(0) = (u_{\gamma,h}-P_h{u_\gamma})(0)+(P_h{u_\gamma}-{u_\gamma})(0):=\zeta_h(0)+\rho(0). \end{equation*}$$

From the approximation property (4.2) and lemma 3.2, we obtain

$$\begin{equation*} \|\rho(0)\|_{L^2(\Omega)}\leqslant ch^2\|{u_\gamma}(0)\|_{\dot H^{2}(\Omega)}\leqslant ch^2 \gamma^{-1}\|u_0\|_{L^2(\Omega)}. \end{equation*}$$

Now we turn to the bound of $\zeta_h(t)$. Using the fact $A_h(t)R_h(t)v = P_hA(t)v$, we observe that

$$\begin{align} &{\partial_t^\alpha} \zeta_h(t)+ A_h(t)\zeta_h(t) = A_h(t)(R_h(t)-P_h){u_\gamma}(t)~~\text{for}~~t\in(0,T], \nonumber\\ & \quad\quad\text{with}~~\gamma \zeta_h(0)+\zeta_h(T) = 0. \end{align} \tag{ 4.12 }$$

For any $t_*\in (0,T]$, we have the solution representation from (4.7) that

$$\begin{align*} \zeta_h(t) &= F_h(t;t_*)\zeta_h(0)+\int_0^t E_h(t-s;t_*)A_h(s)(R_h(s)-P_h){u_\gamma}(s)\, {\,\mathrm{d}}s \\ &\qquad +\int_0^t E_h(t-s;t_*)(A_h(t_*)-A_h(s))\zeta_h(s)\,{\,\mathrm{d}}s. \end{align*}$$

Then with $t = t_* = T$ we apply $\gamma \zeta_h(0)+\zeta_h(T) = 0$ to derive

$$\begin{equation*} \begin{aligned} \zeta_h(0) &= (\gamma I +F_h(T;T))^{-1}\int_0^TE_h(T-s;T)(A_h(s)(P_h-R_h(s)){u_\gamma}(s)\,{\,\mathrm{d}}s \\ &\qquad -(\gamma I+F_h(T;T))^{-1}\int_0^T E_h(T-s;T)(A_h(T)-A_h(s))\zeta(s)\,{\,\mathrm{d}}s. \end{aligned} \end{equation*}$$

Now we apply lemmas 4.4 and 4.6 to obtain

$$\begin{align*} \|\zeta_h(0)\|_{L^2(\Omega)}&\leqslant c\gamma^{-1}\int_0^T \|E_h(T-s;T)A_h(s)(P_h-R_h(s)){u_\gamma}(s)\|_{L^2(\Omega)}{\,\mathrm{d}}s\\ &+\|F_h(T;T)^{-1}\int_0^TE_h(T-s;T)(A_h(T)-A_h(s))\zeta(s) {\,\mathrm{d}}s\|_{L^2(\Omega)}\\ &\leqslant {c_T}h^2\gamma^{-1} \|{u_\gamma}(0)\|_{L^2(\Omega)} +c(1+T^{\,\alpha})\int_0^T \|A_h(T;T)E_h(T-s;T)\\ & \quad (A_h(T)-A_h(s))\zeta_h(s)\|_{L^2(\Omega)} {\,\mathrm{d}}s. \end{align*}$$

Next, we split $\zeta_h(s)$ into homogeneous part and inhomogeneous part. Let $\zeta_h(t) : = \zeta_1(t)+\zeta_2(t)$ where

$$\begin{align*} {\partial_t^\alpha} \zeta_1(t)+A_h(t)\zeta_1(t) &= 0~~\text{for}~~t\in(0,T], \quad\text{with}~~\zeta_1(0) = \zeta_h(0),\\ {\partial_t^\alpha} \zeta_2(t)+A_h(t)\zeta_2(t) &= A_h(t)(R_h(t)-P_h){u_\gamma}(t)~~\text{for}~~t\in(0,T], \quad\text{with}~~\zeta_2(0) = 0. \end{align*}$$

First of all, we fixed $t_*\in (0,T]$ and apply the solution representation in (4.7) and lemmas 4.1, 4.3 and 4.6, and hence derive

$$\begin{align*} \|\zeta_2(t_*)\|_{L^2(\Omega)}&\leqslant \int_0^{t_*} \|E_h(t_*-s;t_*)A_h(s)(R_h(s)-P_h){u_\gamma}(s)\|_{L^2(\Omega)}{\,\mathrm{d}}s\\ &+ \int_0^{t_*} \|E_h(t_*-s;t_*)(A_h(t_*)-A_h(s))\zeta_2(s)\|_{L^2(\Omega)}{\,\mathrm{d}}s\\ &\leqslant ct_*^{-\alpha}h^2\|{u_\gamma}(0)\|_{L^2(\Omega)} + \int_0^{t_*} \|\zeta_2(s)\|_{L^2(\Omega)}{\,\mathrm{d}}s. \end{align*}$$

Then Gronwall's inequality leads to

$$\begin{equation} \|\zeta_2(t)\|_{L^2(\Omega)}\leqslant ch^2e^{ct}t^{-\alpha}\|{u_\gamma}(0)\|_{L^2(\Omega)}. \end{equation} \tag{ 4.13 }$$

For $\zeta_1(t)$, we apply the similar argument in lemma 2.3 to obtain

$$\begin{align*} &\| \zeta_1(t) \|_{L^2(\Omega)} \leqslant c \min(1,t^{-\alpha}) \| {\zeta_h(0)} \|_{L^2(\Omega)}~~\text{and}\\ & \|A_h(T) \zeta_1(t)\|_{L^2(\Omega)} \leqslant c \,e^{ct} t^{-\alpha} \| {\zeta_h(0)} \|_{L^2(\Omega)} \quad \text{for all}~~t\gt0. \end{align*}$$

All the positive constants c in above estimates are independent of t and T. As a result, we have

$$\begin{align*} \|\zeta_h(0)\|_{L^2(\Omega)}&\leqslant ch^2\gamma^{-1} \|{u_\gamma}(0)\|_{L^2(\Omega)} +\sum_{i=1}^2c(1+T^{\,\alpha})\int_0^T \|A_h(T;T)E_h(T-s;T)(A_h(T)\\ &\quad -A_h(s))\zeta_i(s) {\,\mathrm{d}}s\|_{L^2(\Omega)}\\ &\leqslant ch^2\gamma^{-1}\|{u_\gamma}(0)\|_{L^2(\Omega)}+c(1+T^{\,\alpha})\int_0^T \|A_h(T)E_h(T-s;T)\\ & \quad (A_h(T)-A_h(s))\zeta_1(s)\|_{L^2(\Omega)} {\,\mathrm{d}}s. \end{align*}$$

Applying the argument in theorems 2.2 and 2.4, we conclude that $\|\zeta_h(0)\|_{L^2(\Omega)} \leqslant c h^2 \gamma^{-1} \|{u_\gamma}(0)\|_{L^2(\Omega)}$. This completes the proof of the lemma.

4.2. Fully discrete scheme and error analysis

To begin with, we introduce the fully discrete scheme for the direct problem. We divide the time interval $[0,T]$ into a uniform grid, with $t_n = n\tau$, $n = 0,\ldots,N$, and $\tau = T/N$ being the time step size. Then we approximate the fractional derivative by using the CQ-BE (with $\varphi^{\,j} = \varphi(t_j)$) [10, 23]:

$$\begin{equation*} \bar\partial_\tau^\alpha \varphi^n = \sum_{j=0}^n \omega_{n-j}^{(\alpha)} (\varphi^{\,j} - \varphi^0), \quad\mbox{with }~\omega_j^{(\alpha)}= (-1)^{\,j}\frac{\Gamma(\alpha+1)}{\Gamma(\alpha-j+1)\Gamma(\,j+1)}. \end{equation*}$$

The fully discrete scheme for problem (4.4) reads: find ${U_h^n}\in X_h$ such that

$$\begin{equation} \bar\partial_\tau^\alpha U_h^n +A_h(t_n)U_h^n= P_h\, f(t_n),\quad n=1,2,\ldots,N,\quad \text{with}~~U_h^0 = P_hu_0. \end{equation} \tag{ 4.14 }$$

By means of Laplace transform and perturbation argument, with $1\leqslant n_*\leqslant N$, the fully discrete solution Uⁿ can be written as [12, 41]

$$\begin{equation} U_h^n = F_{h, \tau}^n(n_*) U_h^0+ \tau \sum_{k=1}^n E_{h, \tau}^{n-k}(n_*)P_h\,f(t_k)+\tau\sum_{k=1}^nE_{h, \tau}^{n-k}(n_*)(A_h(t_{n_*})-A_h(t_k)) U_h^k \end{equation} \tag{ 4.15 }$$

with $n = 1,2,\cdots,N$. Here the fully discrete operators $F_{h, \tau}^n(n_*)$ and $E_{h, \tau}^n(n_*)$ are defined by

$$\begin{equation} \begin{aligned} F_{h, \tau}^n(n_*) &= \frac{1}{2\pi\mathrm{i}}\int_{\Gamma_{\theta,\sigma}^\tau } e^{zt_n} {e^{-z\tau}} \delta_\tau(e^{-z\tau})^{\alpha-1}({\delta_\tau(e^{-z\tau})^\alpha}+A_h(t_{n_*}))^{-1}\,{\,\mathrm{d}}z,\\ E_{h, \tau}^n(n_*) &= \frac{1}{2\pi\mathrm{i}}\int_{\Gamma_{\theta,\sigma}^\tau } e^{zt_n} ({\delta_\tau(e^{-z\tau})^\alpha}+A_h(t_{n_*}))^{-1}\,{\,\mathrm{d}}z , \end{aligned} \end{equation} \tag{ 4.16 }$$

with $\delta_\tau(\xi) = (1-\xi)/\tau$ and the contour $\Gamma_{\theta,\sigma}^\tau : = \{z\in \Gamma_{\theta,\sigma}:|\Im(z)|\leqslant {\pi}/{\tau} \}$ where $\theta\in(\pi/2,\pi)$ is close to $\pi/2$. (oriented with an increasing imaginary part). The next lemma gives elementary properties of the kernel $\delta_\tau(e^{-z\tau})$. The detailed proof has been given in [10, lemma B.1].

Lemma 4.8. For a fixed $\theta^{^{\prime}}\in(\pi/2,\pi/\alpha)$, there exists $\theta\in(\pi/2,\pi)$ positive constants $c,c_1,c_2$ $($independent of τ) such that for all $z\in \Gamma_{\theta,\sigma}^\tau$

$$\begin{equation*} \begin{aligned} & c_1|z|\leqslant |\delta_\tau(e^{-z\tau})|\leqslant c_2|z|,\qquad \delta_\tau(e^{-z\tau})\in \Sigma_{\theta^{^{\prime}}}, \\ & |\delta_\tau(e^{-z\tau})-z|\leqslant c\tau |z|^{2},\qquad |\delta_\tau(e^{-z\tau})^\alpha-z^\alpha|\leqslant c\tau |z|^{1+\alpha}. \end{aligned} \end{equation*}$$

The next lemma provides some approximation properties of solution operators $F_{h,\tau}^n(n_*)$ and $E_{h,\tau}^n(n_*)$. See [40, lemma 4.2] and [9, theorem 3.5] for the proof of the first estimate, and [12, lemma 4.5] for the second estimate.

Lemma 4.9. For the operator $F_h^\tau$ and $E_h^\tau$ defined in (4.16), we have

$$\begin{align*} \| A_h(t_{n_*})^\beta (F_{h,\tau}^n(n_*)- F_{h}(t_n; t_{n_*}))\|_{L^2(\Omega)}\leqslant c\tau t_n^{-1-\beta\alpha},\\ \Big\|\tau A_h^\beta E_{h,\tau}^{{n_*}-k}(n_*)-\int_{t_{k-1}}^{t_k} A_h^\beta(t_{n_*}) E_h(t_{n_*}-s;t_{n_*}))ds\Big\|\leqslant c\tau^2(t_{n_*}-t_k+\tau)^{-(2-(1-\beta)\alpha)} \end{align*}$$

for any $\beta\in[0,1]$.

Note that the solution operators $F_{h,\tau}^n(n_*)$ and $E_{h,\tau}^n(n_*)$ satisfy the following smoothing properties, whose proof is identical to that of lemma 2.1.

Lemma 4.10. Let $F_h^\tau(n;n)$ and $E_h^\tau(n;n)$ be the operators defined in (4.16). Then they satisfy the following properties for any $n\geqslant 1$ and $v_h\in X_h$,

• (a)  
  $\| A_h(t_*) F_{h,\tau}^n (n_*)v_h\|_{L^2(\Omega)} + t_n^{1-(2-k)\alpha} \| A_h(t_*)^k E_{h,\tau}^n (n_*)v_h \|_{L^2(\Omega)} \leqslant c t_{n+1}^{-\alpha} \|v_h\|_{L^2(\Omega)},$ $k = 1,2$;
• (b)  
  $\|F_{h,\tau}^n (n_*)v_h\|_{L^2(\Omega)}+ t_n^{1-\alpha}\|E_{h,\tau}^n (n_*)v_h\|_{L^2(\Omega)} \leqslant c \min(1, t_n^{-\alpha}) \|v_h\|_{L^2(\Omega)}$;
• (c)  
  $\| F_{h,\tau}^n (n_*)^{-1} v_h \|_{L^2(\Omega)} \leqslant c (1+ t_n^\alpha) \| A_h(t_*) v _h\|_{L^2(\Omega)}$.

Next we introduce some a priori estimate for the discrete solution $U_h^n$ in (4.15), analogue to lemma 2.3 for the continuous problem. We provide the proof in appendix B for completeness.

Lemma 4.11. Let $U_h^n$ be the solution to (4.14), then we have the following a-priori estimate ($f\equiv 0$)

$$\begin{align*} &\|U_h^n\|_{L^2(\Omega)}\leqslant c\min(1,t_n^{-\alpha})\|U_h^0\|_{L^2(\Omega)} ~~\text{and}\\&\|A_h(T)U_h^n\|_{L^2(\Omega)}\leqslant ce^{ct_n}t_n^{-\alpha}\|U_h^0\|_{L^2(\Omega)}~~\text{for}~~n\geqslant1. \end{align*}$$

Moreover, for any $\epsilon\in(0,1/\alpha-1)$, there holds

$$\begin{equation*} \|A_h(T)U_h^n\|_{L^2(\Omega)} \leqslant c t_n ^{-(1-\epsilon)\alpha}\|U_h^0\|_{L^2(\Omega)}~~\text{for all}~~n\geqslant1. \end{equation*}$$

All the constants in above estimates are independent of h, n, N, τ and T.

Now we introduce the fully discrete scheme for solving the backward problem: find $U_{h,\gamma}^{n,\delta} \in X_h$ for $n = 1,\ldots,N$ such that

$$\begin{equation} \bar\partial_\tau^\alpha U_{h,\gamma}^{n,\delta} +A_h(t_n)U_{h,\gamma}^{n,\delta} = 0~~\text{for all}~~1\leqslant n\leqslant N, \quad \text{with}~~ \gamma U_{h,\gamma}^{0,\delta} + U_{h,\gamma}^{N,\delta} = P_h g_\delta. \end{equation} \tag{ 4.17 }$$

Then $U_{h,\gamma}^{n,\delta}$ can be written as

$$\begin{equation} U_{h,\gamma}^{0,\delta} = (\gamma I+F_{h,\tau}^N(N))^{-1}\Big[P_hg_\delta-\tau \sum_{k=1}^N {E_{h,\tau}^{N-k}}(N)(A_h(T)-A_h(t_k))U_{h,\gamma}^{k,\delta}\Big]. \end{equation} \tag{ 4.18 }$$

The following lemma provides a useful estimate of the discrete operator $(\gamma I+F_{h,\tau}^N(N))^{-1}$; see a detailed proof in [42, lemma 4.4].

Lemma 4.12. Let $F_h^\tau(n;n_*)$ and $E_h^\tau(n;n_*)$ be the operators defined in (4.16). Then there holds

$$\begin{align*} & \|(\gamma I+F_h^\tau(N;N))^{-1}v_h\|_{L^2(\Omega)}\leqslant c\gamma^{-1}\|v_h\|_{L^2(\Omega)}~~\text{and}\\ & \|F_h^\tau(N;N)(\gamma I+F_h^\tau(N;N))^{-1} v_h\|_{L^ 2(\Omega)}\leqslant c \end{align*}$$

where c is uniform in T, h, τ and γ.

To show the error between $U_{h,\gamma}^{N,\delta}$ and u₀, we introduce an auxiliary function $\bar U_{h,\gamma}^n \in X_h$ such that

$$\begin{equation} \bar\partial_\tau^\alpha \bar U_{h,\gamma}^n +A_h(t_n)\bar U_{h,\gamma}^n = P_hf(t_n) ~~\text{for all}~~1\leqslant n\leqslant N, \quad \text{with}~~\bar U_{h,\gamma}^n = u_{\gamma,h}(0). \end{equation} \tag{ 4.19 }$$

Then we have the following error estimate for the direct problem, according to [12, theorem 4.1].

Lemma 4.13. Let $u_{\gamma,h}(t)$ and $\bar U_{h,\gamma}^n$ be the solution to (4.4) and (4.19) with $f\equiv 0$, then we have

$$\begin{equation*} \| A_h(0) (\bar U_{h,\gamma}^n-u_{\gamma,h}(t_n))\|_{L^2(\Omega)} \leqslant c \tau \log (n+1) \max(t_n^{-\alpha-1}, t_n^{-\alpha} ) \|u_{\gamma,h}(0)\|_{L^2(\Omega)}. \end{equation*}$$

Proof. Let $e_n = \bar U_{h,\gamma}^n-u_{\gamma,h}(t_n)$. First of all, we recall [12, theorem 4.1] that

$$\begin{equation} \| e_n \|_{L^2(\Omega)} \leqslant c\tau t_n^{-1}\log(n+1) \|u_{\gamma,h}(0)\|_{L^2(\Omega)}. \end{equation} \tag{ 4.20 }$$

We then use the solution representation (4.15) to obtain

$$\begin{equation*} \bar U_{h,\gamma}^n = F_{h,\tau}^{n}(n_*) u_{\gamma,h}(0) + \tau \sum_{k=1}^n E_{h,\tau}^{n-k}(n_*)(A_h(t_{n_*})-A_h(t_k)) \bar U_{h,\gamma}^k. \end{equation*}$$

Then by means of (4.7), we have for fixed $t_{n_*}$

$$\begin{align*} A_h(t_{n_*})e_{n_*} &= A_h(t_{n_*})(F_{h,\tau}^n(n_*)-F_h(t_n;t_{n_*}))u_{\gamma,h}(0) \\ & \quad + \tau\sum_{k=1}^{n_*} A_h(t_{n_*})E_{h,\tau}^{n_*-k}(n_*)(A_h(t_{n_*})-A_h(t_k)) \bar U_{h,\gamma}^k \\ &\quad - \int_0^{t_{n_*}}A_h(t_{n_*}) E_h(t_{n_*}-s;t_{n_*})(A_h(t_{n_*})-A_h(s))u_{\gamma,h}(s) \, {\,\mathrm{d}}s\\ & \quad +\tau\sum_{k=1}^{n_*}A_h(t_{n_*})E_{h,\tau}^{n_*-k}(n_*)(A_h(t_{n_*})-A_h(t_k))(\bar U_{h,\gamma}^k-u_{\gamma,h}(t_k))\\ & = I_1+I_2+I_3. \end{align*}$$

Lemma 4.9 immediately implies the bound for I₁:

$$\begin{align*}\|I_1\|_{L^2(\Omega)}\leqslant c\tau t_{n_*}^{-1-\alpha} \| u_{\gamma,h}(0)\|_{L^2(\Omega)}.\end{align*}$$

A slightly modification of [12, lemma 4.4] leads to a bound for I₂. In particular, we observe

$$\begin{align*} I_2 &= \sum_{k=1}^{n_*} \Big[\tau A_h(t_{n_*})E_{h,\tau}^{n_*-k}(n_*)(A_h(t_{n_*})-A_h(t_k))\\ & \quad - \int_{t_{k-1}}^{t_{k}}A_h(t_{n_*}) E_h(t_{n_*}-s;t_{n_*})(A_h(t_{n_*})-A_h(s)){\,\mathrm{d}}s\Big]\bar U_{h,\gamma}^k\\ & + \sum_{k=1}^{n_*} \int_{t_{k-1}}^{t_k}A_h(t_{n_*}) E_h(t_{n_*}-s;t_{n_*})(A_h(t_{n_*})-A_h(s)) \, {\,\mathrm{d}}s\, e_k\\ & + \sum_{k=1}^{n_*} \int_{t_{k-1}}^{t_k}A_h(t_{n_*}) E_h(t_{n_*}-s;t_{n_*})(A_h(t_{n_*})-A_h(s))\\ & \quad (u_{\gamma,h}(t_k)-u_{\gamma,h}(s)) \, {\,\mathrm{d}}s:= I_{2,1}+I_{2,2}+I_{2,3}. \end{align*}$$

For $I_{2,1}$, by means of lemma 4.9 with β = 1, 4.1 and 4.10 (a) with the solution representation (4.15), we arrive at

$$\begin{align*}\begin{array}{l} {I_{2,1}}{_{{L^2}(\Omega )}}\sum\limits_{k = 1}^{{n_*}} {\left\| {\left[ {\tau {A_h}({t_{{n_*}}})E_{h,\tau }^{{n_*} - k}({n_*})(I - {A_h}({t_k}){A_h}{{({t_{{n_*}}})}^{ - 1}})} \right.} \right.} \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, - \int_{{t_{k - 1}}}^{{t_k}} {{A_h}} ({t_{{n_*}}}){E_h}({t_{{n_*}}} - s;{t_{{n_*}}})(I - {A_h}(s){A_h}{({t_{{n_*}}})^{ - 1}})\left. {\left. {ds} \right]} \right\|\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{\mkern 1mu} {A_h}({t_{{n_*}}})\bar U_{h,\gamma }^k{_{{L^2}(\Omega )}}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,c\sum\limits_{k = 1}^{{n_*}} {{\tau ^2}} {({t_{{n_*}}} - {t_k} + \tau )^{ - 1}}t_k^{ - \alpha }{u_{\gamma ,h}}(0){_{{L^2}(\Omega )}}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,c\tau \log ({n_*} + 1)t_n^{ - \alpha }{u_{\gamma ,h}}(0){_{{L^2}(\Omega )}}. \end{array} \end{align*}$$

For $I_{2,2}$ we apply lemma 4.3 (a) with k = 2, lemma 4.1 and a priori estimate (4.20) to derive

$$\begin{align*} \|I_{2,2}\|_{L^2(\Omega)}&\leqslant c\tau t_{n_*}^{-\alpha-1}\log(n_*+1)\|u_{\gamma,h}(0)\|_{L^2(\Omega)}. \end{align*}$$

Last, for the erm $I_{2,3}$, we denote

$$\begin{equation*} Q_k = \int_{t_{k-1}}^{t_k}A_h(t_{n_*}) E_h(t_{n_*}-s;t_{n_*})(A_h(t_{n_*})-A_h(s))(u_{\gamma,h}(t_k)-u_{\gamma,h}(s)). \end{equation*}$$

For k = 1, we apply lemmas 4.1 and 4.3 to derive the bound

$$\begin{align*} \|Q_1\|_{L^2(\Omega)}&\leqslant \|\int_0^\tau A_h(t_{n_*}) E_h(t_{n_*}-s;t_{n_*})(A_h(t_{n_*})-A_h(s))u_{\gamma,h}(\tau)\,{\,\mathrm{d}}s\|_{L^2(\Omega)}\\ &\qquad +\|\int_0^\tau A_h(t_{n_*}) E_h(t_{n_*}-s;t_{n_*})(A_h(t_{n_*})-A_h(s))u_{\gamma,h}(s)\,{\,\mathrm{d}}s\|_{L^2(\Omega)}\\ & \leqslant c\int_0^\tau (t_{n_*}-s)^{-\alpha}{\,\mathrm{d}}s\|u_{\gamma,h}(0)\|_{L^2(\Omega)} \leqslant c\tau t_{n_*}^{-\alpha} \|u_{\gamma,h}(0)\|_{L^2(\Omega)}. \end{align*}$$

Meanwhile, for $k\geqslant 2$, there holds that

$$\begin{align*} Q_k = \int_{t_{k-1}}^{t_k} A_h(t_{n_*}) E_h(t_{n_*}-s;t_{n_*})(A_h(t_{n_*})-A_h(s))\int_{s}^{t_k} u_{\gamma,h}^{^{\prime}}(\xi){\,\mathrm{d}}\xi\,{\,\mathrm{d}}s. \end{align*}$$

The discrete analogue to theorem 2.1 (a) (see detail proof in [12, theorem 2.3(i)]), $u_{\gamma,h}^{^{\prime}}(t)$ can be bounded by

$$\begin{equation} \|u_{\gamma,h}^{^{\prime}}(t)\|_{L^2(\Omega)}\leqslant ct^{-1}\|u_{\gamma,h}(0)\|_{L^2(\Omega)}. \end{equation} \tag{ 4.21 }$$

Then by lemmas 4.1, 4.3 and regularity estimate (4.21) there holds

$$\begin{align*} \|Q_k\|_{L^2(\Omega)}&\leqslant c \int_{t_{k-1}}^{t_k} \|A_h(t_{n_*}) E_h(t_{n_*}-s;t_{n_*})(A_h(t_{n_*})\\ & \quad -A_h(s))\|\,\int_{s}^{t_k} \xi^{-1}{\,\mathrm{d}}\xi\,{\,\mathrm{d}}s\|u_{\gamma,h}(0)\|_{L^2(\Omega)}\\ &\leqslant c\int_{t_{k-1}}^{t_k} (t_{n_*}-s)^{-\alpha}\int_s^{t_k} \xi^{-1}{\,\mathrm{d}}\xi\,{\,\mathrm{d}}s \|u_{\gamma,h}(0)\|_{L^2(\Omega)}\\ &\leqslant c \tau \int_{t_{k-1}}^{t_k} (t_{n_*}-s)^{-\alpha}s^{-1}{\,\mathrm{d}}s \|u_{\gamma,h}(0)\|_{L^2(\Omega)}. \end{align*}$$

Summing those terms from k = 2 to $k = n_*$, we obtain

$$\begin{align*} \sum_{k=2}^{n_*} \|Q_k\|_{L^2(\Omega)}&\leqslant c\tau \|u_{\gamma,h}(0)\|_{L^2(\Omega)} \int_\tau^{t_{n_*}}(t_{n_*}-s)^{-\alpha}s^{-1}{\,\mathrm{d}}s\\ &\leqslant c\tau t_{n_*}^{-\alpha-1}\log(n_*+1)\|u_{\gamma,h}(0)\|_{L^2(\Omega)}. \end{align*}$$

As a result, we arrive at

$$\begin{equation*} \|I_2\|_{L^2(\Omega)} \leqslant c\tau t_{n_*}^{-\alpha-1} \log (n_*+1) \|u_{\gamma,h}(0)\|_{L^2(\Omega)}. \end{equation*}$$

Finally, lemmas 4.1, 4.10 and the estimate (4.20) imply that

$$\begin{align*} \|I_3\|_{L^2(\Omega)}&\leqslant c \tau \sum_{k=1}^{n_*} \|A_h(t_{n_*})^2 E_{h,\tau}^{{n_*}-k}(n_*)\|\,\|I-A_h(t_k)A_h(t_{n_*})^{-1}\|\, \| e_k\|_{L^2(\Omega)}\\ &\leqslant c\tau \sum_{k=1}^{n_*}(t_{n_*} - t_k)^{-\alpha}\| e_k\|_{L^2(\Omega)} \\ &\leqslant c \tau^2 \sum_{k=1}^{n_*}(t_{n_*} - t_k)^{-\alpha}t_k^{-1}\log(k+1) \|u_{\gamma,h}(0)\|_{L^2(\Omega)}\\ &\leqslant c\tau \log(n_*+1) t_n^{-\alpha}\|u_{\gamma,h}(0)\|_{L^2(\Omega)}. \end{align*}$$

This completes the proof of the lemma.

Next, we introduce an auxiliary function

$$\begin{equation} \bar\partial_\tau^\alpha U_{h,\gamma}^n +A_h(t_n)U_{h,\gamma}^n = 0 ~~\text{for all}~~1\leqslant n\leqslant N, \quad \text{with}~~ \gamma U_{h,\gamma}^0+U_{h,\gamma}^N = P_hg. \end{equation} \tag{ 4.22 }$$

Then $U_{h,\gamma}^{0}$ can be written as

$$\begin{equation} U_{h,\gamma}^{0} = (\gamma I+F_{h,\tau}^N(N))^{-1}\Big[P_hg-\tau \sum_{k=1}^N F_{h,\tau}^{N-k}(N)(A_h(T)-A_h(t_k))U_{h,\gamma}^{k}\Big]. \end{equation} \tag{ 4.23 }$$

Then the next lemma provides an estimate for $U_{h,\gamma}^{0,\delta}-U_{h,\gamma}^{0}$.

Lemma 4.14. Let $U_{h,\gamma}^{n,\delta}$ and $U_{h,\gamma}^{n}$ be the solution to problems (4.17) and (4.22) respectively. Suppose assumption 2.5 is valid. Then there holds

$$\begin{equation*} \|U_{h,\gamma}^{0,\delta}-U_{h,\gamma}^{0}\|_{L^2(\Omega)}\leqslant c\delta\gamma^{-1}, \end{equation*}$$

where the constant c is independent on γ, h, τ and t.

Proof. Let $e_n = U_{h,\gamma}^{n,\delta} - U_{h,\gamma}^{n}$. Then e_n satisfies the relation that

$$\begin{equation} \bar\partial_\tau^\alpha e_n +A_h(t_n) e_n =0~~\text{for all}~~1\leqslant n\leqslant N, \quad \text{with}~~ \gamma e_0+ e_N = P_h(g_\delta -g). \end{equation} \tag{ 4.24 }$$

Using representations (4.18) and (4.23) we obtain

$$\begin{equation*} \begin{aligned} e_0&=(\gamma I+F_{h,\tau}^{N}(N))^{-1} \Big[P_h(g_\delta-g) - \tau \sum_{k=1}^N E_{h,\tau}^{N-k}(N)(A_h(T)-A_h(t_k)) e_k\Big]. \end{aligned} \end{equation*}$$

Now we apply lemmas 4.10 and 4.11 to obtain

$$\begin{equation*} \begin{aligned} &\| e_0\|_{L^2(\Omega)}\leqslant c\delta\gamma^{-1} + \|F_h^\tau(N;N)^{-1}\tau\sum_{k=1}^NE_h^\tau(N-k;N)(A_h(T)-A_h(t_k)) e_k\|_{L^2(\Omega)}\\ &\leqslant c\delta\gamma^{-1} + c(1+T^{\,\alpha})\sum_{k=1}^N\|\tau A_h(T)E_h^\tau(N-t_k;N)(A_h(T)-A_h(t_k)) e_k\|_{L^2(\Omega)}. \end{aligned} \end{equation*}$$

Then the desired results follows immediately from the a priori estimate in lemma 4.11 and the same argument in theorems 2.2 and 2.4.

Time discretization would give the following fully error estimate.

Lemma 4.15. Let $u_{\gamma,h}(t)$ and $U_{h,\gamma}^{n}$ be the solutions to (4.9) and (4.22) respectively. Suppose assumption 2.5 is valid. Then there holds

$$\begin{align*} \|u_{\gamma,h}(0)- U_{h,\gamma}^{0} \| \leqslant c \tau |\log\tau| (h^2\gamma^{-1}+1) \| u_0 \|_{L^2(\Omega)},\end{align*}$$

where the constant c is independent on γ, h and t.

Proof. Let $\bar U_{h,\gamma}^n$ be the solution to (4.19) and $e_n = \bar U_{h,\gamma}^n- U_{h,\gamma}^n$, which satisfies the following equation

$$\begin{equation} \bar\partial_\tau^\alpha e_n +A_h(t_n)e_n =0~~\text{for all}~~1\leqslant n\leqslant N, \quad \text{with}~~ \gamma e_0 + e_N = \bar U_{h,\gamma}^N - u_{\gamma,h}(T)=:Q. \end{equation} \tag{ 4.25 }$$

Then we apply the representation of fully discrete scheme to derive

$$\begin{equation} e_0 = (\gamma I+F_{h,\tau}^N(N)))^{-1}\Big[Q-\sum_{k=1}^N \tau E_{h,\tau}^{N-k}(N)(A_h(T)-A_h(t_k))e_k\Big]. \end{equation} \tag{ 4.26 }$$

Lemmas 4.10 and 4.12 give that

$$\begin{align*} \|e_0\|_{L^2(\Omega)} &\leqslant \Big\| F_{h,\tau}^N(N)^{-1} \Big[Q-\sum_{k=1}^N \tau E_h^\tau (N-k;N)(A_h(T)-A_h(t_k))e_k\Big]\Big\|_{L^2(\Omega)}\\ &\leqslant c_T\|A_h(T)Q\|_{L^2(\Omega)} + c(1+T^{\,\alpha})\|\sum_{k=1}^N\tau A_h(T) E_{h,\tau}^{N-k}(N) (A_h(T)-A_h(t_k))e_k\|_{L^2(\Omega)}. \end{align*}$$

This combined with lemma 4.13 leads to

$$\begin{align*}\|e_0\|_{L^2(\Omega)}&\leqslant c_T \tau |\log\tau|\|u_{\gamma,h}(0)\|_{L^2(\Omega)}\\&\quad + c(1+T^{\,\alpha})\left\|\sum_{k=1}^N\tau A_h(T) E_{h,\tau}^{N-k}(N) (A_h(T)-A_h(t_k))e_k\right\|_{L^2(\Omega)}. \end{align*}$$

Then by applying the a priori estimate in lemma 4.11 and the same argument in theorems 2.2 and 2.4, we derive

$$\begin{align*}\|e_0\|_{L^2(\Omega)} \leqslant c_T \tau |\log\tau| \|u_{\gamma,h}(0)\|_{L^2(\Omega)}.\end{align*}$$

Finally, the lemmas 3.2 and 4.7 leads to the desired result.

Now we are ready to state the main theorem showing the error of the numerical reconstruction from noisy data. The proof is a direct result of lemma 3.3, 4.7, 4.14 and 4.15.

Theorem 4.1. Let $U_{h,\gamma}^{0,\delta}$ be the numerical reconstructed initial data using the fully discrete scheme (4.17), and u₀ be the exact initial data. Suppose assumption 2.5 is valid. Then if $\| u_0 \|_{\dot H^{q}(\Omega)}\leqslant c$ with $q\in (0,2]$ there holds

$$\begin{equation*} \|U_{h,\gamma}^{0,\delta}-u_0\|_{L^2(\Omega)}\leqslant c\Big(\gamma^\frac q2+ \delta\gamma^{-1} + h^2\gamma^{-1}+\tau {|\log\tau|} (h^2\gamma^{-1}+1) \Big) \end{equation*}$$

Moreover, for $u_0\in L^2(\Omega)$, there holds

$$\begin{equation*} \| U_{h,\gamma}^{0,\delta} - u(0)\|_{L^2(\Omega)} \rightarrow 0\quad \text{as}~~\delta,\gamma,h,\tau\rightarrow0, ~~\frac\delta\gamma\rightarrow0~~\text{and}~~\frac{h^2}{\gamma} \rightarrow0. \end{equation*}$$

The a priori error estimate in theorem 4.1 gives a useful guideline to choose the regularization parameter γ and the discretization parameters h and τ according to the noise level δ. In particular, if $u_0 \in \dot H^{q}(\Omega)$, by choosing

$$\begin{equation*} \gamma \sim \delta^{\frac{2}{q+2}},~~h\sim \delta^\frac 12~~\text{and}~~\tau |\log \tau| \sim \delta^{\frac{q}{q+2}}, \end{equation*}$$

we obtain the optimal approximation error

$$\begin{equation*} \|U_{h,\gamma}^{0,\delta} - u(0)\|_{L^2(\Omega)}\leqslant c\delta^{\frac{q}{q+2}}. \end{equation*}$$

Now we test several two dimensional examples with $\Omega = (0,1)^2$ in order to illustrate our theoretical results. Throughout the section, we apply the standard Galerkin piecewise linear FEM with uniform mesh size $h = 1/(M+1)$ for the space discretization, and the CQ-BE method with uniform mesh size $\tau = T/N$ for time discretization. We solve the direct problem to obtain the exact observation data by using fine meshes, i.e. $h = 1/100$, $\tau = T/500$. Then we compute the noisy observational data by

$$\begin{equation*} g_\delta = u(T)+ \varepsilon\delta\sup_{x\in\Omega} u(x,T) \end{equation*}$$

where ε is generated from standard Gaussian distribution and δ denotes the related noisy level.

We begin with the following time-dependent diffusion coefficient:

$$\begin{equation*} a_1(x,y,t) = \begin{pmatrix} y\sin((1+t)^{0.5})+2 & -0.1 \\ -0.1 & \sin(\pi x)(t+1.2)^{-0.8}+2 \end{pmatrix}, \end{equation*}$$

satisfying conditions (1.2), (1.3) and assumption 2.3. We solve the linear system (4.17) by using the conjugate gradient method.

Example 5.1 (Smooth initial data). We begin with a smooth initial data:

$$\begin{align*}u_0 = \sin(2\pi x)\sin(2\pi y) \in \dot H^{2}(\Omega).\end{align*}$$

According to theorem 4.1, we compute $U_{h,\gamma}^{0,\delta}$ with $\gamma \sim \sqrt{\delta}$ and $h,\tau \sim \sqrt{\delta}$, and expect a convergence of order $O(\sqrt{\delta})$. Numerical results presented in figure 1 fully support the theoretical result. On the other hand, our numerical results indicate that the recovery is stable for all T, might be neither very large nor very small. This interesting phenomenon warrants further investigation in the future. In figure 2, we present profiles of solutions and errors with different noise level.

Example 5.2 (Nonsmooth initial data). In this example we consider the following nonsmooth initial condition

$$\begin{equation*} u_0 = \begin{cases} 1,~~\text{if}~0.5 \leqslant x \leqslant 1,\\ 0,~~\text{otherwise}. \end{cases} \end{equation*}$$

Note that $u_0 \in \dot H^{\frac12-\varepsilon}(\Omega)$ for any $\varepsilon\in (0,\frac 12)$. Then theorem 4.1 indicate that the optimal convergence rate is almost $O(\delta^{0.2})$ provided that $\gamma = O(\delta^{0.8})$, $h = O(\sqrt{\delta})$ and $\tau = O(\delta^{0.2})$. This is fully supported by the numerical results presented in figure 3. In figure 4 we plot the profiles of solutions and errors, which also confirm that the numerical recovery is reliable.

Example 5.3 (Diffusion coefficient violating assumption 2.3). We also test the following diffusion coefficient

$$\begin{equation*} a_2(x,y,t) = \begin{pmatrix} e^{-x}\cos(t)+2 & (1.5-(t+1)^{-0.2})/10\\ (1.5-(t+1)^{-0.2})/10 & \cos(\pi y)\sin(t)+2 \end{pmatrix}. \end{equation*}$$

Note that a₂ satisfies conditions (1.2) and (1.3), but assumption 2.3 is not fulfilled.

Numerical experiments show that the numerical reconstruction via the fully discrete scheme (4.17) still converges under proper parameter choices. For example, we test the smooth initial data $u_0 = \sin(2\pi x)\sin(2\pi y)$ and large terminal time T = 10. We choose $\gamma,h,\tau\sim\sqrt{\delta}$, and observe a convergence rate around $O(\sqrt{\delta})$, see figure 5. We will continue to consider the general case in our future studies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The research of Z Zhou is partially supported by Hong Kong Research Grants Council (Project No. 15304420) and an internal grant of Hong Kong Polytechnic University (Project ID: P0031041, Work Programme: ZZKS).

No new data were created or analysed in this study.

For p = 0, conditions (1.2), (1.3) and assumption 2.3 imply

$$\begin{align*} \|(A(t)-A(s))v\|_{L^2(\Omega)} &\leqslant c\big(\|\nabla(a(t)-a(s))\|_{L^\infty(\Omega)}+\|a(t)-a(s)\|_{L^\infty(\Omega)}\big)\|v\|_{\dot H^{2}(\Omega)}\\ &\leqslant c \min(1,\min(t,s)^{-\kappa}|t-s|)\|v\|_{\dot H^{2}(\Omega)}. \end{align*}$$

For $p = -2$, from using the duality argument, we have

$$\begin{align*} \|(A(t)-A(s))v\|_{\dot H^{-2}(\Omega)} &= \sup_{\varphi\in \dot H^{2}(\Omega)} \frac{\langle(A(t)-A(s))v,\varphi\rangle}{\|\varphi\|_{\dot H^{2}(\Omega)}} = \sup_{\varphi\in \dot H^{2}(\Omega)} \frac{(v,(A(t)-A(s))\varphi)}{\|\varphi\|_{\dot H^{2}(\Omega)}} \\ & \leqslant \sup_{\varphi\in \dot H^{2}(\Omega)} \frac{\| v \|_{L^2(\Omega)} \| (A(t)-A(s))\varphi \|_{L^2(\Omega)}}{\|\varphi\|_{\dot H^{2}(\Omega)}}\\ &\leqslant c \min(1,\min(t,s)^{-\kappa}|t-s|) \sup_{\varphi\in \dot H^{2}(\Omega)}\| v \|_{L^2(\Omega)} \end{align*}$$

This completes the proof of the lemma.

Recalling the fact that [15, lemma 3.3]

$$\begin{align*} U_h^n\bar\partial_\tau^\alpha U_h^n \geqslant \frac 12 \bar\partial_\tau^\alpha |U_h^n|^2. \end{align*}$$

Therefore like lemma 2.3 we define an operator $\underline{A_h} = -c_0 \Delta_h$. Condition (1.2) gives that the operator $A_h(t) - \underline{A_h}$ is selfadjoint and positive semidefinite for all $n\geqslant 1$. Rewrite the equation (4.14) as

$$\begin{align*} \bar\partial_\tau^\alpha (U_h^n-U_h^0) + \underline{A_h} U_n = (\underline{A_h}-A_h(t)) U_h^n \quad \text{for all}~~1\leqslant n\leqslant N. \end{align*}$$

Taking inner product with $U_h^n$ on the above equation and by definition of $-\Delta_h$ and $A_h(t)$, we obtain

$$\begin{align*} (\bar\partial_\tau^\alpha (U_h^n-U_h^0), U_h^n) + c_0 \| \nabla U_h^n \|_{L^2(\Omega)}^2 = \big((c_0-a(\cdot,t)) \nabla U_h^n, \nabla U_h^n\big)\leqslant 0 \quad \text{for all}~~1\leqslant n\leqslant N. \end{align*}$$

Using the above inequality and Poincaré inequality we arrive at

$$\begin{align*} &\bar\partial_\tau^\alpha(\|U_h^n\|_{L^2(\Omega)}-\|U_h^0\|_{L^2(\Omega)})+c\|U_h^n\|_{L^2(\Omega)}\\ &\leqslant \bar\partial_\tau^\alpha\Big[( \|U_h^n\|_{L^(\Omega)}-\|U_h^0\|_{L^2(\Omega)})(1+\|U_h^0\|_{L^2(\Omega)}/\|U_h^n\|_{L^2(\Omega)}) \Big] + c \|U_h^n\|_{L^2(\Omega)} \\ &\leqslant 0 \quad \text{for all}~~n\geqslant 1, \end{align*}$$

for a constant c uniform in t_n . Then the comparison principle for discrete fractional ODEs [19] leads to

$$\begin{align*} \|U_h^n\|_{L^2(\Omega)} \leqslant F_\tau^n(c)\| U_h^0 \|_{L^2(\Omega)} \leqslant \frac{c}{1+ct_n^\alpha} \| U_h^0\|_{L^2(\Omega)} \end{align*}$$

where the definition of $F_\tau^n(c)$ can be found in [41, 42]. This immediately leads to the desired result.

Next by solution representation (4.15) we have

$$\begin{align*} &\|A_h(t_{n_*})U_h^n\|_{L^2(\Omega)}\\ & \leqslant \|A_h(t_{n_*})F_{h,\tau}^n (n_*)U_h^0\|_{L^2(\Omega)} + \tau\| \sum_{k=1}^n A_h(t_{n_*})E_{h,\tau}^{n-k} (n_*)(I-A_h(t_k)A_h(t_{n_*})^{-1})\\ & \quad A_h(t_{n_*})U_h^k\|_{L^2(\Omega)}\\ &\leqslant c t_n^{-\alpha} \|U_h^0\|_{L^2(\Omega)} + \sum_{k=1}^n \|\tau A_h(t_{n_*})E_{h,\tau}^{n-k} (n_*)(I-A_h(t_k)A_h(t_{n_*})^{-1})\|\,\\ & \quad \|A_h(t_{n_*})U_h^k\|_{L^2(\Omega)}, \end{align*}$$

lemmas 4.1 and 4.10 show that

$$\begin{align*} \|A_h(t_{n_*})U_h^n\|_{L^2(\Omega)}\leqslant c t_n^{-\alpha} \|U_h^0\|_{L^2(\Omega)} + \sum_{k=1}^n c\tau\|A_h(t_{n_*})U_h^k\|_{L^2(\Omega)}, \end{align*}$$

the discrete version of Gronwall's inequality [33, lemma 10.5] gives that

$$\begin{equation*} \|A_h(t_{n_*})U_h^n\|_{L^2(\Omega)}\leqslant c \exp(ct_n) t_n^{-\alpha} \|U_h^0\|_{L^2(\Omega)} \end{equation*}$$

here c is uniform in n, τ and t_n .

Meanwhile, in the other hand $\| I - A(t_*)^{-1} A(s) \| \leqslant c |t_*-s|^{\beta}$ for any $\beta \in [0,1]$. Then if $\beta = (1+\epsilon)\alpha$ with $\epsilon\in(0,1/\alpha-1)$ we can derive that

$$\begin{align*} &\|A_h(t_{n_*})U_h^n\|_{L^2(\Omega)} \leqslant \|A_h(t_{n_*})F_{h,\tau}^n (n_*)U_h^0\|_{L^2(\Omega)} \\ &+ \tau\| \sum_{k=1}^n A_h(t_{n_*})E_{h,\tau}^{n-k} (n_*)(I-A_h(t_k)A_h(t_{n_*})^{-1})A_h(t_{n_*})U_h^k\|_{L^2(\Omega)}\\ &\leqslant c t_n^{-\alpha} \|U_h^0\|_{L^2(\Omega)} + \sum_{k=1}^n \|\tau A_h^2(t_{n_*})E_{h,\tau}^{n-k} (n_*)(I-A_h(t_{n_*})^{-1}A_h(t_k))\|\,\|U_h^k\|_{L^2(\Omega)}\\ &\leqslant c t_n^{-\alpha} \|U_h^0\|_{L^2(\Omega)} +c\tau \sum_{k=1}^n (t_{n_*}-t_k)^{-1+\varepsilon \alpha} t_k^{-\alpha} \|U_h^0\|_{L^2(\Omega)}\\ &\leqslant c t_n^{-\alpha} \|U_h^0\|_{L^2(\Omega)} +c\int_0^{t_n}(t_{n_*}-s)^{-1+\epsilon \alpha}s^{-\alpha}{\,\mathrm{d}}s\\ & \quad \|U_h^0\|_{L^2(\Omega)}\leqslant ct_n^{-(1-\epsilon)\alpha} \|U_h^0\|_{L^2(\Omega)}. \end{align*}$$