Numerical reconstruction of the spatial component in the source term of a time-fractional diffusion equation

TL;DR

This paper develops a numerical method using Tikhonov regularization and iterative thresholding to reconstruct the spatial component of a source term in a time-fractional diffusion equation, demonstrating stability, convergence, and efficiency.

Contribution

It introduces a stabilized nonlinear minimization approach with a novel convergence proof for the discrete scheme in reconstructing the source term.

Findings

The method is stable and convergent under the proposed scheme.

Numerical experiments confirm the efficiency and accuracy of the iterative thresholding algorithm.

The approach effectively reconstructs the spatial component in fractional diffusion problems.

Abstract

In this article, we are concerned with the analysis on the numerical reconstruction of the spatial component in the source term of a time-fractional diffusion equation. This ill-posed problem is solved through a stabilized nonlinear minimization system by an appropriately selected Tikhonov regularization. The existence and the stability of the optimization system are demonstrated. The nonlinear optimization problem is approximated by a fully discrete scheme, whose convergence is established under a novel result verified in this study that the $H^1$-norm of the solution to the discrete forward system is uniformly bounded. The iterative thresholding algorithm is proposed to solve the discrete minimization, and several numerical experiments are presented to show the efficiency and the accuracy of the algorithm.

Tables

Table 1. Table 1 : Numerical performance in Example 5.2 under various combinations of the observation subdomain ω 𝜔 \omega and the noise level δ 𝛿 \delta .

$\delta$	$\omega$	$\mathrm{err}$	$K$
$1\%$	$\Omega \setminus [1/10,9/10]$	$3.75\%$	$13$
$1\%$	$\Omega \setminus [1/20,19/20]$	$3.38\%$	$18$
$1\%$	$\Omega \setminus [1/40,39/40]$	$4.17\%$	$27$
$0.5\%$	$\Omega \setminus [1/20,19/20]$	$2.96\%$	$16$
$2\%$	$\Omega \setminus [1/20,19/20]$	$5.53\%$	$20$
$4\%$	$\Omega \setminus [1/20,19/20]$	$11.17\%$	$24$

Table 2. Table 2 : Numerical performance in Example 5.4 under various combinations of the observation subdomain ω 𝜔 \omega and the noise level δ 𝛿 \delta .

$\delta$	$\omega$	$\mathrm{err}$	$K$
$1\%$	$\Omega \setminus {[1/10,9/10]}^2$	$2.63\%$	$13$
$1\%$	$\Omega \setminus {[1/20,19/20]}^2$	$4.94\%$	$27$
$1\%$	$\Omega \setminus [0,9/10]\times [1/10,9/10]$	$4.36\%$	$9$
$0.5\%$	$\Omega \setminus {[1/10,9/10]}^2$	$2.47\%$	$24$
$1\%$	$\Omega \setminus {[1/10,9/10]}^2$	$3.61\%$	$15$
$2\%$	$\Omega \setminus {[1/10,9/10]}^2$	$5.06\%$	$7$
$4\%$	$\Omega \setminus {[1/10,9/10]}^2$	$6.58\%$	$5$

Equations

$${\partial }_t^{\alpha }g(t):=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{g^{\prime }(s)}{(t-s{)}^{\alpha }}\mathrm{d}s,$$

$$\{\begin{array}{ll}({\partial }_t^{\alpha }u-△+1)u(x,t)=f(x)\mu (t),& (x,t)\in Q:=\Omega \times (0,T),\\ u(x,0)=0,& x\in \Omega ,\\ {\partial }_{\nu }u(x,t)=0,& (x,t)\in \partial \Omega \times (0,T),\end{array}$$

$${\partial }_t^{\alpha }u(x,t)-\sum _{i,j=1}^d\frac{\partial }{\partial x_i}\left(a_{ij}(x)\frac{\partial u(x,t)}{\partial x_j}\right)+c(x)u(x,t)=f(x)\mu (t),$$

$$z_n\rightharpoonup z\text{mbox}inX\text{mbox}asn\to \infty .$$

$${}_0{H}^{\alpha }(0,T):=\{\begin{array}{rlrl}& H^{\alpha }(0,T),& & 0<\alpha <\frac{1}{2},\\ & \left\{g\in H^{\frac{1}{2}}(0,T);{\int }_0^T\frac{|g(t){|}^2}{t}\mathrm{d}t<\infty \right\},& & \alpha =\frac{1}{2},\\ & \{g\in H^{\alpha }(0,T);g(0)=0\},& & \frac{1}{2}<\alpha <1.\end{array}$$

$$J_{0+}^{\alpha }g(t):=\frac{1}{\Gamma (\alpha )}{\int }_0^t\frac{g(\tau )}{(t-\tau {)}^{1-\alpha }}\mathrm{d}\tau ,J_{T-}^{\alpha }g(t)=\frac{1}{\Gamma (\alpha )}{\int }_t^T\frac{g(\tau )}{(\tau -t{)}^{1-\alpha }}\mathrm{d}\tau$$

$${\int }_0^T(J_{0+}^{\alpha }{g}_1)g_2\mathrm{d}t={\int }_0^T{g}_1{J}_{T-}^{\alpha }{g}_2\mathrm{d}t,\forall g_1,g_2\in L^2(0,T).$$

$${\int }_0^T({\partial }_t^{\alpha }g)h\mathrm{d}t=-g(0)J_{T-}^{1-\alpha }h(0)-{\int }_0^{T}g(J_{T-}^{1-\alpha }h{)}^{\prime }\mathrm{d}t.$$

$${\int }_0^T({\partial }_t^{\alpha }{g}_n)h\mathrm{d}t$$

$$=-g_n(0)J_{T-}^{1-\alpha }h(0)-{\int }_0^T{g}_n(J_{T-}^{1-\alpha }h{)}^{\prime }\mathrm{d}t,$$

$$\Vert u(f){\Vert }_{H^{\alpha }(0,T;L^2(\Omega ))}+\Vert u(f){\Vert }_{C([0,T];L^2(\Omega ))}+\Vert u(f){\Vert }_{L^2(0,T;H^2(\Omega ))}\le C\Vert f{\Vert }_{L^2(\Omega )}.$$

$${\Vert u^{\delta }-u(f^{\ast })\Vert }_{L^2(\omega \times (0,T))}\le \delta ,$$

$$\underset{f\in L^2(\Omega )}{\min }J(f),J(f):={\Vert u(f)-u^{\delta }\Vert }_{L^2(\omega \times (0,T))}^2+\beta \Vert f{\Vert }_{L^2(\Omega )}^2,$$

$${\int }_0^T{\int }_{\Omega }({\partial }_t^{\alpha }u(f)\phi +\nabla u(f)\cdot \nabla \phi +u(f)\phi )\mathrm{d}x\mathrm{d}t={\int }_0^T{\int }_{\Omega }f\mu \phi \mathrm{d}x\mathrm{d}t.$$

$$\underset{n\to \infty }{\lim }J(f^n)=\underset{f\in L^2(\Omega )}{\inf }J(f).$$

$$f^n\rightharpoonup f^{\ast }\text{mbox}in{L}^2(\Omega )\text{mbox}asn\to \infty .$$

$$u(f^n)\rightharpoonup u^{\ast }\text{mbox}in{}_0{H}^{\alpha }(0,T;L^2(\Omega ))\cap L^2(0,T;H^2(\Omega )).$$

$${\int }_0^T{\int }_{\Omega }({\partial }_t^{\alpha }u(f^n)\phi +\nabla u(f^n)\cdot \nabla \phi +u(f^n)\phi )\mathrm{d}x\mathrm{d}t={\int }_0^T{\int }_{\Omega }{f}^n\mu \phi \mathrm{d}x\mathrm{d}t$$

$${\partial }_t^{\alpha }u(f^n)\rightharpoonup {\partial }_t^{\alpha }{u}^{\ast },\nabla u(f^n)\rightharpoonup \nabla u^{\ast }\text{mbox}in{L}^2(Q),$$

$${\int }_0^T{\int }_{\Omega }({\partial }_t^{\alpha }{u}^{\ast }\phi +\nabla u^{\ast }\cdot \nabla \phi +u^{\ast }\phi )\mathrm{d}x\mathrm{d}t={\int }_0^T{\int }_{\Omega }{f}^{\ast }\mu \phi \mathrm{d}x\mathrm{d}t$$

$${\int }_{\Omega }{\int }_0^T{\partial }_t^{\alpha }u(f^n)\psi v\mathrm{d}t\mathrm{d}x=-{\int }_{\Omega }u(f^n)(\cdot ,0)J_{T^{-}}^{1-\alpha }\psi (0)v\mathrm{d}x-{\int }_{\Omega }{\int }_0^{T}u(f^n)\frac{\mathrm{d}}{\mathrm{d}t}{J}_{T^{-}}^{1-\alpha }\psi v\mathrm{d}t\mathrm{d}x.$$

$${\int }_{\Omega }{\int }_0^T{\partial }_t^{\alpha }{u}^{\ast }\psi v\mathrm{d}t\mathrm{d}x=-{\int }_{\Omega }{\int }_0^T{u}^{\ast }\frac{\mathrm{d}}{\mathrm{d}t}{J}_{T^{-}}^{1-\alpha }\psi v\mathrm{d}t\mathrm{d}x.$$

$${\int }_{\Omega }{\int }_0^T{\partial }_t^{\alpha }{u}^{\ast }\psi v\mathrm{d}t\mathrm{d}x=-{\int }_{\Omega }{u}^{\ast }(\cdot ,0)J_{T^{-}}^{1-\alpha }\psi (0)v\mathrm{d}x-{\int }_{\Omega }{\int }_0^T{u}^{\ast }\frac{\mathrm{d}}{\mathrm{d}t}{J}_{T^{-}}^{1-\alpha }\psi v\mathrm{d}t\mathrm{d}x$$

$$J(f^{\ast })$$

$$\le \underset{n\to \infty }{\mathrm{liminf}}{\Vert u(f^n)-u^{\delta }\Vert }_{L^2(\omega \times (0,T))}^2+\beta \underset{n\to \infty }{\mathrm{liminf}}\Vert f^n{\Vert }_{L^2(\Omega )}^2$$

$$\le \underset{n\to \infty }{\mathrm{liminf}}J(f^n)=\underset{f\in L^2(\Omega )}{\inf }J(f),$$

$$u_{\ell }^{\delta }\to u^{\delta }\text{mbox}in{L}^2(0,T;L^2(\omega ))\text{mbox}as\ell \to \infty ,$$

$$\underset{f\in L^2(\Omega )}{\min }{J}_{\ell }(f),J_{\ell }(f):={\Vert u(f)-u_{\ell }^{\delta }\Vert }_{L^2(\omega \times (0,T))}^2+\beta \Vert f{\Vert }_{L^2(\Omega )}^2,\ell =1,2,\dots .$$

$$J_{\ell }(f^{\ell })\le J_{\ell }(f),\forall f\in L^2(\Omega ),$$

$$f^{\ell }\rightharpoonup f^{\ast }\text{mbox}in{L}^2(\Omega )\text{mbox}as\ell \to \infty .$$

$$u(f^{\ell })\rightharpoonup u(f^{\ast })\text{mbox}in{}_0{H}^{\alpha }(0,T;L^2(\Omega ))\cap L^2(0,T;H^2(\Omega ))\text{mbox}as\ell \to \infty ,$$

$$u(f^{\ell })-u_{\ell }^{\delta }\rightharpoonup u(f^{\ast })-u^{\delta }\text{mbox}in{L}^2(0,T;L^2(\omega ))\text{mbox}as\ell \to \infty .$$

Full text

Numerical Reconstruction of the Spatial Component in the Source Term of a Time-Fractional Diffusion Equation

Daijun Jiang111School of Mathematics and Statistics & Hubei Key Laboratory of Mathematical Sciences, Central China Normal University, Wuhan 430079, China. E-mail: jiangdaijun@mail.ccnu.edu.cn

  

Yikan Liu222Corresponding author. Research Center of Mathematics for Social Creativity, Research Institute for Electronic Science, Hokkaido University, N12W7, Kita-Ward, Sapporo 060-0812, Japan. E-mail: ykliu@es.hokudai.ac.jp

  

Dongling Wang333Department of Mathematics and Center for Nonlinear Studies, Northwest University, Xi’an, Shaanxi 710075, China. E-mail: wdymath@nwu.edu.cn

ABSTRACT

In this article, we are concerned with the analysis on the numerical reconstruction of the spatial component in the source term of a time-fractional diffusion equation. This ill-posed problem is solved through a stabilized nonlinear minimization system by an appropriately selected Tikhonov regularization. The existence and the stability of the optimization system are demonstrated. The nonlinear optimization problem is approximated by a fully discrete scheme, whose convergence is established under a novel result verified in this study that the $H^{1}$-norm of the solution to the discrete forward system is uniformly bounded. The iterative thresholding algorithm is proposed to solve the discrete minimization, and several numerical experiments are presented to show the efficiency and the accuracy of the algorithm.

Keywords  Time-fractional diffusion equation, Inverse source problem, Finite element method, Iterative thresholding algorithm

Mathematics subject classification (2010)  35R11, 65M32, 41A35

1 Introduction

Let $T>0$ and $\Omega\subset{\mathbb{R}}^{d}$ ($d=1,2,3$) be a bounded convex polygonal domain. For $0<\alpha<1$, by $\partial_{t}^{\alpha}$ we denote the Caputo derivative defined as (see, e.g., [10, p.​ 91])

[TABLE]

where $\Gamma(\,\cdot\,)$ is the usual Gamma function. In this paper, we consider the following initial-boundary value problem for a time-fractional diffusion equation (TFDE) with the homogeneous Neumann boundary condition

[TABLE]

where $\partial_{\nu}u:=\nabla u\cdot\nu$ and $\nu=(\nu_{1},\ldots,\nu_{d})$ denotes the unit outward normal derivative on the boundary $\partial\Omega$. The source term in (1.1) takes the form of separated variables, in which the spatial component $f(x)$ models the spatial distribution e.g.​ of the contaminant source, and the temporal component $\mu(t)$ describes the time evolution pattern.

The governing equation in (1.1) is a typical representative of a wide range of TFDEs, which were proposed as a powerful candidate for describing anomalous diffusion phenomena in heterogenous media. In the past decades, TFDEs have been investigated thoroughly from both theoretical and numerical aspects. It reveals from [19, 16, 5] and other literature that TFDEs resemble their integer counterpart (i.e., $\alpha=1$) in such qualitative aspects like analyticity and maximum principle, but show certain difference in the senses of limited smoothing effect, short-/long-time asymptotic behavior and weak unique continuation. For the numerical simulation of TFDEs, we refer e.g.​ to [14, 11, 9, 6, 21].

Simultaneously, various kinds of inverse problems for TFDEs have also gathered consistent popularity within the last decade. Owing to the limited smoothing property of the forward problem, usually the ill-posedness of inverse problems for TFDEs is less severe than that for the case of $\alpha=1$, which can be witnessed from the backward problems. Nevertheless, due to the lack of techniques for analysis, for some inverse problems, one should relate TFDEs with other types of equations, so that one can indirectly obtain uniqueness and sometimes stability. For the inverse source problems on determining the spatial component, it turns out that the majority of existing literature dealt with the final observation data $u(\,\cdot\,,T)$ except for [29]. We refer to Sakamoto and Yamamoto [20] for the generic well-posedness, and [24, 26] etc.​ for the numerical reconstruction by various regularization methods. Recently, by a newly established unique continuation property for TFDEs, Jiang et al. [5] proved the uniqueness for the same problem by the partial interior observation and developed an iterative thresholding algorithm. On reconstructing the temporal component in the source term, we refer e.g.​ to [18, 25, 17]. For a complete bibliography on inverse problems for TFDEs, see the topical review articles [8, 15, 13, 12].

In the same direction of [5], in this paper we mainly focus on the numerical aspect of the following inverse source problem.

Problem 1.1**.**

Let $\omega\subset\Omega$ be a nonempty open subdomain and $u(f)$ be the solution to (1.1). Provided that $\mu$ is known on $[0,T]$, determine the spatial component $f$ in $\Omega$ by the partial interior observation of $u(f)$ in $\omega\times(0,T)$.

In fact, here we can also consider a more general formulation such as

[TABLE]

where we assume certain regularity of $a_{ij}$ and $c$, and $(a_{ij})_{d\times d}$ is a strictly positive-definite matrix on $\overline{\Omega}\,$. However, we restrict ourselves to the formulation (1.1) not only for its simplicity in the numerical simulation, but also due to the belief that the underlying ill-posedness are essentially the same.

On the other hand, we notice that the iteration method proposed in [5] lacks detailed analysis especially on the convergence issue, although the numerical performance was discussed. From the practical viewpoint, it is obligatory to reformulate Problem 1.1 in a discrete setting and investigate whether the resulting system inherits the corresponding properties of the continuous one. Moreover, we should verify the convergence of the discretized solution to the continuous one in some sense. This motivates the central topic of the present article, which turns out to be a very important and necessary supplementation to the investigation of Problem 1.1.

To fully discretize system (1.1), we shall employ the standard Galarkin method with piecewise linear finite element in space and the L1 scheme in time. This method is one of the most popular and successful numerical scheme for discretizing the subdiffusion problem, and it has been analyzed from various aspects. The optimal error estimate in $L^{2}(\Omega)$ norm with respect to the regularity of the problem solution is established in [6] for uniform time step size. In order to deal with the weak singularity of time-fractional diffusion equations, the L1-type scheme on graded time meshes is employed in the discrete scheme, and the optimal error estimate in $L^{\infty}(\Omega)$ norm is established for $d=1$ in [21] and $d=2,3$ in [9]. For the convergence analysis, in this paper we further establish the uniform $H^{1}{(\Omega)}$ error estimate on $[0,T]$ by using the Poincaré-Friedrichs inequality and the discrete fractional Grönwall’s inequality for the L1 method, which turns out to be novel to the best of our knowledge.

The remainder of this paper is organized as follows. Recalling key ingredients in theory, in Section 2 we interpret Problem 1.1 as an optimization problem with Tikhonov regularization, and demonstrate its regularizing effects. In Section 3, we discretize the optimization problem by a fully discrete finite element approximation and study its basic properties. Next, Section 4 is devoted to the convergence analysis of the solution to the discrete optimization problem. Finally, we propose the iterative thresholding algorithm to solve the discrete problem in Section 5, and implement several numerical examples to show the efficiency and accuracy of the algorithm.

2 Preliminary and Tikhonov regularization

In this section, we shall make general preparations and formulate Problem 1.1 stated in Section 1 as a stabilized minimization system, and establish the unique existence of the solution as well as the stability of the minimization formulation.

Let $L^{2}(\Omega)$ be a usual $L^{2}$-space equipped with the inner product $(\,\cdot\,,\,\cdot\,)$, and $H^{\alpha}(0,T)$, $H^{\gamma}(\Omega)$ ($\gamma\in{\mathbb{R}}$) etc.​ denote Sobolev spaces (e.g., Adams [1]). Throughout this paper, $C>0$ stands for generic constants which may change line by line. In a Banach space $X$, we denote the weak convergence of a sequence $\{z_{n}\}$ to $z$ by

[TABLE]

We first revisit some basic facts concerning the initial-boundary value problem (1.1). For the solution regularity, we define the fractional Sobolev spaces ${{}_{0}}H^{\alpha}(0,T)$ as

[TABLE]

We also define the Riemann-Liouville integral operator $J_{0+}^{\alpha}$ and the backward Riemann-Liouville integral operator $J^{\alpha}_{T-}$ respectively as

[TABLE]

for $g\in L^{2}(0,T)$. We know that the Caputo derivative can be rewritten as $\partial_{t}^{\alpha}g=J_{0+}^{1-\alpha}\frac{\mathrm{d}}{\mathrm{d}t}g$, and the following relation holds (see [5, Lemma 4.1])

[TABLE]

For later use, we show a technical lemma on the basis of (2.1)

Lemma 2.1**.**

Let $g\in{{}_{0}}H^{\alpha}(0,T)\cap C[0,T]$ and $h\in C^{1}[0,T]$ with $J_{T-}^{1-\alpha}h(T)=0$. Then

[TABLE]

Proof.

Pick a sequence $\{g_{n}\}\subset C^{\infty}[0,T]$ such that $g_{n}\to g$ in ${{}_{0}}H^{\alpha}(0,T)\cap C[0,T]$. Then for each $n=1,2,\ldots$, we utilize (2.1) and calculate

[TABLE]

where the last equality follows from the assumption $J_{T-}^{1-\alpha}h(T)=0$. Then the proof is completed by passing $n\to\infty$ on both sides of the above identity. ∎

On the basis of [5, Lemma 2.4] and [19, Theorem 2.2(i)], we summarize the well-posedness of the forward problem (1.1) as follows.

Lemma 2.2**.**

Let $f\in L^{2}(\Omega)$ and $\mu\in L^{\infty}(0,T)$. Then the initial-boundary value problem (1.1) admits a unique solution $u(f)\in{{}_{0}}H^{\alpha}(0,T;L^{2}(\Omega))\cap C([0,T];L^{2}(\Omega))\cap L^{2}(0,T;H^{2}(\Omega))$. Moreover, there exists a constant $C>0$ depending on $\Omega,T,\alpha$ and $\mu$ such that

[TABLE]

In [19, Theorem 2.2(i)], it was proved that the solution $u\in C([0,T];H^{-2\gamma}(\Omega))$ for any $\gamma>\frac{d}{4}-1$ provided that the source term belongs to $L^{\infty}(0,T;L^{2}(\Omega))$. Since we restrict the spatial dimensions to $d=1,2,3$ throughout this paper, one can choose negative $\gamma$ and thus $u(f)$ makes sense in $C([0,T];L^{2}(\Omega))$ under the assumptions of Lemma 2.2.

For the theoretical aspect of Problem 1.1, we recall the uniqueness result stated in [5, Theorem 2.6].

Lemma 2.3**.**

Let $\omega\subset\Omega$ be an arbitrarily chosen open subdomain and $u(f)$ be the solution to (1.1). Assume that $f\in L^{2}(\Omega)$ and $\mu\in C^{1}[0,T]$ with $\mu(0)\neq 0$. Then $u=0$ in $\omega\times(0,T)$ implies $f=0$ in $\Omega$.

Suppose that we are given the noisy observation data $u^{\delta}\in L^{2}(\omega\times(0,T))$ in practice. Usually, $u^{\delta}$ satisfies

[TABLE]

where $f^{*}\in L^{2}(\Omega)$ and $\delta>0$ stand for the true source term and the noise level, respectively. To deal with the ill-posedness of Problem 1.1, we still adopt a classical Tikhonov regularization methodology as that in [5] to consider the following optimization problem with Tikhonov regularization

[TABLE]

where $\beta>0$ is the regularization parameter, and $u=u(f)$ satisfies the initial condition $u(x,0)=0$ and the variational system associated with the equation (1.1): for any $\varphi\in H^{\alpha}(0,T;L^{2}(\Omega))\cap L^{2}(0,T;H^{1}(\Omega))$,

[TABLE]

Regarding the minimizer of (2.2), we first show the following unique existence result.

Theorem 2.1**.**

For any $u^{\delta}\in L^{2}(\omega\times(0,T))$, there exists a unique minimizer $f^{*}\in L^{2}(\Omega)$ to the optimization problem (2.2).

Proof.

Since $J(f)$ is nonnegative, we know that $\inf_{f\in L^{2}(\Omega)}J(f)$ is finite. Thus there exists a minimizing sequence $\{f^{n}\}\subset L^{2}(\Omega)$ such that

[TABLE]

Then by the definition of $J(f^{n})$, it is obvious that $\{f^{n}\}$ is uniformly bounded in $L^{2}(\Omega)$. Therefore, there exist $f^{*}\in L^{2}(\Omega)$ and a subsequence of $\{f^{n}\}$, still denoted by $\{f^{n}\}$, such that

[TABLE]

We shall prove that $f^{*}$ is indeed the unique minimizer of (2.2).

Since each $f^{n}$ corresponds with a solution $u(f^{n})$ to (1.1) with $f=f^{n}$, it follows immediately from Lemma 2.2 that the sequence $\{u(f^{n})\}$ is also uniformly bounded in ${{}_{0}}H^{\alpha}(0,T;L^{2}(\Omega))\cap L^{2}(0,T;H^{2}(\Omega))$. This, together with the compact embedding $H^{2}(\Omega)\subset\subset L^{2}(\Omega)$, indicates the existence of some $u^{*}\in{{}_{0}}H^{\alpha}(0,T;L^{2}(\Omega)\cap L^{2}(0,T;H^{2}(\Omega))$ and a subsequence of $\{u(f^{n})\}$, again still denoted by $\{u(f^{n})\}$, such that

[TABLE]

We claim $u^{*}=u(f^{*})$. Actually, we utilize the fact that

[TABLE]

for all $\varphi\in L^{2}(0,T;H^{1}(\Omega))$. Since (2.4) implies

[TABLE]

we pass $n\to\infty$ in (2.5) to obtain

[TABLE]

for all $\varphi\in L^{2}(0,T;H^{1}(\Omega))$. Next, we shall prove $u^{*}(\,\cdot\,,0)=0$, which with (2.6) and the definition (2.3) of weak solutions implies $u^{*}=u(f^{*})$.

Indeed, it follows from Lemma 2.2 that $u(f^{n})\in{{}_{0}}H^{\alpha}(0,T;L^{2}(\Omega))\cap C([0,T];L^{2}(\Omega))$, which satisfies the assumption of Lemma 2.1. Then by taking $\psi\in C^{1}[0,T]$ with $J^{1-\alpha}_{T^{-}}\psi(T)=0$ and $v\in L^{2}(\Omega)$ arbitrarily, we employ Lemma 2.1 to deduce

[TABLE]

Noting that $u(f^{n})(x,0)=0$ and letting $n\to\infty$ in the above equation, we obtain

[TABLE]

In a parallel manner, we also have

[TABLE]

for any $\psi\in C^{1}[0,T]$ with $J^{1-\alpha}_{T^{-}}\psi(T)=0$ and $v\in L^{2}(\Omega)$, which, together with (2.7), implies that $u^{*}(\,\cdot\,,0)=0$.

By $f^{n}\rightharpoonup f^{*}$ in $L^{2}(\Omega)$ and (2.4), we employ the lower semi-continuity of the $L^{2}$-norm to conclude

[TABLE]

indicating that $f^{*}$ is indeed a minimizer to the optimization problem (2.2). Furthermore, the uniqueness of $f^{*}$ is readily seen from the convexity of $J(f)$. ∎

Next, we justify the stability of (2.2), namely, the minimization system (2.2) is indeed a stabilization for Problem 1.1 with respect to the perturbation in observation data.

Theorem 2.2**.**

Let $\{u^{\delta}_{\ell}\}\subset L^{2}(\omega\times(0,T))$ be a sequence such that

[TABLE]

and $\{f^{\ell}\}$ be a sequence of minimizers of problems

[TABLE]

Then $\{f^{\ell}\}$ converges strongly in $L^{2}(\Omega)$ to the minimizer of (2.2).

Proof.

The unique existence of each $f^{\ell}$ is guaranteed by Theorem 2.1. By definition, we have

[TABLE]

which implies the uniform boundedness of $f^{\ell}$ in $L^{2}(\Omega)$. Hence, there exist $f^{*}\in L^{2}(\Omega)$ and a subsequence of $\{f^{\ell}\}$, still denoted by $\{f^{\ell}\}$, such that

[TABLE]

Now it suffices to show that $f^{*}$ is indeed the unique minimizer of (2.2). Actually, repeating the same argument as that in the proof of Theorem 2.1, we can derive

[TABLE]

up to taking a further subsequence. Combining the above convergence with (2.8), we obtain

[TABLE]

Therefore, we get

[TABLE]

For any $f\in L^{2}(\Omega)$, again we take advantage of the lower semi-continuity of the $L^{2}$-norm to deduce

[TABLE]

which verifies that $f^{*}$ is the minimizer of (2.2).

Next, we shall prove $\{f^{\ell}\}$ converges to $f^{*}$ strongly in $L^{2}(\Omega)$ by contradiction. Assuming that it is not true, then we know that $\{\|f^{\ell}\|_{L^{2}(\Omega)}\}$ does not converge to $\|f^{*}\|_{L^{2}(\Omega)}$. As $f^{\ell}\rightharpoonup f^{*}$ in $L^{2}(\Omega)$, by the weak lower semi-continuity of the norm, we have

[TABLE]

Hence, setting $A:=\limsup_{\ell\to\infty}\|f^{\ell}\|_{L^{2}(\Omega)}$, we get

[TABLE]

and there exists a subsequence $\{f^{m}\}$ of $\{f^{\ell}\}$ such that

[TABLE]

Now taking $f=f^{*}$ in (2.10), we find that

[TABLE]

This, together with (2.11), implies that

[TABLE]

which contradicts with (2.9). The proof of Theorem 2.2 is completed. ∎

3 Fully Discrete Finite Element Approximation

In this section we propose a fully discrete finite element method for approximating the nonlinear optimization problem (2.2). We first introduce some appropriate time and space discretization. For the space discretization, we consider a shape regular triangulation $\mathcal{T}_{h}$ of $\Omega$ with a mesh size $h$, consisting of tetrahedral elements (see [2]). Then we introduce the standard nodal finite element spaces of piecewise linear functions

[TABLE]

where $P_{1}(A)$ is the space of linear polynomials on $A$. It is obvious that $V_{h}\subset S_{h}$.

To fully discretize the minimization problem (2.2), we also need the time discretization. To this end, we divide the time interval $[0,T]$ into $M$ subintervals by the equidistant nodal points

[TABLE]

with $t_{m}=m\,\tau$, where $\tau=T/M$ is the step length. Hereinafter, we additionally assume that the noisy observation data $u^{\delta}$ satisfies $u^{\delta}\in C([0,T];L^{2}(\omega))$, so that each $u^{\delta}(\,\cdot\,,t_{m})$ makes sense in $L^{2}(\omega)$.

To discretize the Caputo derivative in time, we refer to [14, 7] and use the following approximation:

[TABLE]

where $d_{j}:=(j+1)^{1-\alpha}-j^{1-\alpha}$ for $j=0,1,\ldots,m$, and $\bar{\partial}_{\tau}^{\alpha}u(x,t_{m+1})$ and $r_{\tau}^{m+1}$ denote the L1 numerical approximation444For historical reasons, the numerical scheme in (3.1) is often referred to as L1 numerical approximation for Caputo fractional derivative. One possible explanation is that in this scheme we replace the first derivative with a first order linear interpolation in each cell $[t_{j},t_{j+1}]$, and is called L1 scheme. This name has nothing to do with the usual notation $L^{1}(\Omega)$ for integrable function in $\Omega$. and local truncation error respectively. The coefficients $\gamma_{j}$ ($j=0,1,\ldots,m+1$) in the discrete convolution quadrature (3.1) are given by

[TABLE]

Following the same line as that in [7], now we arrive at a fully discrete scheme for the forward problem (1.1): given $u_{h}^{0}=0$, find $u_{h}^{m+1}\in V_{h}$ ($m=0,\ldots,M-1$) such that

[TABLE]

or equivalently, find $u_{h}^{m}\in V_{h}$ ($m=1,\ldots,M$) such that

[TABLE]

holds for all $\chi_{h}\in V_{h}$, where $b_{0}:=\Gamma(2-\alpha)\,\tau^{\alpha}$. Here and hereinafter, we understand $f_{h}\in S_{h}$ as some approximation of $f$ and $\mu^{m}=\mu(t_{m})$. To emphasize the dependency, we also denote the solution to (3.2) or (3.3) as $u_{h}^{m}(f_{h})$.

On the basis of (3.2), now we are well prepared to propose the fully discrete finite element approximation of the nonlinear optimization (2.2):

[TABLE]

Here $c_{m}$ ($m=0,1,\ldots,M$) are the coefficients of the composite trapezoidal rule for the time integration over $[0,T]$, i.e., $c_{0}=c_{M}=\frac{1}{2}$ and $c_{m}=1$ for $m=1,\ldots,M-1$.

Before verifying the existence of a minimizer to the discrete minimization problem (3.4), we shall derive some useful a priori estimates for the discrete solution $u_{h}^{m}(f_{h})$ to (3.2) or (3.3). To do so, we need the following lemmas which allows us to use the energy-like methods in the analysis.

Lemma 3.1** (see [23]).**

For the L1 numerical approximation (3.1) of the Caputo derivative, there holds

[TABLE]

Lemma 3.2** (see [11]).**

Suppose that the nonnegative sequences $\{v^{n}\}$ and $\{g^{n}\}$ satisfy

[TABLE]

where $\lambda_{1},\lambda_{2}$ are given positive constants independent of the time step $\tau$. Then there exists a positive constant $\tau_{*}>0$ such that

[TABLE]

for $0<\tau<\tau_{*}$, where $E_{\alpha}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma(1+k\alpha)}$ is the Mittag-Leffler function and $\lambda=\lambda_{1}+\frac{\lambda_{2}}{2-2^{1-\alpha}}$.

Just as the exponential function plays a fundamental role in the integral order equation, the Mittag-Leffler function naturally appears in the fractional order calculus, and it has many nice properties similar to the exponential function. In fact, we have $E_{1}(\pm z)=\mathrm{e}^{\pm z}$ for all $z\in\mathbb{C}$. Let us mention a couple of properties that we will use later, and refer to the recent monograph [4] for further details.

For $\alpha>0$, the Mittag-Leffler function $E_{\alpha}(z)$ is an entire function. In particular, for real variable $\xi\in{\mathbb{R}}$ and $\alpha\in(0,1)$, $E_{\alpha}(\xi)$ is strictly positive and monotone increasing. Moreover, we have $E_{\alpha}(\xi)>0$, $E_{\alpha}(0)=1$,

[TABLE]

where $E_{\alpha,\alpha}(\xi):=\sum_{k=0}^{\infty}\frac{\xi^{k}}{\Gamma(\alpha+k\alpha)}$ is the two parameter Mittag-Leffler function. These properties can be derived from the definition of $E_{\alpha}(\xi)$ and the basic properties of gamma function $\Gamma(\xi)$ for $\xi>0$. In view of the above properties for Mittag-Leffler functions, we know that the function $E_{\alpha}(2\lambda\,t_{n}^{\alpha})$ involved in Lemma 3.2 is greater than one.

Now we proceed to the proof the uniform $H^{1}$-norm estimate for the solution $u_{h}^{m}(f_{h})$ to the fully discrete scheme (3.2) or (3.3).

Lemma 3.3**.**

The fully discrete scheme (3.2) or (3.3) is unconditional stable. Moreover, the following estimates hold for $m=1,\ldots,M$:

[TABLE]

where $C_{0},C_{1}>0$ are constants independent of $h$ and $\tau$.

Proof.

We adopt an inductive argument to prove $(i)$ in (3.6). First, for $m=1$, it follows immediately from (3.2) and $u_{h}^{0}=0$ that

[TABLE]

Substituting $\chi_{h}=u_{h}^{1}(f_{h})$ into the above equation, we obtain

[TABLE]

which implies

[TABLE]

Next, assuming that $(i)$ in (3.6) holds for some $m\geq 1$, we shall prove that it also holds for $m+1$. Choosing $\chi_{h}=u_{h}^{m+1}(f_{h})$ in (3.2) and employing the monotone decreasing property of the sequence $\{d_{j}\}$, we obtain

[TABLE]

By the same argument as that in the proof of [7, Lemma 4.1], we can suitably choose $C_{0}>0$ such that

[TABLE]

holds uniformly for all $m=1,2,\ldots$. Then we immediately obtain $\|u_{h}^{m+1}(f_{h})\|_{L^{2}(\Omega)}\leq C_{0}\|f_{h}\|_{L^{2}(\Omega)}$.

We now prove the assertion $(ii)$ in (3.6). Taking $\chi_{h}=2\bar{\partial}_{\tau}^{\alpha}u_{h}^{m}$ in (3.3) yields

[TABLE]

where

[TABLE]

by (3.5). Then (3.8) leads to

[TABLE]

Hence, we have that

[TABLE]

where we made use of the Poincaré-Friedrichs inequality

[TABLE]

Now applying Lemma 3.2 and (3.7), we obtain that

[TABLE]

This completes the proof. ∎

The next theorem provides the existence of a solution to the discrete system (3.4).

Theorem 3.1**.**

For each fixed $\tau>0$ and $h>0$, there exists at least a minimizer to the discrete system (3.4).

Proof.

It is readily seen from the non-negativity of $J_{h,\tau}(f_{h})$ that $\inf J_{h,\tau}(f_{h})$ is finite. Then there exists a minimizing sequence $\{f_{h}^{n}\}\subset S_{h}$ such that

[TABLE]

Then $\{f_{h}^{n}\}$ is uniformly bounded in $S_{h}\subset L^{2}(\Omega)$. Therefore, there exist $f_{h}^{*}\in S_{h}$ and a subsequence of $\{f_{h}^{n}\}$, still denoted by $\{f_{h}^{n}\}$, such that

[TABLE]

Similarly to the proof of Theorem 2.1, we shall prove that $f_{h}^{*}$ is a minimizer of (3.4).

Since $\tau>0$ is now a fixed constant, if follows from Lemma 3.3 that for each $m=1,\ldots,M$, the sequence $\{u^{m}_{h}(f_{h}^{n})\}$ is uniformly bounded in $V_{h}\subset L^{2}(\Omega)$ with respect to $n$. Again, for each $m=1,\ldots,M$, this indicates the existence of $u^{m}_{h,*}\in V_{h}$ and a subsequence of $\{u^{m}_{h}(f_{h}^{n})\}$, still denoted by $\{u_{h}^{m}(f_{h}^{n})\}$, such that

[TABLE]

In view of the norm equivalence in finite-dimensional spaces, the above two weak convergence results are actually strong, i.e.,

[TABLE]

Now it suffices to show $u^{m}_{h,*}=u_{h}^{m}(f_{h}^{*})$. In fact, by definition (3.2) and up to taking a further subsequence, we see that $\{u_{h}^{m}(f_{h}^{n})\}_{m=1}^{M}$ satisfies

[TABLE]

for all $\chi_{h}\in V_{h}$. Passing $n\to\infty$ in the above equation and using (3.9), we obtain

[TABLE]

for all $\chi_{h}\in V_{h}$. By definition (3.2) again, we conclude $u^{m}_{h,*}=u_{h}^{m}(f_{h}^{*})$ by setting $u_{h,*}^{0}=0$ artificially.

Finally, collecting the above results, we again employ the lower semi-continuity to deduce

[TABLE]

Consequently, $f_{h}^{*}$ is indeed a minimizer of $J_{h,\tau}(f_{h})$ over $S_{h}$. ∎

4 Convergence of the Fully Discrete Approximation

This section is devoted to the convergence analysis of the fully discrete finite element approximation (3.4). In order to relate (3.4) with the continuous minimization problem (2.2), we start with the introduction of a few auxiliary tools and useful results.

First, we recall the standard $L^{2}$ projection $P_{h}$: $L^{2}(\Omega)\to V_{h}$ and Ritz projection $R_{h}$: $H^{1}(\Omega)\to V_{h}$, defined respectively by

[TABLE]

The operators $P_{h}$ and $R_{h}$ satisfy the following approximation properties (see [2, 22]):

[TABLE]

The next lemma is a useful classical approximation result (see [27, 28]).

Lemma 4.1**.**

For a Banach space $X$ and a function $g\in C([0,T];X)$, define a step function approximation by

[TABLE]

where $\mathbf{1}_{(t_{m-1},t_{m}]}$ is the characteristic function of $(t_{m-1},t_{m}]$. Then we have the convergence

[TABLE]

In order to demonstrate the convergence of the finite element approximation (3.4) to the continuous minimization problem (2.2), we shall turn to the following key error estimate, which is a straightforward consequence of [6, Theorem 3.6].

Lemma 4.2**.**

Let $u(f)$ be the solution to system (1.1), where $f\in L^{2}(\Omega)$ and $\mu\in C^{1}[0,T]$. Then for the solution $\{u_{h}^{m}(P_{h}f)\}$ to the fully discrete scheme (3.2) with $f_{h}=P_{h}f$, we have the following error estimate: for any $m=1,\ldots,M$, there holds

[TABLE]

Remark 4.1**.**

According to [6, Theorem 3.6], the source term $f\,\mu$ should satisfy

[TABLE]

for $0<t\leq T$. This is automatically guaranteed by the assumption $f\in L^{2}(\Omega)$ and $\mu\in C^{1}[0,T]$. Therefore, the term including (4.4) with $t=t_{m}$ in [6, Theorem 3.6] is of order $O(\tau)$, which is absorbed in the last term in (4.3). On the other hand, since $t_{m}^{\alpha-1}\tau\leq\tau^{\alpha}\to 0$ as $\tau\to 0$, (4.3) immediately implies

[TABLE]

The following lemma provides a crucial strong convergence.

Lemma 4.3**.**

Let $\{f_{h}\}_{h>0}$ be a sequence in $S_{h}$ which converges weakly to some $f\in L^{2}(\Omega)$ as $h\to 0$. Let $\{u_{h}^{m}(f_{h})\}$ and $u(f)$ be the solutions of (3.2) and (1.1), respectively. Then the following convergence holds:

[TABLE]

Proof.

For simplicity, in this proof we abbreviate $u^{m}:=u(f)(\,\cdot\,,t_{m})$.

Since we have $u(f)\in C([0,T];L^{2}(\Omega))$ by Lemma 2.2, first we can directly apply Lemma 4.1 to obtain

[TABLE]

Using this convergence, (4.5) follows immediately if we can demonstrate

[TABLE]

In order to show (4.6), we utilize the a priori estimates of $u_{h}^{m}(f_{h})$ and $u^{m}$ to deduce

[TABLE]

Clearly, the convergence in (4.6) follows if there holds

[TABLE]

To this end, we split the error $u_{h}^{m}(f_{h})-u^{m}$ into two parts:

[TABLE]

From the estimate (4.3) and Remark 4.1, we can easily see that $\|\theta_{h}^{m}\|_{L^{2}(\Omega)}\to 0$ as $h,\tau\to 0$.

Now we concentrate on the estimate of $\eta_{h}^{m}$. By the definitions of $u_{h}^{m}(f_{h})$ and $u_{h}^{m}(P_{h}f)$, it is readily seen that $\eta_{h}^{m}$ satisfies

[TABLE]

for all $\chi_{h}\in V_{h}$. By the assumption of Lemma 4.3 and (4.1), we have

[TABLE]

which indicates $f_{h}-P_{h}f\rightharpoonup 0$ in $L^{2}(\Omega)$ as $h\to 0$. Meanwhile, by Lemma 3.3 and $f_{h}\,,P_{h}f\in S_{h}$, we estimate

[TABLE]

Then for each $m=1,\ldots,M$, there exist some element $\eta_{*}^{m}\in H^{1}(\Omega)$ and a subsequence of $\{\eta_{h}^{m}\}$, still denoted by $\{\eta_{h}^{m}\}$, such that

[TABLE]

Picking any $\chi\in H^{1}(\Omega)$, we take $\chi_{h}=R_{h}\chi\in V_{h}$ as the test function in (4.8). Then it follows from (4.2) that $\chi_{h}\to\chi$ in $H^{1}(\Omega)$ as $h\to 0$. Together with (4.9), we derive

[TABLE]

Using (4.9) again, we pass $h\to 0$ in equation (4.8) to deduce

[TABLE]

which turns out to be the semidiscrete scheme in time for the original problem (1.1) with $f=0$. Then it follows immediately from [14] that $\|\eta_{*}^{m}\|_{L^{2}(\Omega)}=0$ or equivalently $\eta_{*}^{m}=0$.

Consequently, we obtain $\|\eta_{h}^{m}\|_{L^{2}(\Omega)}\to 0$ as $h\to 0$, and finally conclude

[TABLE]

as $h,\tau\to 0$. This completes the proof of (4.7). ∎

We conclude this section with the convergence of the fully discrete finite element approximation (3.4) to the continuous minimization problem (2.2).

Theorem 4.1**.**

Let $\{f_{h}^{*}\}_{h>0}$ be the minimizers to the discrete minimization problem (3.4). Then there exists a subsequence of $\{f_{h}^{*}\}_{h>0}$ which converges weakly in $L^{2}(\Omega)$ to the minimizer of the continuous problem (2.2) as $h,\tau\to 0$.

Proof.

It is routine to get the boundedness of $\{f_{h}^{*}\}$ in $L^{2}(\Omega)$. Thus there exist $f^{*}\in L^{2}(\Omega)$ and a subsequence of $\{f_{h}^{*}\}$, still denoted by $\{f_{h}^{*}\}$, such that $f_{h}^{*}\rightharpoonup f^{*}$ in $L^{2}(\Omega)$ as $h\to 0$.

Now it suffices to show that $f^{*}$ is the minimizer of the continuous problem (2.2). For any $f\in L^{2}(\Omega)$, we have from (4.1) that

[TABLE]

Since $f_{h}^{*}$ is the minimizer of $J_{h,\tau}(f_{h})$ over $S_{h}$, we obtain

[TABLE]

Then using Lemma 4.3, (4.10), (4.11) and the lower semi-continuity of the $L^{2}$-norm, we deduce

[TABLE]

which implies that $f^{*}$ is indeed a minimizer of the continuous problem (2.2). ∎

5 Iterative Thresholding Algorithm and Numerical Examples

In this section, we develop an efficient iterative thresholding algorithm to solve the discrete minimization problem (3.2)–(3.4) and present several numerical experiments to show its efficiency and accuracy.

Nearly all effective iterative methods for solving nonlinear optimization problems need the information of the derivatives of the concerned objective functional. We shall first derive the derivative of the discrete nonlinear functional $J_{h,\tau}(f_{h})$. Since $u_{h}^{m}(f_{h})$ is linear with respect to $f_{h}$ in view of (3.2), we have $u_{h}^{m}(f_{h})^{\prime}p_{h}=u_{h}^{m}(p_{h})$ for any $p_{h}\in S_{h}$. Hence, it is straightforward to obtain

[TABLE]

Needless to say, it is extremely expensive to use the above formula directly to evaluate the derivatives, since computing the derivative at one fixed point $f_{h}$ needs to solve equation (3.2) once for every direction $p_{h}\in S_{h}$. In order to reduce the computational costs for computing the derivatives, we recall the backward Riemann-Liouville derivative

[TABLE]

and introduce the following adjoint system

[TABLE]

where $\mathbf{1}_{\omega}$ denotes the characterization function of $\omega$. As before, we still denote the solution of (5.2) as $v(f)$ to emphasize its dependency upon $f$.

Lemma 5.1**.**

Let $f,p\in L^{2}(\Omega)$ and $u,v$ be the solutions of systems (1.1) and (5.2) respectively. Then there holds

[TABLE]

Proof.

From the variational system (2.3), we have for any $\varphi\in H^{\alpha}(0,T;L^{2}(\Omega))\cap L^{2}(0,T;H^{1}(\Omega))$ that

[TABLE]

On the other hand, it is routine to get the variational system associated with (5.2) that for any $\psi\in H^{\alpha}(0,T;L^{2}(\Omega))\cap L^{2}(0,T;H^{1}(\Omega))$,

[TABLE]

Now letting $\varphi=v(f)$ and $\psi=u(p)$ in (5.3) and (5.4) respectively, we obtain

[TABLE]

Further, by the integration by parts, the relation (2.1) and noting that

[TABLE]

we have

[TABLE]

Substituting the above equality into (5.5), we arrive at the desired equality. ∎

Writing the above equality into its discrete counterpart as that in previous sections, we have

[TABLE]

Therefore, we can rewrite (5.1) as

[TABLE]

Since $p_{h}\in S_{h}$ was chosen arbitrarily, the above equality gives the following necessary condition for some element $f_{h}^{*}\in S_{h}$ to be a minimizer of the discrete optimization problem (3.4):

[TABLE]

Following the same line as that in [3, 5], this variational equation can be solved by the following iterative thresholding algorithm

[TABLE]

where $L>0$ is a tuning parameter for the convergence. The choice of $L$ should be larger than the operator norm of the forward operator

[TABLE]

and we refer to [5] for details. It is obvious that at each step of (5.6), we only need to solve equation (3.2) for $u_{h}^{m}(f_{h}^{k})$ and then the discrete formulation of system (5.2) for $v_{h}^{m}(f_{h}^{k})$ subsequently, which does not involve heavy computational costs. By the way, although there appears the backward Riemann-Liouville derivative $D_{T-}^{\alpha}$ in (5.2), we know that the solution $v(f)$ coincides with the following problem with a backward Caputo derivative

[TABLE]

where

[TABLE]

Therefore, instead of (5.2) we just discretize system (5.7) to compute $v^{m}_{h}(f_{h})$ by almost the same discrete method for $u^{m}_{h}(f_{h})$.

We are now ready to propose the iterative thresholding algorithm for the reconstruction.

Algorithm 5.1**.**

Choose a tolerance $\epsilon>0$, a regularization parameter $\beta>0$ and a tuning parameter $L>0$. Give an initial guess $f_{h}^{0}$, and set $k=0$.

Compute $f_{h}^{k+1}$ according to the iterative update (5.6). 2. 2.

If $\|f_{h}^{k+1}-f_{h}^{k}\|_{L^{2}(\Omega)}/\|f_{h}^{k}\|_{L^{2}(\Omega)}<\epsilon$, stop the iteration. Otherwise, update $k\leftarrow k+1$ and return to Step 1.

Next we shall present several numerical experiments by applying Algorithm 5.1 to solve Problem 1.1 for $d=1,2$. In general, we set $\Omega=(0,1)^{d}$, $T=1$ and specify various parameters involved in Algorithm 5.1 as follows. Without loss of generality, we always choose the initial guess $f_{h}^{0}$ as a constant, e.g., $f_{h}^{0}\equiv 2$. With the true source term $f^{*}$ and thus the noiseless data $u(f^{*})$, the noisy data $u^{\delta}$ is generated as

[TABLE]

where $\mathrm{rand}(-1,1)$ denotes the random number uniformly distributed on $[-1,1]$. For simplicity, in all examples we set the regularization parameter $\beta=10^{-4}$. Besides the illustrative figures, we mainly evaluate the numerical performance of Algorithm 5.1 by the number $K$ of iterations and the relative $L^{2}$ error

[TABLE]

where we understand $f_{h}^{K}$ as the result of the numerical reconstruction.

In the one-dimensional case, we divide the space-time region $\overline{\Omega}\times[0,T]=[0,1]^{2}$ into $40\times 40$ equidistant meshes. We set the tuning parameter $L$, the stopping criteria $\epsilon$ and the known component $\mu(t)$ as $L=1$, $\epsilon=2\times 10^{-3}$ and $\mu(t)=5+10\,t$ respectively in Algorithm 5.1. Except for the factors mentioned above, we shall test the numerical performance with different choices of exact source terms $f^{*}$, fractional orders $\alpha$, noise levels $\delta$ and observation subdomains $\omega$.

Example 5.1**.**

In this example, we fix the observation subdomain $\omega$ and the noise level $\delta$ as $\omega=\Omega\setminus[1/20,19/20]$ and $\delta=1\%$, respectively. We test Algorithm 5.1 with different fractional orders $\alpha$ and exact source terms $f^{*}$ as follows:

1. (a)

$\alpha=0.3$, $f^{*}(x)=\sin(\pi x/2)+x^{2}+1$. 2. (b)

$\alpha=0.5$, $f^{*}(x)=\sin(\pi x)-2$.

Figure 1 compares the true source terms $f^{*}$ with the corresponding reconstructions $f_{h}^{K}$, and also shows the iteration steps $K$ and the relative errors in the caption.

Example 5.2**.**

In this example, we fix $\alpha=0.8$ and the true source term $f^{*}(x)=-\sin(\pi x)+x+4$, and test different combinations of the observation subdomain $\omega$ and the noise level $\delta$ to observe their influence on the numerical performance. We first fix the noise level $\delta=1\%$ and shrink the size of $\omega$ from $\Omega\setminus[1/10,9/10]$, $\Omega\setminus[1/20,19/20]$ to $\Omega\setminus[1/40,39/40]$. Next, we keep $\omega=\Omega\setminus[1/20,19/20]$ and enlarge the noise level $\delta$ from $0.5\%$, $1\%$, $2\%$ to $4\%$. The resulting numerical performance of the reconstructions is listed in Table 1.

We can see from Figures 1 that with different fractional orders $\alpha$ and $1\%$ noise in the observation data, the numerical reconstructions $f_{h}^{K}$ appear to be quite satisfactory in view of the ill-posedness of Problem 1.1, regardless of the casual choice of the initial guess $f_{h}^{0}$. In addition, we can observe from Table 1 that Algorithm 5.1 has two important advantages, namely, its strong robustness against the oscillating noise in observation data, and its insensitivity to the smallness of the observation subdomain $\omega$.

Now we consider the more challenging two-dimensional case, where we divide the space-time region $\overline{\Omega}\times[0,T]=[0,1]^{3}$ into $40\times 40\times 20$ equidistant meshes. Here we set the tuning parameter $L$ and the known component $\mu(t)$ in Algorithm 5.1 as $L=2$ and $\mu(t)=1+10\pi\,t^{2}$ respectively. Analogously to the one-dimensional counterpart, we evaluate the numerical performance of Algorithm 5.1 from various aspects, including different combinations of the true source terms, noise levels and observation subdomains.

Example 5.3**.**

Parallel to Example 5.1, first we fix the observation subdomain $\omega$ and the noise level $\delta$ as $\omega=\Omega\setminus[1/10,9/10]^{2}$ and $\delta=1\%$, respectively. We take the stopping criteria $\epsilon=\delta/3$ and specify two pairs of fractional orders and true source terms as follows:

1. (a)

$\alpha=0.3$, $f^{*}(x)=\sin(x_{1})+\sin(x_{2})+1$. 2. (b)

$\alpha=0.5$, $f^{*}(x)=\cos(\pi x_{1})\cos(\pi x_{2})+2$.

In Figure 2, we illustrate the true source terms $f^{*}$ with the corresponding reconstructions $f_{h}^{K}$, and the iteration steps $K$ and the relative errors are shown in the caption.

Example 5.4**.**

In this example we fix $\alpha=0.8$, $f^{*}(x)=\exp((x_{1}+x_{2})/4)+1$ and test Algorithm 5.1 with various choices of noise levels $\delta$ and observation subdomains $\omega$ in a similar manner as that in Example 5.2. Here we set $\epsilon=\delta/5$ as the stopping criteria. For the choice of $\omega$, we not only adjust its size, but also change its coverage of the boundary $\partial\Omega$. The resulting numerical performance of the reconstructions is listed in Table 2.

It is readily seen from Examples 5.3–5.4 that Algorithm 5.1 also works efficiently and accurately in the two-dimensional case. It inherits almost all advantages witnessed in the one-dimensional tests in view of its strong robustness against oscillating noise as well as its insensitivity to the smallness of the observation subdomain.

Acknowledgement  The authors thank Professor Bangti Jin (University College London) for his constructive discussions, and appreciate the valuable comments by the anonymous referees.

Funding information  The first author is supported by National Natural Science Foundation of China (NSFC, Nos. 11871240, 11401241, 11571265) and NSFC-RGC (China-Hong Kong, No. 11661161017). The second author is supported by Japan Society for the Promotion of Science KAKENHI Grant Number JP15H05740. The last author is supported by NSFC (Nos. 11871057 and 11501447).