Reconstruction of the Temporal Component in the Source Term of a (Time-Fractional) Diffusion Equation

TL;DR

This paper develops a stable method for reconstructing a time-dependent source term in a time-fractional diffusion equation from single-point observations, especially when the observation point is outside the source support.

Contribution

It introduces a novel stability analysis and a convergent fixed point iteration for numerical reconstruction of the source term in fractional diffusion equations with limited observation data.

Findings

Established multiple logarithmic stability for the inverse problem.

Developed a convergent fixed point iterative reconstruction method.

Validated the approach with numerical examples.

Abstract

In this article, we consider the reconstruction of $\rho(t)$ in the (time-fractional) diffusion equation $(\partial_t^\alpha-\triangle)u(x,t)=\rho(t)g(x)$ ($0<\alpha \le 1$) by the observation at a single point $x_0$. We are mainly concerned with the situation of $x_0 \notin$ supp g, which is practically important but has not been well investigated in literature. Assuming the finite sign changes of $\rho$ and an extra observation interval, we establish the multiple logarithmic stability for the problem based on a reverse convolution inequality and a lower estimate for positive solutions. Meanwhile, we develop a fixed point iteration for the numerical reconstruction and prove its convergence. The performance of the proposed method is illustrated by several numerical examples.

Equations

$${\partial }_t^{\alpha }f(t):=\{\begin{array}{rlrl}& \frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{f^{\prime }(s)}{(t-s{)}^{\alpha }}\mathrm{d}s,& & 0<\alpha <1,\\ & f^{\prime }(t),& & \alpha =1,\end{array}$$

$$\{\begin{array}{ll}({\partial }_t^{\alpha }-△)u(x,t)=\rho (t)g(x),& x\in {\mathbb{R}}^d,0<t<T,\\ u(x,0)=0,& x\in {\mathbb{R}}^d\end{array}$$

$$\{\begin{array}{ll}({\partial }_t^{\alpha }-△)u(x,t)=\rho (t)g(x),& x\in \Omega ,0<t<T,\\ u(x,0)=0,& x\in \Omega ,\\ u(x,t)=0,& x\in \partial \Omega ,0<t<T.\end{array}$$

$$\Vert \rho {\Vert }_{C[0,T]}\le C{\left|\log \Vert u_{\rho }(x_0,\cdot ){\Vert }_{L^2(0,\infty )}\right|}^{-2}$$

$$\Vert \rho {\Vert }_{C[0,T]}\le C\Vert {\partial }_t^{\alpha }{u}_{\rho }(x_0,\cdot ){\Vert }_{C[0,T]}$$

$$\mathcal{D}((-△{)}^{\gamma })=\left\{f\in L^2(\Omega );\sum _{n=1}^{\infty }|{\lambda }_n^{\gamma }(f,{\phi }_n){|}^2<\infty \right\},(-△{)}^{\gamma }f:=\sum _{n=1}^{\infty }{\lambda }_n^{\gamma }(f,{\phi }_n){\phi }_n,$$

$$\Vert f{\Vert }_{\mathcal{D}((-△{)}^{\gamma })}={\left(\sum _{n=1}^{\infty }|{\lambda }_n^{\gamma }(f,{\phi }_n){|}^2\right)}^{\frac{1}{2}}.$$

$$J^{\beta }f(t):=\{\begin{array}{rlrl}& f(t),& & \beta =0,\\ & \frac{1}{\Gamma (\beta )}{\int }_0^t\frac{f(s)}{(t-s{)}^{1-\beta }}\mathrm{d}s,& & 0<\beta \le 1,\end{array}f\in C[0,\infty ).$$

$$D_t^{\beta }f(t):=(J^{1-\beta }f{)}^{\prime }(t)=\{\begin{array}{rlrl}& \frac{1}{\Gamma (1-\beta )}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^t\frac{f(s)}{(t-s{)}^{\beta }}\mathrm{d}s,& & 0\le \beta <1,\\ & f^{\prime }(t),& & \beta =1.\end{array}$$

$$D_t^{\beta }f(t)=\frac{f(0)}{\Gamma (1-\beta )}{t}^{-\beta }+{\partial }_t^{\beta }f(t),f=J^{\beta }(D_t^{\beta }f),f^{\prime }=D_t^{1-\beta }(D_t^{\beta }f).$$

$${\partial }_t^{{\beta }^{\prime }}{\partial }_t^{\beta }f={\partial }_t^{\beta +{\beta }^{\prime }}f\text{mbox}for0\le {\beta }^{\prime }\le 1-\beta .$$

$$\{\begin{array}{ll}g:\text{mbox}boundedcontinuous& \text{mbox}ifd=1,\\ g:\text{mbox}locallyH\ddot{o}ldercontinuous& \text{mbox}ifd=2,3.\end{array}$$

$$u_{\rho }(x,t)={\int }_0^t{D}_s^{1-\alpha }\rho (s)v(x,t-s)\mathrm{d}s,t>0,$$

$$\{\begin{array}{ll}({\partial }_t^{\alpha }-△)v=0& \text{mbox}in{\mathbb{R}}^d\times (0,\infty ),\\ v=g& \text{mbox}in{\mathbb{R}}^d\times \{0\}.\end{array}$$

$$g\in \mathcal{D}((-△{)}^{\epsilon })\text{mbox}withsome\epsilon >0.$$

$$\{\begin{array}{ll}({\partial }_t^{\alpha }-△)v=0& \text{mbox}in\Omega \times (0,\infty ),\\ v=g& \text{mbox}in\Omega \times \{0\},\\ v=0& \text{mbox}on\partial \Omega \times (0,\infty ).\end{array}$$

$$v(x,t)={\int }_{{\mathbb{R}}^d}{K}_{\alpha }(x-x^{\prime },t)g(x^{\prime })\mathrm{d}{x}^{\prime },0<\alpha \le 1,$$

$$|K_{\alpha }(x,t)|\le \{\begin{array}{ll}{c}_1{t}^{-\frac{\alpha }{2}},& d=1,\\ c_1{t}^{-\alpha }(|\log (t^{-\alpha }|x{|}^2)|+1),& d=2,\\ c_1{t}^{-\alpha }|x{|}^{-1},& d=3.\end{array}$$

$$|K_{\alpha }(x,t)|\le c_1{t}^{-\frac{\alpha d}{2}}\exp \left(-c_2{t}^{-\frac{\alpha }{2-\alpha }}|x{|}^{\frac{2}{2-\alpha }}\right).$$

$$K_{\alpha }(x,t)\sim c_3{t}^{-\frac{\alpha d}{2(2-\alpha )}}|x{|}^{-\frac{d(1-\alpha )}{2-\alpha }}\exp \left(-c_4{t}^{-\frac{\alpha }{2-\alpha }}|x{|}^{\frac{2}{2-\alpha }}\right).$$

$$K_1(x,t)=(4\pi t{)}^{-\frac{d}{2}}\exp \left(-\frac{|x{|}^2}{4t}\right).$$

$$E_{\alpha ,\beta }(z):=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+\beta )},z\in \mathbb{C},\alpha >0,\beta \in \mathbb{R}.$$

$$|E_{\alpha ,\beta }(-\eta )|\le \frac{c_0}{1+\eta },\forall \eta \ge 0.$$

$$v(\cdot ,t)=\sum _{n=1}^{\infty }{E}_{\alpha ,1}(-{\lambda }_n{t}^{\alpha })(g,{\phi }_n){\phi }_n.$$

$$v(x,t)$$

$$+{\int }_0^t(4\pi (t-s){)}^{-\frac{d}{2}}{\int }_{\partial \Omega }\exp \left(-\frac{|x-y{|}^2}{4(t-s)}\right){\partial }_{\nu }v(y,s)\mathrm{d}y\mathrm{d}s,$$

$${\partial }_{t}v(\cdot ,t)=-t^{\alpha -1}\sum _{n=1}^{\infty }{\lambda }_n{E}_{\alpha ,\alpha }(-{\lambda }_n{t}^{\alpha })(g,{\phi }_n){\phi }_n,t>0.$$

$$\Vert {\partial }_{t}v(\cdot ,t){\Vert }_{\mathcal{D}((-△{)}^{\gamma })}^2$$

$$\le t^{\alpha (1-2\gamma )-2}\sum _{n=1}^{\infty }\frac{c_0({\lambda }_n{t}^{\alpha }{)}^{\gamma +\frac{1}{2}}}{1+{\lambda }_n{t}^{\alpha }}{\left|{\lambda }_n^{\frac{1}{2}}(g,{\phi }_n)\right|}^2\le {\left(c_0\Vert g{\Vert }_{H^1(\Omega )}{t}^{\alpha (\frac{1}{2}-\gamma )-1}\right)}^2.$$

$$|{\partial }_{t}v(x_0,\cdot )|\le \Vert {\partial }_{t}v(\cdot ,t){\Vert }_{C(\overline{\Omega })}\le c_{\gamma }{c}_0\Vert g{\Vert }_{H^1(\Omega )}{t}^{\alpha (\frac{1}{2}-\gamma )-1},t>0,$$

$${\mathcal{U}}_N:=\{f\in C^1[0,T];$$

Full text

[TABLE]

**AMS Subject Classifications ** 35R11, 35R30, 26A33, 26D10, 65M32

††footnotetext:

Manuscript last updated: .

∗

Graduate School of Mathematical Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8914, Japan. E-mail: ykliu@ms.u-tokyo.ac.jp

†

Corresponding author. Department of Mathematics, Texas A&M University, College Station, TX 77843, USA. E-mail: zhidong@math.tamu.edu

@hangfrom@seccntformat

sectionIntroduction

@xsect

2.3ex plus.2ex

Let $0<\alpha\leq 1$ and $\partial_{t}^{\alpha}$ denote the Caputo derivative defined as

[TABLE]

where $\Gamma(\,\cdot\,)$ is the usual Gamma function. Let $d=1,2,3$ be the spatial dimensions and $\Omega\subset\mathbb{R}^{d}$ be an open bounded domain with a smooth boundary (e.g., of $C^{2}$-class). In this paper, we consider the Cauchy problem for a (time-fractional) diffusion equation

[TABLE]

as well as the corresponding initial-boundary value problem

[TABLE]

Here $\triangle:=\sum_{j=1}^{d}\partial_{x_{j}}^{2}$ denotes the Laplacian in space. The conditions on the temporal component $\rho$ and the spatial component $g$ in the source term will be specified later.

With $\alpha=1$, the governing equation in (1.2) and (1.3) is well known as an orthodox model for diffusion and thermal conduction processes. Since the 1990s, its fractional counterpart with $0<\alpha<1$ has witnessed its significance as a candidate for modeling the anomalous diffusion phenomena in heterogeneous media (see, e.g., [6, 8]). From then on, the time-fractional parabolic operator $\partial_{t}^{\alpha}-\triangle$ and its further variants have gathered increasing attentions within mathematicians. Here we do not intend to provide a complete list of bibliographies, but just refer to [4, 18, 23, 7, 14] which reveal the fundamental properties of time-fractional diffusion equations. It turns out that the fractional case resembles its integer prototype to a certain extent, whereas diverges remarkably in such senses as weaker smoothing effect and unique continuation property. Nevertheless, concerning the strong maximum principle, the strict positivity of the solution with realistic initial data is shared for any $0<\alpha\leq 1$ according to [17, 9]. As a direct application of the established maximum principle, the uniqueness of the following inverse source problem was also immediately proved.

Although [17, 9] only dealt with the initial-boundary value problem (1.3) with $0<\alpha<1$, the same argument simply works for all other cases, and we omit the details here. As natural sequels after the uniqueness, we are interested in the theoretical stability as well as the numerical treatment of Problem 1.1, which are the main focuses of this paper.

As was explained in [17, 16], the formulation of Problem 1.1 stands for the identification of the time evolution pattern $\rho(t)$ of the (anomalous) diffusion by the observation data taken at a single monitoring point $x_{0}$. The time-independent component $g(x)$ of the source term models the spatial distribution of the contaminant source, which is assumed to be given. Such a problem setting is applicable e.g.​ for the Chernobyl disaster and the Fukushima Daiichi nuclear disaster, in which the location of the source is exactly known, but the decay of the radiative strength in time is unknown. In these cases, the contaminant source is harmful and it is extremely dangerous or even impossible to measure inside the support of $g$. Therefore, a more desirable situation in practice is to carry out the observation far away from the location of the source, that is, $x_{0}\not\in\mathrm{supp}\,g$.

Unfortunately, the above practical setting causes considerable difficulties in mathematical analysis, and the related literature is rather limited even in the case of $\alpha=1$. Actually, if $x_{0}\in\mathrm{supp}\,g$ or simply $g(x_{0})\neq 0$, one can reduce Problem 1.1 into a Volterra equation of the second kind, by which the unique solvability follows immediately (see [21, §1.5]). In the case of $x_{0}\not\in\mathrm{supp}\,g$, Cannon and Pérez Esteva [2] proved the following logarithmic stability

[TABLE]

for $\|u_{\rho}(x_{0},\,\cdot\,)\|_{L^{2}(0,\infty)}<1$ and $d=1,2$, where the known component $g$ was restricted as the characterization function $\chi_{\omega}$ of a subdomain $\omega\subset\Omega$. However, the data should be taken in the unrealistic $(0,\infty)$ interval. By assuming the finite sign changes of $\rho$, Saitoh, Tuan and Yamamoto [22] asserted the Hölder stability for Problem 1.1 with a small extra interval $(T,T+\delta)$ of observation data. However, there seems to be several fatal flaws in the arguments in [22], so that the essential ill-posedness of the problem was much underestimated. Following the same line, in this paper we attempt to fill the gaps in [22] and generalize the result to the fractional case. Regarding other inverse problems on determining time-dependent coefficients in parabolic equations, we refer to [1, 3].

Concerning inverse problems for time-fractional diffusion equations, publications for time-dependent unknown coefficients are rather minor compared with that for space-dependent ones. On Problem 1.1, Sakamoto and Yamamoto [23] proved

[TABLE]

by assuming $g(x_{0})\neq 0$. Allowing $g=g(x,t)$ with $|g(x_{0},\,\cdot\,)|\geq g_{0}$ a.e.​ in $(0,T)$ for a constant $g_{0}>0$, Fujishiro and Kian [5] obtained the similar Lipschitz stability as (1.4) in the $L^{p}(0,T)$ norm. Recently, a numerical reconstruction method for $\rho$ was proposed in [24] by the partial boundary measurements. It turns out that all the above mentioned results rely heavily on some non-vanishing assumptions such as $x_{0}\in\mathrm{supp}\,g$, and the requirements on the data regularity were quite strong. In the case of $x_{0}\not\in\mathrm{supp}\,g$, only [17, 16] showed the uniqueness by the strong maximum principle, but there seems no published results on the stability and the numerical inversion to the best knowledge of the authors.

The present paper is mainly motivated by the settings in [22], namely, the assumption on finite sign changes of $\rho$ and an extra interval of observation data. As a correction to the flaws therein, we prove a new reverse convolution inequality (see Lemma 3.1), which is the key to the stability. Then we employ the fundamental solutions of (time-fractional) diffusion equations to give lower estimates of positive solutions to homogeneous problems, which enables us to establish the extremely weak stability of multiple logarithmic type (see Theorem 2.6). Numerically, we develop a fixed point iteration for the numerical reconstruction, and provide a convergence result.

The rest of this paper is organized as follows. Recalling some facts in fractional calculus and (time-fractional) diffusion equations, in Section Reconstruction of the Temporal Component in the Source Term of a (Time-Fractional) Diffusion Equation we state the main stability results in Theorems 2.5 and 2.6 concerning Problem 1.1. Then Section Reconstruction of the Temporal Component in the Source Term of a (Time-Fractional) Diffusion Equation is devoted to the proofs of a key lemma and the above theorems. In Section Reconstruction of the Temporal Component in the Source Term of a (Time-Fractional) Diffusion Equation, we propose an iteration method (4.10) for the numerical treatment of Problem 1.1 and prove Theorem 4.1 for its convergence. We implement several numerical examples in Section Reconstruction of the Temporal Component in the Source Term of a (Time-Fractional) Diffusion Equation to illustrate the performance of the proposed method, and summarize the paper in Section Reconstruction of the Temporal Component in the Source Term of a (Time-Fractional) Diffusion Equation with some concluding remarks.

@hangfrom@seccntformat

sectionPreliminaries and Main Results

@xsect

2.3ex plus.2ex

We start from the introduction of general notations and terminologies. First, we introduce the usual $L^{2}(\Omega)$ space with the inner product $(\,\cdot\,,\,\cdot\,)$ and the Sobolev spaces $H_{0}^{1}(\Omega)$, $H^{2}(\Omega)$, etc.. Next, let $\{(\lambda_{n},\varphi_{n})\}_{n=1}^{\infty}$ be the eigensystem of the elliptic operator $-\triangle:H^{2}(\Omega)\cap H_{0}^{1}(\Omega)\longrightarrow L^{2}(\Omega)$ such that $0<\lambda_{1}<\lambda_{2}\leq\cdots$, $\lambda_{n}\uparrow\infty$ ($n\to\infty$) and $\{\varphi_{n}\}\subset H^{2}(\Omega)\cap H_{0}^{1}(\Omega)$ forms a complete orthonormal basis of $L^{2}(\Omega)$. Then the fractional power $(-\triangle)^{\gamma}$ can be further induced for $\gamma\geq 0$ as (see, e.g., [20])

[TABLE]

and $\mathcal{D}((-\triangle)^{\gamma})$ is a Hilbert space equipped with the norm

[TABLE]

Further, we know $\mathcal{D}((-\triangle)^{\gamma})\subset H^{2\gamma}(\Omega)$ for $\gamma>0$ and especially $\mathcal{D}((-\triangle)^{\frac{1}{2}})=H_{0}^{1}(\Omega)$.

For later use, we recall the Riemann-Liouville integral operator

[TABLE]

Then by definition, the Caputo derivative $\partial_{t}^{\beta}$ can be rewritten as $\partial_{t}^{\beta}f=J^{1-\beta}(f^{\prime})$ for $0<\beta\leq 1$. Further, for $f\in C^{1}[0,\infty)$ the Riemann-Liouville derivative is defined by

[TABLE]

The following lemma collects basic properties of fractional derivatives, which can be simply verified by direct calculations.

[TABLE]

Especially, $D_{t}^{\beta}f=\partial_{t}^{\beta}f$ if $f(0)=0$. Moreover, if $f(0)=\partial_{t}^{\beta}f(0)=0$, then

[TABLE]

(a)* Let $u_{\rho}$ be the solution to (1.2), where we assume*

[TABLE]

Then $u_{\rho}$ allows the representation

[TABLE]

where $v$ solves the homogeneous problem

[TABLE]

(b)* Let $u_{\rho}$ be the solution to (1.3), where we assume*

[TABLE]

Then $u_{\rho}$ also allows the representation (2.3), where $v$ solves the homogeneous problem

[TABLE]

For integers $\alpha=1,2,\ldots$, the above lemma is a fundamental property for evolution equations. For $0<\alpha<1$, there appears the Riemann-Liouville derivative of the temporal component $\rho$ in the convolution (2.3), which automatically results from direct calculations. On the other hand, conditions (2.2) and (2.5) on the spatial component $g$ are technical and we refer to [4, 17] for further details.

Now that the inhomogeneous problems (1.2) and (1.3) are related to the homogeneous ones (2.4) and (2.6) respectively via Lemma 2.2, we shall also recall some basic properties of the homogeneous problems for later use. First we start with the whole space case (2.4).

[TABLE]

where the fundamental solution $K_{\alpha}(x,t)$ satisfies the following properties.

(a)* If $0<|x|^{2}\leq t^{\alpha}$, then there exists a constant $c_{1}>0$ depending on $\alpha,d$ such that*

[TABLE]

(b)* If $|x|^{2}\geq t^{\alpha}$, then there exist constants $c_{1}>0$ and $c_{2}>0$ depending on $\alpha,d$ such that*

[TABLE]

(c)* If $|x|>r$ for some fixed $r>0$, then there exist constants $c_{3}>0$ and $c_{4}>0$ depending on $\alpha,d,r$ such that the following asymptotic behavior holds for $t\downarrow 0$:*

[TABLE]

The fundamental function $K_{\alpha}(x,t)$ involves the Fox $H$-function (see [19]), which takes a rather complicated form. For conciseness, we will not scrutinize the $H$-function in detail, but just exploit some of its properties such as (2.9) and (2.10). Moreover, $K_{\alpha}(x,t)$ generalizes the usual heat kernel with $\alpha=1$, i.e.,

[TABLE]

To express the solution to the initial-boundary value problem (2.6), we recall the Mittag-Leffler function defined by

[TABLE]

If $0<\alpha<2$, then there exists a constant $c_{0}>0$ depending only on $\alpha,\beta$ such that

[TABLE]

Now we proceed to the bounded domain case (2.6).

[TABLE]

(b)* Especially, for $\alpha=1$, $v$ also allows the representation*

[TABLE]

where $\partial_{\nu}v:=\nabla v\cdot\nu$ is the normal derivative and $\nu=\nu(y)$ denotes the unit outward normal vector at $y\in\partial\Omega$.

(c)* If the spatial dimension $d=1$ and $g\in H_{0}^{1}(\Omega)$, then $\partial_{t}v(x_{0},\,\cdot\,)\in L^{1}(0,T)$ for any $x_{0}\in\Omega$.*

[TABLE]

Since $g\in H_{0}^{1}(\Omega)=\mathcal{D}((-\triangle)^{\frac{1}{2}})$, we fix any $\gamma\in(\frac{1}{4},\frac{1}{2})$ and use (2.11) to estimate

[TABLE]

By the one-dimensional Sobolev embedding $\mathcal{D}((-\triangle)^{\gamma})\subset C(\overline{\Omega})$ with $\gamma>\frac{1}{4}$, there exists a constant $c_{\gamma}>0$ such that

[TABLE]

indicating $\partial_{t}v(x_{0},\,\cdot\,)\in L^{1}(0,T)$ immediately. ∎

Next we prepare for the statement of the main theorems. Let $0<\alpha\leq 1$ and $T>0$. For given constants $M>0$ and $N=0,1,2,\ldots$, define the admissible set

[TABLE]

The restriction on the number of sign changes follows the same line as that in [22], where the same problem as Problem 1.1 was investigated for heat equations. Nevertheless, here we assume $\rho\in C^{1}[0,T]$ for $D_{t}^{1-\alpha}\rho$ to make sense in $L^{1}(0,T)$.

Now we state the first main result of this paper.

(a)* If $N=0$, then for any $\delta\in(0,T)$, there exists a constant $B_{\delta}>0$ depending only on $x_{0},g$ such that $B_{\delta}\uparrow\infty\ (\delta\downarrow 0)$ and*

[TABLE]

(b)* If $N>0$, then there exist constants $C_{0}>0$ such that for any sufficiently small $\delta>0$, there holds*

[TABLE]

where $B_{\delta}$ is the same constant as that in (a).

By the assumptions in Theorem 2.5 and Lemmas 2.3–2.4, we see that the observation data $u_{\rho}(x_{0},\,\cdot\,)$ makes sense in $L^{1}(0,T)$. Parallel to that in [22], both estimates in Theorem 2.5 are obtained at the cost of a small extra interval $(T-\delta,T)$ of measurements. Such an interval cannot be arbitrarily small because the constant $B_{\delta}\uparrow\infty$ as $\delta\downarrow 0$. Actually, the choice of $B_{\delta}$ will be

[TABLE]

in the proof, where $v$ stands for the solution to the homogeneous problem (2.4) or (2.6).

Let us explain the results in Theorem 2.5 in detail. First, in comparison with [23, 5] where the same problem as Problem 1.1 were treated, the restriction $x_{0}\in\mathrm{supp}\,g$ is removed, indicating that the monitoring point is allowed to be far away from the support of $g$. As was mentioned in Section Reconstruction of the Temporal Component in the Source Term of a (Time-Fractional) Diffusion Equation, such a relaxation is significant in some practical applications. Moreover, estimates (2.15) and (2.16) only require the $L^{1}$-norm of the observation data, while that in [23, 5] involves the Caputo derivative. Obviously, the former is more desirable in treating the noisy data.

Second, Theorem 2.5(a) gives a Lipschitz stability estimate (2.15) for the inverse source problem under the assumption that $\rho$ does not change its sign on $(0,T)$. In this case, the maximum principle (see Luchko [18]) indicates that $u_{\rho}(x_{0},\,\cdot\,)$ also keeps the sign. This is unrealistic because $u_{\rho}(x_{0},\,\cdot\,)$ stands for the measurement error in the context of stability analysis. Nevertheless, in Section Reconstruction of the Temporal Component in the Source Term of a (Time-Fractional) Diffusion Equation we will propose an iterative method producing a monotone sequence converging to the true solution, where the stability of Lipschitz type becomes remarkable for the fast convergence.

Third, in the general case of $N=1,2,\ldots$, i.e., $\rho$ changes signs finite times on $(0,T)$, the estimate (2.16) in Theorem 2.5(b) looks complicated with some extra terms. Especially, the constant multiplier of $\|u_{\rho}(x_{0},\,\cdot\,)\|_{L^{1}(0,T)}$ increases with larger $N$. In order to further treat the estimate (2.16), according to (2.17) we have to give a lower estimate of $\|v(x_{0},\,\cdot\,)\|_{L^{1}(0,\delta)}$ with respect to $\delta$. This is possible in some cases and we collect the result as follows.

(a)* If $x_{0}\in\mathrm{supp}\,g$, further assume $g\in C(\overline{\Omega})$ in the case of (1.3). Then there exist a sufficiently small constant $\delta_{*}>0$ and a constant $C_{1}>0$ such that*

[TABLE]

(b)* In the case of (1.2), additionally assume that $g$ is compactly supported and $x_{0}\not\in\overline{\mathrm{supp}\,g}$. In the case of (1.3), additionally require $\alpha=1$ and*

[TABLE]

Then there exist a sufficiently small $\delta_{*}>0$ and a constant $C_{2}>0$ such that

[TABLE]

Obviously, estimates (2.18) and (2.20) indicate a weak Hölder stability and an extremely weak stability of the multiple logarithmic type, respectively. In this sense, Theorem 2.6 reflects the severe ill-posedness of Problem 1.1 as the observation data $u_{\rho}(x_{0},\,\cdot\,)\to 0$ in $L^{1}(0,T)$. Especially, the stability becomes worse with a larger number $N$ of the sign changes of $\rho$. Roughly speaking, both (2.18) and (2.20) result from a minimization process based on the estimate (2.16), and the minimizer $\delta_{*}>0$ more or less takes the place of an optimal choice of the regularization parameter.

In Theorem 2.6(a), we reconsider the situation of $x_{0}\in\mathrm{supp}\,g$ as the majority of related literature restricted. In this case, [23, Theorem 4.4] and [5, Theorem 1.1] achieve stronger stabilities than our Hölder one (2.18) at the cost of higher regularity assumptions on the observation data.

Remarkably, Theorem 2.6(b) allows us to choose $x_{0}\not\in\mathrm{supp}\,g$ under certain technical conditions. Owing to the power $-2(\frac{2}{\alpha}-1)$ in (2.20), we observe that the stability is slightly improved with smaller $\alpha$, though such an improvement does not change the essential weakness of the stability. This indicates the similarity of the cases of $\alpha=1$ and $0<\alpha<1$ in our problem, unlike the significant difference occurring in many other problems. In fact, such a similarity originates from the exponential decay of $\|v(x_{0},\,\cdot\,)\|_{L^{1}(0,\delta)}$ as $\delta\downarrow 0$ for all $0<\alpha\leq 1$ under the conditions in Theorem 2.6. For this part, we have to take advantage of the fundamental solutions introduced before. Especially, in the bounded domain case, we require $\alpha=1$ because the fundamental solution for (2.6) with $0<\alpha<1$ is not yet available.

@hangfrom@seccntformat

sectionProofs of Theorems 2.5 and 2.6

@xsect

2.3ex plus.2ex

In this section, we prove the above theorems based on the representation (2.3) in Lemma 2.2. Following the argument in the proof of [10, Theorem 2.6], we take the Riemann-Liouville integral operator $J^{1-\alpha}$ in (2.3) to deduce

[TABLE]

Moreover, as long as $u_{\rho}(x_{0},\,\cdot\,)\in L^{1}(0,T)$, Young’s inequality immediately gives

[TABLE]

Therefore, the problem is now reduced to estimating $\rho$ in the convolution (3.1) with $x=x_{0}$. To this end, the following lemma plays a fundamental role.

(a)* If $f_{1}$ keeps its sign on $(\eta,T_{0}+\delta)$, then*

[TABLE]

(b)* If $f_{1}$ only keeps its sign on $(\eta,T_{0})$, then*

[TABLE]

[TABLE]

Without loss of generality, we can assume $f_{1}\geq 0$ on $(\eta,T_{0})$, or otherwise we may consider $-f_{1}$ instead. Similarly to the proof of [22, Theorem 3.1], we estimate $I$ from below as

[TABLE]

(a) If $f_{1}\geq 0$ throughout $(\eta,T_{0}+\delta)$, then obviously

[TABLE]

and thus (3.5) is exactly (3.3).

(b) Considering the possibility of $f_{1}\leq 0$ on $(T_{0},T_{0}+\delta)$, we split $I$ into three parts and continue to estimate

[TABLE]

Since the term $\int_{T_{0}}^{t}|f_{1}(s)|\,f_{2}(t-s)\,\mathrm{d}s$ is the convolution of $|f_{1}|$ and $f_{2}$ starting from $t=T_{0}$, it follows immediately from Young’s inequality that

[TABLE]

and consequently

[TABLE]

The combination of the above inequality with (3.5) finishes the proof of (3.4). ∎

First we mention that both $v(x_{0},\,\cdot\,)$ and $u_{\rho}(x_{0},\,\cdot\,)$ make sense in $L^{1}(0,T)$. Actually, for the whole space case (2.4), $v$ is pointwisely defined according to Lemma 2.3. For the bounded domain case (2.6), it follows from [13, Theorem 2.1] that $v\in L^{1}(0,T;H^{2}(\Omega)\cap H_{0}^{1}(\Omega))$. Since we restrict $d\leq 3$, the Sobolev embedding theorem implies $v(x_{0},\,\cdot\,)\in L^{1}(0,T)$. On the other hand, we incorporate Lemma 2.1 with $\|\rho\|_{C^{1}[0,T]}\leq M$ to estimate for $0<t\leq T$ that

[TABLE]

implying $D_{t}^{1-\alpha}\rho\in L^{1}(0,T)$. Hence, Young’s inequality immediately gives $u_{\rho}(x_{0},\,\cdot\,)\in L^{1}(0,T)$.

Next, since the initial value $g$ of $v$ is assumed to be non-negative and non-vanishing, we know $v>0$ in $\mathbb{R}^{d}\times(0,\infty)$ in the case of (2.4) by Lemma 2.3, and $v>0$ a.e.​ in $\Omega\times(0,\infty)$ in the case of (2.6) by the strong maximum principle established in [17, Theorem 1.2]. Now that $v(x_{0},\,\cdot\,)>0$ a.e.​ in $(0,T)$, it is readily seen that the constant $B_{\delta}$ defined in (2.17) is strictly positive and well-defined for any $\delta>0$.

Now we are in a position to substitute $f_{1}:=\rho$ and $f_{2}:=v(x_{0},\,\cdot\,)$ into Lemma 3.1 and discuss the case of $N=0$ and $N>0$ respectively.

(a) As $\rho\in\mathcal{U}_{0}$ means that $\rho$ keeps the sign on $(0,T)$, we simply take $\eta=0$ and $T_{0}=T-\delta$ ($0<\delta<T$) in (3.3) to obtain

[TABLE]

where we applied (3.2) in the last inequality.

(b) Now let $N=1,2,\ldots$ and $\delta\in(0,1)$ be sufficiently small. Since $\rho\in\mathcal{U}_{N}$, by definition we can assume that $\rho$ changes signs at $0<t_{1}<t_{2}<\cdots<t_{N^{\prime}}<T-\delta$ ($N^{\prime}\leq N$) on $(0,T-\delta)$. For later convenience, additionally set $t_{0}:=0$ and $t_{N^{\prime}+1}:=T-\delta$. Taking $x=x_{0}$ in (3.1) and rearranging according to the sign of $\rho$, we can write

[TABLE]

where

[TABLE]

Applying (3.4) in Lemma 3.1 to the above equality with $\eta=t_{k}$ and $T_{0}=t_{k+1}$, we obtain

[TABLE]

For $\|\rho\|_{L^{1}(t_{k+1},t_{k+1}+\delta)}$, it follows from $\|\rho\|_{C^{1}[0,T]}\leq M$ and the fact that $\rho(t_{k+1})=0$ that

[TABLE]

For $\|I_{j}\|_{L^{1}(t_{k},t_{k+1}+\delta)}$, we change the order of integrals to deduce

[TABLE]

To deal with $\|v(x_{0},\,\cdot\,)\|_{L^{1}(t_{k}-s,t_{k+1}+\delta-s)}$, we shall turn to the explicit solution to give a pointwise estimate of $|v(x_{0},t)|$ for $0<t<T$. To this end, we discuss the whole space case (2.4) and the bounded domain case (2.6) separately.

For the Cauchy problem (2.4), we recall the formula (2.7) in Lemma 2.3 and the assumption (2.2) to bound

[TABLE]

where we apply (2.9) to treat the case of $|x-x_{0}|>t^{\frac{\alpha}{2}}$ and define

[TABLE]

We first use (2.8) to treat $J_{1}$. For $d=1$, we immediately get $J_{1}\leq 2c_{1}$. For $d=2$, we change to the polar coordinate to estimate

[TABLE]

Similarly, for $d=3$ we have

[TABLE]

By a parallel argument, we can estimate $J_{2}$ as

[TABLE]

where $V_{d}$ denotes the volume of a $d$-dimensional unit ball, and we used the definition of Gamma functions. Substituting all the above estimates into (3.9), we see that $|v(x_{0},t)|$ is dominated by a single constant $c_{1}^{\prime}>0$ for all $t>0$, and hence

[TABLE]

For the initial-boundary value problem (2.6), we invoke the explicit solution by the eigensystem expansion in Lemma 2.4(a). Since we assume $g\in\mathcal{D}((-\triangle)^{\varepsilon})$ with some $\varepsilon>0$ in (2.5), we argue similarly as that in [13] to estimate

[TABLE]

where the Mittag-Leffler function $E_{\alpha,1}(-\lambda_{n}t^{\alpha})$ is treated by (2.11). Then it follows from the embedding $\mathcal{D}((-\triangle)^{\frac{d}{4}+\varepsilon})\subset H^{\frac{d}{2}+2\varepsilon}(\Omega)\subset\subset C(\overline{\Omega})$ that

[TABLE]

implying

[TABLE]

Plugging (3.10) and (3.11) into (3.8), we conclude

[TABLE]

Consequently, we substitute the above inequality and (3.7) into (3.6) to deduce

[TABLE]

for $k=0,1,\ldots,N^{\prime}$. Especially, taking $k=0$ in (3.12) yields

[TABLE]

Now we employ an inductive argument to show for $k=1,2,\ldots,N^{\prime}+1$ that

[TABLE]

In fact, (3.14) with $k=1$ is exactly (3.13). Supposing that (3.14) is valid for some $k\leq N^{\prime}$, we can apply (3.12) to bound $\|\rho\|_{L^{1}(0,t_{k+1})}$ as

[TABLE]

which is indeed (3.14) with $k+1$. Taking $k=N^{\prime}+1$ in (3.14) and recalling $t_{N^{\prime}+1}=T-\delta$, we finally arrive at the desired inequality (2.16) with the aid of (3.2). ∎

Next, under the smallness assumption of the observation data, we proceed to the proof of Theorem 2.6 on the basis of key estimates (3.12) and (3.13). The main idea is to choose different parameters $\delta$ to minimize the estimate on each interval where $D_{t}^{1-\alpha}\rho$ keeps the sign.

On the other hand, as before we assume that $D_{t}^{1-\alpha}\rho$ changes signs at $0<t_{1}<t_{2}<\cdots<t_{N^{\prime}}<T$ ($N^{\prime}\leq N$) on $(0,T)$ and additionally set $t_{0}:=0$ and $t_{N^{\prime}+1}:=T-\delta_{*}$. Here $\delta_{*}>0$ is a sufficiently small constant to be determined later.

Write $D:=\|u_{\rho}(x_{0},\,\cdot\,)\|_{L^{1}(0,T)}$ for simplicity, and recall that $D$ was assumed to be sufficiently small. Since the uniqueness (i.e., the case of $D=0$) for Problem 1.1 was proved in [17], henceforth we may assume $D>0$ without loss of generality.

(a) As $x_{0}\in\mathrm{supp}\,g$ and $g\geq 0$ is assumed to be continuous in $\mathbb{R}^{d}$ or $\overline{\Omega}\,$, we know $g(x_{0})>0$. Due to the continuity of the solution, there holds $v(x_{0},t)\geq\frac{g(x_{0})}{2}$ for sufficiently small $t>0$ and thus

[TABLE]

On the first interval $(0,t_{1})$, we substitute (3.15) into (3.13) and use (3.2) to obtain

[TABLE]

where $a_{0}:=\max(\frac{2\,T^{1-\alpha}}{\Gamma(2-\alpha)\,g(x_{0})},M)$. Since $D$ is sufficiently small, we balance the right-hand side of the above inequality by choosing $\delta_{1}:=D^{\frac{1}{3}}$, which gives

[TABLE]

Although the choice of $\delta_{1}$ here differs from the minimizer up to a multiplier, it gives the same asymptotic behavior with respect to $D$, and henceforth we will adopt such a balancing argument for simplicity.

Now we utilize an inductive argument to show for $k=0,1,\ldots,N^{\prime}$ that

[TABLE]

where $a_{k}$ is defined inductively by

[TABLE]

In fact, (3.16) is exactly the case of $k=0$. Supposing that (3.17) holds for some $k-1$, we apply (3.2) and substitute (3.17) with $k-1$ and (3.15) into (3.12) to treat the interval $(t_{k},t_{k+1})$ as

[TABLE]

where we used the smallness of $D$ to bound $\|u_{\rho}(x_{0},\,\cdot\,)\|_{L^{1}(t_{k},t_{k+1}+\delta)}\leq D\leq D^{(\frac{2}{3})^{k}}$, and $a_{k}$ takes the form of (3.18). Here we choose

[TABLE]

to balance the right-hand side of the above inequality, so that

[TABLE]

as claimed in (3.17).

Finally, we take $k=N^{\prime}$ in (3.17) and recall $t_{N^{\prime}+1}=T-\delta_{*}$ to conclude

[TABLE]

According to (3.19), this is achieved by choosing

[TABLE]

and it follows from the smallness of $D$ that $\delta_{*}$ is also sufficiently small. However, since the sign changes of $\rho$ is independent of $\delta_{*}$, it is possible that the interval $(T-\delta_{*},T)$ includes several of the points $t_{k}$ ($1\leq k\leq N^{\prime}$). In this case, the induction for (3.17) automatically stops before taking $k=N^{\prime}$, indicating that the above estimate is indeed the worst case.

(b) If we do not impose $x_{0}\in\mathrm{supp}\,g$, as before we shall discuss the whole space case (2.4) and the bounded domain case (2.6) separately to estimate $B_{\delta}$.

First we treat the Cauchy problem (2.4). Since we assume the compactness of $\mathrm{supp}\,g$, it is readily seen that $\kappa_{0}=\sup_{x\in\mathrm{supp}\,g}|x-x_{0}|>0$ is finite. Meanwhile, due to $x_{0}\not\in\overline{\mathrm{supp}\,g}$, for sufficiently small $t>0$ we can substitute the asymptotic behavior (2.10) into (2.7) to deduce

[TABLE]

where

[TABLE]

Then we calculate

[TABLE]

To further estimate the above integral from below, we invoke L’Hospital’s rule to see

[TABLE]

for $\beta>0$ and sufficiently large $R>0$. This implies

[TABLE]

which eventually leads us to the lower estimate

[TABLE]

Next we deal with the initial-boundary value problem (2.6) with $\alpha=1$ and the assumption (2.19). Taking $x=x_{0}$ in the representation (2.13) in Lemma 2.4(b), we have

[TABLE]

where

[TABLE]

For $v_{1}(t)$, the same argument as that for (3.20) yields for sufficiently small $t>0$ that

[TABLE]

For $v_{2}(t)$, the Hopf lemma asserts $\partial_{\nu}v<0$ on $\partial\Omega\times(0,T)$ and thus $v_{2}(t)<0$ by definition. Then it is necessary to give an upper bound for $\partial_{\nu}v$ in order to estimate $|v_{2}(t)|$. To this end, we turn to the explicit solution (2.12) again and the trace theorem to dominate

[TABLE]

Therefore, we obtain

[TABLE]

where $\kappa_{1}^{\prime}:=\frac{1}{4}\kappa_{1}^{2}$ and we applied Young’s inequality to the involved convolution. For sufficiently small $t>0$, we use again L’Hospital’s rule to obtain

[TABLE]

indicating

[TABLE]

Combining the above inequality with (3.22) yields

[TABLE]

for sufficiently small $t>0$. Thanks to the assumption (2.19), there holds $\kappa_{0}^{\prime}\leq\kappa_{1}^{\prime}$, implying that the term $t^{-\frac{d}{2}}\exp(-\kappa_{0}^{\prime}\,t^{-1})$ is dominating and thus

[TABLE]

for sufficiently small $t>0$. Repeating the same argument using L’Hospital rule, we conclude

[TABLE]

Noting that the above inequality coincides well with (3.21) for the whole space case, finally we obtain a uniform lower bound for $\|v(x_{0},\,\cdot\,)\|_{L^{1}(0,\delta)}$ in both cases, which implies the following estimate for $B_{\delta}$:

[TABLE]

Since the growth rate of the term $\exp(\kappa^{\prime}\,\delta^{-\frac{\alpha}{2-\alpha}})$ is much faster than any negative power of $\delta$ as $\delta\downarrow 0$, there exists a constant $\kappa>\kappa^{\prime}$ such that for sufficiently small $\delta>0$, there holds

[TABLE]

Now we can follow the same line as that for the case of $x_{0}\in\mathrm{supp}\,g$ to prove (2.20). We start from the first interval $(0,t_{1})$ by using (3.13) and (3.23) to estimate

[TABLE]

where $b_{0}:=(2(\frac{2}{\alpha}-1))^{2(\frac{2}{\alpha}-1)}\max(\frac{C_{1}^{\prime}\,T^{1-\alpha}}{\Gamma(2-\alpha)},M)$. Since $D$ is sufficiently small, we can choose

[TABLE]

to balance the right-hand side of the above inequality. Then $\delta_{1}$ is also sufficiently small, and consequently

[TABLE]

Now we shall take advantage of the induction to show for $k=0,1,\ldots,N^{\prime}$ that

[TABLE]

where $b_{k}$ is defined inductively by

[TABLE]

In fact, (3.24) is exactly the case of $k=0$. Now suppose that (3.25) and (3.26) hold for some $k-1$, and we are in a position to proceed to the case of $k$. Applying (3.12), (3.23) and (3.2), we obtain

[TABLE]

Then we utilize (3.26) with $k-1$ and the smallness of $D$ to further bound

[TABLE]

This indicates

[TABLE]

where $b_{k}$ is exactly (3.27). As before, we choose

[TABLE]

to balance the right-hand side of the above inequality, which yields

[TABLE]

or equivalently (3.25). In addition, the combination of (3.25) with $k$ and (3.26) with $k-1$ implies (3.26) with $k$ immediately.

As a final step, we take $k=N^{\prime}$ in (3.26) and recall $t_{N^{\prime}+1}=T-\delta_{*}$ to conclude

[TABLE]

where

[TABLE]

Meanwhile, now it is clear that the choice of $\delta_{*}$ should be (3.28) with $k=N^{\prime}$, that is,

[TABLE]

Consequently, noting the fact that $N^{\prime}\leq N$, we arrive at the desired estimate (2.20) and the proof is completed. ∎

@hangfrom@seccntformat

sectionNumerical Reconstruction Method

@xsect

2.3ex plus.2ex

This section is devoted to the establishment of an iteration method for the numerical treatment of Problem 1.1. Henceforth, we will only deal with problem (1.3) in a bounded domain $\Omega$ with $0<\alpha<1$ for technical convenience. We denote the true solution of Problem 1.1 by $\rho_{*}$, and write $u_{*}(t):=u_{\rho_{*}}(x_{0},t)$ as the noiseless observation data. The choice of $x_{0}\in\Omega$ can be arbitrary, i.e., we do not require $x_{0}\in\mathrm{supp}\,g$.

First we mention the coincidence of the current problem with the classical deconvolution problem, which is readily seen by taking $x=x_{0}$ in (3.1):

[TABLE]

Indeed, the deconvolution appears frequently as a standard example of ill-posed problems (see, e.g., [12]). However, regardless of the plenary amount of existing approaches to the deconvolution problem, here we attempt to develop a specialized method taking advantage of the underlying fractional diffusion equation.

Practically, we are only given the discrete measurement data of $u_{*}$ and have to work on the discretized setting. For simplicity, we divide the time interval $[0,T]$ into an equidistant partition $0=s_{0}<s_{1}<\cdots<s_{L}=T$, i.e., $s_{\ell}=\frac{T}{L}\ell$ ($\ell=0,1,\ldots,L$). Then the reconstruction of $\rho_{*}$ reduces to the determination of a vector $(\rho_{*}(s_{0}),\rho_{*}(s_{1}),\ldots,\rho_{*}(s_{L}))\in\mathbb{R}^{L+1}$. As a constraint on this vector, we assume

[TABLE]

where $M>0$ is a given constant. In other words, we require that the true solution $\rho_{*}$ is bounded and vanishes at the origin, which is reasonable in practice. In order to prove the convergence of the iteration method proposed later, we shall interpret $\rho_{*}$ as an infinitely differentiable function. More precisely, we assume

[TABLE]

This is always possible by applying appropriate interpolation methods to the discrete pairs $\{(s_{\ell},\rho_{*}(s_{\ell}))\}_{\ell=0}^{L}$. Correspondingly, the resulting solution $u_{\rho_{*}}$ is also smooth and especially $u_{*}=u_{\rho_{*}}(x_{0},\,\cdot\,)\in C^{\infty}[0,T]$. On the other hand, since the known spatial component $g$ is also given on grid points, we can similarly regard $g$ as a sufficiently smooth function, so that the solution $v$ to the homogeneous problem (2.6) is continuous. This, together with the non-negativity of $g$, gives a constant $K>0$ such that

[TABLE]

Now we are well prepared to propose the fixed point iteration

[TABLE]

where $K>0$ refers to the same constant in (4.3). The convergence of (4.4) is guaranteed by the following theorem.

[TABLE]

[TABLE]

Actually, taking the usual time derivative in (4.1), by definition and $\rho(0)=0$ we have

[TABLE]

Due to the homogeneous initial condition of $u_{\rho}$, Lemma 2.1 indicates $\partial_{t}^{\alpha}u_{\rho}(x_{0},t)=D_{t}^{\alpha}u_{\rho}(x_{0},t)$ and thus (4.5) immediately.

On the other hand, investigating the residue $\rho_{*}-\rho_{m}$ produced by (4.4), formally we can write

[TABLE]

We shall show by induction that

[TABLE]

where

[TABLE]

We start from $m=1$. Thanks to the assumption (4.2), we can employ (4.5) to deduce

[TABLE]

Now assume that (4.6) holds true for some $m=1,2,\ldots$, namely,

[TABLE]

We proceed to the case of $m+1$. By the assumption $\rho_{*}^{(m)}(0)=0$, we calculate

[TABLE]

As $(\rho_{*}-\rho_{m})(0)=0$, again we take advantage of (4.5) to obtain

[TABLE]

Since $\rho_{*}^{(m)}(0)=0$ also gives $\rho_{*}^{(m)}(s)=\int_{0}^{s}\rho_{*}^{(m+1)}(\tau)\,\mathrm{d}\tau$ for any $s>0$, we further calculate

[TABLE]

where $\Phi_{m+1}$ takes the exact form of (4.7) with $m+1$. Here we changed the order of integration in (4.8), and applied the obvious identity

[TABLE]

in (4.9). Consequently, the claim (4.6) holds for all $m=1,2,\ldots$. Meanwhile, the assumption (4.3) implies

[TABLE]

According to the inductive relation (4.7), it is straightforward to obtain

[TABLE]

Applying the above inequality and (4.2) to (4.6), for any $t\in[0,T]$ we estimate

[TABLE]

Therefore, we conclude $\|\rho_{*}-\rho_{m}\|_{C[0,T]}\to 0$ as $m\to\infty$ and the proof is completed. ∎

In practice, we implement the iteration (4.4) in the following way. First, we solve the homogeneous problem (2.6) with the given spatial component $g$. This not only provides us with the data of $v(x_{0},t)$, but also gives the upper bound $K$ in (4.3) which will be used repeatedly in the iteration. Next, from the observation data $u_{*}=u_{\rho_{*}}(x_{0},\,\cdot\,)$ we calculate its Caputo derivative $\partial_{t}^{\alpha}u_{*}$. Since $\rho_{0}=0$ implies $\rho_{1}=\frac{1}{K}\partial_{t}^{\alpha}u_{*}$, we can start the iteration from $m=1$ directly. Owing to the representation (4.5), eventually we can rewrite (4.4) as

[TABLE]

Then it becomes obvious that at each iterative update, the main computational costs only involve the numerical differentiation of $\rho_{m-1}$ and its convolution with $v(x_{0},\,\cdot\,)$. Therefore, throughout the reconstruction procedure, it suffices to solve only one fractional diffusion equation (2.6) offline, and the remaining iteration is accomplished by repeating one-dimensional data manipulations. As a result, one can expect that the proposed iterative update (4.10) is extremely efficient. The treatment for the ill-posedness in the numerical differentiation will be explained in Section Reconstruction of the Temporal Component in the Source Term of a (Time-Fractional) Diffusion Equation.

We conclude this section with stating the main algorithm for the numerical reconstruction of Problem 1.1.

Solve the homogeneous problem (2.6) to obtain $v(x_{0},t)$ and fix a constant $K>0$ satisfying (4.3). 2. 2.

Compute $\rho_{m+1}$ by the iterative update (4.10). 3. 3.

If $\|\rho_{m+1}-\rho_{m}\|_{L^{2}(0,T)}\leq\varepsilon$, then stop the iteration. Otherwise, update $m\leftarrow m+1$ and return to Step 2.

@hangfrom@seccntformat

sectionNumerical Experiments

@xsect

2.3ex plus.2ex

In this section, we implement Algorithm 4.2 proposed above to evaluate its numerical performance by several examples.

We start from the general settings of the numerical reconstruction. Throughout this section, basically we fix the fractional order $\alpha=0.9$ unless specified otherwise. Since both unknown function and observation data are variables in time, the spatial dimension is unimportant, so that it suffices to consider $d=1$ and even take $\Omega=(0,1)$ without loss of generality. We set $T=1$ and choose the spatial component of the source term as

[TABLE]

Then it is readily seen that $g\in H_{0}^{1}(0,1)$. Especially, we take the observation point as $x_{0}=\frac{1}{8}$ so that $x_{0}\not\in\mathrm{supp}\,g$. In this case, we can integrate by parts to further simplify (4.10) as

[TABLE]

which even circumvents the numerical differentiation at each step. Note that in this case, $\partial_{t}v(x_{0},\,\cdot\,)$ makes sense in $L^{1}(0,T)$ owing to Lemma 2.4(c). On the other hand, we set the stopping criteria $\varepsilon=10^{-5}$ in the algorithm, and denotes the reconstruction result as $\overline{\rho}$.

In the numerical experiments, we test Algorithm 4.2 with two choices of true solutions to Problem 1.1, namely, a smooth one

[TABLE]

and a non-smooth one

[TABLE]

To solve the forward problems for the observation data $u_{*}(t)$ and the solution $v$ to the homogeneous problem (2.6), we discretize the space-time region $\overline{\Omega}\times[0,T]=[0,1]^{2}$ into $64\times 128$ equidistant meshes, and utilize the L1 time stepping method (see [11, 15]). After computing the value of $v(x_{0},t)$, we choose $K=0.2$ according to (4.3).

First we evaluate the numerical performance of Algorithm 4.2 with the noiseless data $u_{*}$ generated by the true solutions. In this case, we do not need any special treatments for the numerical differentiation involved in (4.10) or (5.1). The numerical results for the true solutions defined in (5.2) and (5.3) are displayed in Figure 1, which are obtained from the iteration (4.10) after $1348$ and $1305$ iterations respectively. For the smooth solution (5.2), Figure 1 (left) shows the convergence of the sequence $\{\rho_{m}\}_{m=0}^{\infty}$, which confirm Theorem 4.1. Furthermore, we attempt the true solution (5.3) which does not belong to $C^{\infty}[0,T]$. However, Figure 1 (right) illustrates that Algorithm 4.2 still works, which means the smoothness restriction on $\rho$ can be suitably weakened in the numerical reconstruction.

With the noiseless data, it turns out from Figure 1 that Algorithm 4.2 produces surprisingly accurate approximations to the respective true solutions which are almost indistinguishable from each other. However, though we do not add any artificial noise, there should exist certain error in the numerical solution of the forward problems. In view of the multiple logarithmic stability asserted in Theorem 2.6, such tiny error in the observation data may result in tremendous error in the reconstruction. Therefore, the above numerical results suggest the possibility to improve the theoretical stability of Problem 1.1. On the other hand, we also observe an interesting monotonicity property in the iteration. More precisely, the sequences $\{\rho_{m}\}_{m=0}^{\infty}$ generated by Algorithm 4.2 are always monotonically increasing and converge to $\rho_{*}$ in both cases of (5.2) and (5.3). This phenomenon is demonstrated by Figure 2, where we plot the first four steps in the iteration of both cases with different fractional orders $\alpha$. If one can rigorously show the monotonicity of $\{\rho_{m}\}$ under certain conditions, then the accuracy of numerical results can be explained by the Lipschitz stability stated in Theorem 2.5(a), because $\rho_{m+1}-\rho_{m}$ ($m=0,1,\ldots$) does not change sign. Unfortunately, we could not prove such a monotonicity property, and the problem still remains open.

Next, we proceed to deal with the noisy data for the practical purpose. Although the theoretical stability was discussed in an $L^{1}(0,T)$ framework, here for simplicity we just add uniform random noises to the noiseless data to produce the noisy data as

[TABLE]

where $\sigma>0$ is the relative noise level and $\mathrm{rand}(-1,1)$ denotes the random number uniformly distributed in $[-1,1]$. According to the iterative update (4.10), the oscillation in $w_{\sigma}$ will definitely cause huge error in computing its Caputo derivative to obtain $\rho_{1}$, and further influence the subsequent numerical differentiation appearing in the iteration. To this end, we have to slightly reform Algorithm 4.2 in the following way. For $f\in C[0,T]$ and $\epsilon\in(0,T)$, define the extension $\widetilde{f}$ of $f$ on $[-\epsilon,T+\epsilon]$ as

[TABLE]

Then we can introduce the mollification

[TABLE]

where the mollifier $\zeta_{\epsilon}$ is defined by

[TABLE]

Now it is obviously that $\zeta_{\epsilon}\geq 0$, $\int_{-\infty}^{\infty}\zeta_{\epsilon}(t)\,\mathrm{d}t=1$ and $f^{\varepsilon}\in C^{1}[0,T]$. Replacing $u_{*}$ and $\rho_{m}$ by $w_{\sigma}^{\epsilon}$ and $\rho_{m}^{\epsilon}$ respectively in (4.10), we arrive at the regularized iteration

[TABLE]

Substituting (4.10) with (5.4) in Algorithm 4.2, again we evaluate the numerical performance with the true solutions (5.2) and (5.3) and different noise levels $\sigma$. Here we take $\epsilon=\frac{5}{128}$ in the mollifier $\zeta_{\epsilon}$, where $\frac{1}{128}$ is the step length in time in the numerical discretization. The comparisons of true solutions with their reconstructions are shown in Figure 3. In all examples, the regularized iteration (5.4) terminates within $40$ steps.

It reveals from Figure 3 that the simple mollification method effectively reduces the ill-posedness of the numerical differentiation in the iteration. Moreover, we can observe that all of the reconstructed solutions $\overline{\rho}(t)$ match the true solutions more accurately for small $t$. In the rear part of the time, the error tends to increase with respect to $t$. This phenomenon provides evidence for the necessity of an extra observation interval $(T-\delta,T)$ in Theorem 2.5, which indicates the increasing instability near $T$.

@hangfrom@seccntformat

sectionConclusions

@xsect

2.3ex plus.2ex

From both theoretical and numerical aspects, in this paper we study the inverse problem on determining $\rho(t)$ in the (time-fractional) diffusion equation

[TABLE]

by the observation taken along $\{x_{0}\}\times(0,T)$. In existing literature, the majority only dealt with the case of $x_{0}\in\mathrm{supp}\,g$, in which the problem turned out to be moderately ill-posed. If $x_{0}\not\in\mathrm{supp}\,g$, it remains a long-standing open problem even for $\alpha=1$ in spite of its practical significance. Motivated by [22], we restrict sign changes of $\rho$ to be finite, and require an extra interval of observation data. Based on the reverse convolution inequality and the lower bound of positive solutions, we prove stability results for the problem in Theorems 2.5 and 2.6. Especially, we establish the multiple logarithmic stability for the case of $x_{0}\not\in\mathrm{supp}\,g$, which, as far as the authors know, seems to be the first affirmative answer to the problem regardless of its weakness.

Numerically, we develop a fixed point iteration for the reconstruction (see Algorithm 4.2), whose convergence is guaranteed by Theorem 4.1. In order to treat noisy data, we apply the mollification method to stabilize the numerical differentiation involved in the iteration. The efficiency and accuracy of our approach are demonstrated by several numerical experiments.

As a representative, we only consider the simplest (time-fractional) parabolic operator $\partial_{t}^{\alpha}-\triangle$ in this paper, but most probably our argument also works for more general formulations, e.g., replacing $\triangle$ by a uniformly non-degenerate elliptic operator. To this end, we shall investigate the corresponding fundamental solutions to give lower estimates for positive solutions. Other future topics include the improvement of theoretical stability under certain assumptions as well as the monotonicity property of the proposed iteration method.

Acknowledgement  This work is supported by Grant-in-Aid for Scientific Research (S) 15H05740, Japan Society for the Promotion of Science (JSPS). The first author is supported by JSPS Postdoctoral Fellowship for Overseas Researchers, Grant-in-Aid for JSPS Fellows 16F16319 and A3 Foresight Program “Modeling and Computation of Applied Inverse Problems”, JSPS. The second author is supported by NSF Grant DMS-1620138.