Numerical analysis of scattered point measurement-based regularization for backward problems for fractional wave equations

Abstract

In this work, our aim is to reconstruct the unknown initial value from terminal data. We develop a numerical framework on nonuniform time grids for fractional wave equations under the lower regularity assumptions. Then, we introduce a regularization method that effectively handles scattered point measurements contaminated with stochastic noise. The optimal error estimates of stochastic convergence not only balance discretization errors, the noise, and the number of observation points, but also propose an a priori choice of regularization parameters. Finally, several numerical experiments are presented to demonstrate the efficiency and accuracy of the algorithm.

Keywords—backward fractional wave, fully discretization, scattered point measurement, regularization method, stochastic error estimates

MSC2020: 35R11, 35R09, 35B40

1 Introduction

Assuming that $\alpha\in(1,2)$ and that $\Omega\subset\mathbb{R}^{d}$, $d=1,2$, is a bounded domain with sufficiently smooth boundary $\partial\Omega$, we consider the following fractional wave equation:

$$\begin{cases}\partial_{t}^{\alpha}u-\Delta u=f(x,t),&(x,t)\in\Omega\times(0,T),\\ u(x,0)=a_{0}(x),&x\in\Omega,\\ \frac{\partial}{\partial t}u(x,0)=a_{1}(x),&x\in\Omega,\\ u(x,t)=0,&(x,t)\in\partial\Omega\times(0,T),\end{cases}$$ (1.1)

where the operator $\partial_{t}^{\alpha}$ is referred to as the Caputo derivative of order $\alpha$, defined by

$$\partial_{t}^{\alpha}\psi(t):=\frac{1}{\Gamma(2-\alpha)}\int_{0}^{t}(t-\tau)^{1-\alpha}\psi^{\prime\prime}(\tau)d\tau,\quad t>0.$$

Equation (1.1) is one of the most famous fractional differential equations. Here, we pay attention to the reconstruction for initial function $a_{1}(x)$ from observation $u(T)$. This is a continuation of our previous paper [2]. Specifically, numerical framework for forward problems and stochastic convergence for scattered point measurement-based regularization are investigated. Without losing generality, let $a_{0}(x)=0$, $f(x,t)=0$. For nonhomogeneous or $a_{0}(x)\neq 0$ cases, it can be obtained by simple change of variables.

Because there are many monographs and papers on this problem, we provide a brief review. There are two types of numerical framework for forward problems: uniform [5, 6, 11] and nonuniform [16, 9, 8] time grids. The most representative papers for sub-diffusion equations are: ‘An analysis of the L1 scheme for the subdiffusion equation with non-smooth data’ [5] and ‘Sharp error estimate of a nonuniform L1 formula for time-fractional reaction subdiffusion equations’ [9]. The main ideas of them are based on the discrete Laplace transform and the discrete complementary convolution kernels technology, respectively. As one of the most famous classical numerical methods, L1 schemes for diffusion-wave equations are presented, see Table 1. $h$ represents the size of the space grids. $N$ is the number of partitions in time grids $0=t_{0}<\cdots<t_{n}<\cdots<t_{N}=T$. And, $\tau=TN^{-1}$ on uniform grids. L1 schemes [17, 15, 12] are proposed based on two kinds of order reduction methods, and they achieve optimal convergence on suitable time grids. The convergence rate of L1 [15] is better than that in [12, 1], but the numerical analysis of it is based on an uncertified assumption. Under the regularity assumptions of $u$, L1 scheme on nonuniform grids with finite difference method and finite element method are presented in [12, 22]. Furthermore, an equivalent integro-differential problem is considered in [13]. The choice of regularization method in the backward problem is based on the regularity of the initial function. In [6], two numerical schemes on uniform grids are proposed under the assumption $a_{1}\in D((-\Delta)^{\gamma/2})\subset H^{\gamma}(\Omega)$, $\gamma\in[0,2]$, which shows that the initial singularity may affect the space convergence when $1+(\gamma-2)\alpha/2<0$. To fill the gap of the above results, we consider investigating the L1 method on nonuniform grids under the lower regularity assumptions $a_{1}\in D((-\Delta)^{q+\epsilon})$, $q\in(\frac{d}{4},1]$ and $0<\epsilon\ll 1$, see Section 2. It shows that taking the suitable observation time $T$, our scheme reaches the optimal convergence rate $O(N^{-(2-\alpha/2)}+h^{2})$.

Table 1: The convergence rates for existing L1 schemes for diffusion-wave equations (1.1).

scheme	rate	time grids	regularity assumption
L1 [17]	$O(\tau^{3-\alpha}+h^{2})$	uniform	$u(x,t)\in C^{(4,3)}_{x,t}(\Omega\times[0,T])$
L1 [23]	$O(\tau^{3-\alpha}+h^{4})$	uniform	$u(x,t)\in C^{(6,3)}_{x,t}(\Omega\times[0,T])$
L1 [15]	$O(N^{-(3-\alpha)}+h^{2})$	nonuniform	$a_{1}\in D((-\Delta)^{3})$
L1 [12]	$O(N^{-(2-\alpha/2)}+h^{2})$	nonuniform	$u(x,\cdot)\in H^{4}(\Omega)$
L1 [1]	$O(N^{-(2-\alpha/2)}+h)$	nonuniform	$a_{1}\in D((-\Delta)^{2})$
Our L1	$O(N^{-(2-\alpha/2)}+t_{n}^{1-\alpha(1-q-\epsilon)}h^{2})$	nonuniform	$a_{1}\in D((-\Delta)^{q+\epsilon})$

It is well known that the Mittag-Leffler functions $E_{\alpha,\beta}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma{(\alpha k+\beta)}}$, $z\in\mathbb{C}$, $\alpha\in(0,1)\cup(1,2)$, $\beta\in\mathbb{C}$, play an important role in investigating the behavior of the solution for fractional differential equations. The potential existence of real roots from the Mittag-Leffler functions in the case $\alpha\in(1,2)$ makes the solution to backward problems for fractional wave equations non-unique. One must make additional assumptions on the terminal time, initial value, or observation data [20, 21, 24, 19, 4]. Notably, when $\alpha$ is in $(1,\frac{4}{3}]$, additional constraints on the terminal time $T$ are not required, which relaxes the conditions for the stability of the backward problem in our previous work. A Tikhonov regularization method based on scattered observations was proposed. Despite the presence of large observation errors, we can still obtain more precise inversion results by increasing the number of observation points $n$, which is difficult for classical regularization algorithms to achieve. Let $a_{1}^{\sigma}$ be the solution found by regularization methods, where $\sigma$ represents the noise level. Taking the optimal regularization parameters $\rho$, the estimates $Err=\|a_{1}-a_{1}^{\sigma}\|_{L^{2}}$ are presented in the following table.

Table 2: The convergence rates for existing regularization methods for backward diffusion-wave equations.

method	optimal estimate	regularity assumption
Quasi-reversibility [20]	$Err\rightarrow 0$ as $\sigma\rightarrow 0$ ($\frac{\rho(\sigma)}{\sigma}\rightarrow 0$)	$a_{1}\in L^{2}(\Omega)$
Quasi-reversibility [21]	$O(\sigma^{\frac{1}{2}})$	$a_{1}\in H_{0}^{2}(\Omega)$
Tikhonov [19]	$O(\sigma^{\frac{\gamma}{\gamma+2}})$	$a_{1}\in D((-\Delta)^{\gamma/2})$
Quasi-boundary [24]	$O(\sigma^{\frac{\gamma}{\gamma+2}})$	$a_{1}\in D((-\Delta)^{\gamma/2})$
Our method	$O\bigg{(}(\sigma n^{-\frac{1}{2}})^{\frac{8(1+q+\epsilon)}{4(1+q+\epsilon)+d}}+n^{-\frac{4(1+q+\epsilon)}{d}}\bigg{)}^{\frac{q+\epsilon}{2(1+q+\epsilon)}}$	$a_{1}\in D((-\Delta)^{q+\epsilon})$

Our previous work focuses on the theoretical analysis of backward problems. As a continuation, we consider such problems in the numerical framework. Our goal is to give an answer to the question: Is it possible to derive an a priori error estimate, showing the way to balance discretization error, the noise, the regularization parameter, and the number of observation points? Specifically, our innovation points are as follows.

1. The optimal error estimates not only balance discretization error, the noise, and the number of observation points, but also propose an a priori choice of regularization parameters.

2. To our best knowledge, it is the first work considering numerical framework on nonuniform grids under the lower regularity assumptions. We also propose optimal choices of mesh parameter $r_{opt}$ for different application cases on graded meshes $t_{n}=T(n/N)^{r}$, $r\geq 1$. Specifically, $r_{opt}=2$, $1+\frac{2}{2-\alpha}$ for $\Delta a_{1}\in L^{2}(\Omega)$ and $a_{1}\in D((-\Delta)^{\gamma})$, $\gamma\in(\frac{d}{4},1)$, respectively.

3. In our previous work [2], there are two cases of observation time $T$ for the stability of the backward problem to the fractional wave equations when $\alpha$ is in $(1,\frac{4}{3}]$ or $(\frac{4}{3},2)$. In our numerical framework, we find that optimal error estimates can be obtained under the same strategy of observation time $T$.

The structure of this paper is as follows. In Section 2, a L1 numerical framework is proposed for the forward problem. Convergence of it is presented under the lower regularity assumptions. The strategy of observation time $T$ for the backward problem is given based on theoretical results. In Section 3, we introduce a scattered point measurement-based regularization method and derive the optimal error estimates of stochastic convergence under the numerical framework. Regularity assumptions used in Section 2 are confirmed in Section 4. Numerical experiments are carried out to verify the theoretical results in Section 5.

2 Numerical framework for forward problems

Actually, weak singularity of the solution has become an important subject in numerical analysis for fractional diffusion wave equations. Recently, a novel order reduction method [12], called symmetric fractional order reduction(SFOR), was proposed such that the analysis techniques [9, 10] on temporal nonuniform mesh work successfully. Although this method has been applied to numerically solve multiple problems, its feasibility has not been fully considered.

2.1 Feasibility of the SFOR method

The main idea is presented in Lemma 2.1. It has certain requirements for the smoothness of the solution.

Lemma 2.1.

For $\alpha\in(1,2)$ and $u(t)\in C^{1}([0,T])\cap C^{2}((0,T])$, we have

$$\partial_{t}^{\alpha}u(t)=\partial_{t}^{\frac{\alpha}{2}}(\partial_{t}^{\frac{\alpha}{2}}u(t))-u^{\prime}(0)\omega_{2-\alpha}(t),$$

where $\omega_{p}(t)=\frac{t^{p-1}}{\Gamma{(p)}}$, $p\geq 0$.

The proof is given in [12]. Here, we just focus on the key step of it as follows

$$\partial_{t}^{\frac{\alpha}{2}}u(t) =\int_{0}^{t}\omega_{1-\frac{\alpha}{2}}(t-s)u^{\prime}(s)ds$$

$$=-u^{\prime}(s)\omega_{2-\frac{\alpha}{2}}(t-s)|_{0}^{t}+\int_{0}^{t}\omega_{2-\frac{\alpha}{2}}(t-s)u^{\prime\prime}(s)ds$$

$$=u^{\prime}(0)\omega_{2-\frac{\alpha}{2}}(t)+\int_{0}^{t}\omega_{2-\frac{\alpha}{2}}(t-s)u^{\prime\prime}(s)ds.$$

Obviously, the integral $\int_{0}^{t}\omega_{2-\frac{\alpha}{2}}(t-s)u^{\prime\prime}(s)ds$ exists when $|u^{\prime\prime}(t)|\leq C(1+t^{\gamma})$, $\gamma>-1$. It implies that the SFOR method works under the conditions $a_{0}=0$ and $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{d}{4},1]$ and $0<\epsilon\ll 1$, see Theorems 10 and 11. And, $|v(0)|=0$ is derived from $u(t)\in C^{1}([0,T])$. Let $\nu=\alpha/2$, model (1.1) is transformed into the following form

$$\begin{cases}\partial_{t}^{\nu}v-\Delta u=a_{1}(x)\omega_{2-\alpha}(t),&(x,t)\in\Omega\times(0,T),\\ v=\partial_{t}^{\nu}u,&(x,t)\in\Omega\times(0,T),\\ u(x,0)=v(x,0)=0,&x\in\Omega,\\ u(x,t)=v(x,t)=0,&(x,t)\in\partial\Omega\times(0,T).\end{cases}$$ (2.1)

In the next subsections, a semi-discrete time scheme is presented for Eq. (2.1). The stability and convergence of them are derived under some reasonable regularity assumptions. For spatial discretization, we adopt standard Galerkin method with continuous piecewise linear finite elements. Let $u_{h}$ be the space semi-discrete solution with a mesh size $h$. Following the idea in [6], the spatial error $u-u_{h}$ is presented in Theorem 1. When $\alpha(1-\gamma-\epsilon)<1$, prefactor $t^{1-\alpha(1-\gamma-\epsilon)}\rightarrow\infty$, $t\rightarrow 0$, which reflects the initial singularity.

Theorem 1.

If $a_{0}=0$, $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{d}{4},1]$ and $0<\epsilon\ll 1$, then

$$\|u(t)-u_{h}(t)\|_{L^{2}}+h\|\nabla(u(t)-u_{h}(t))\|_{L^{2}}\leq Ch^{2}t^{1-\alpha(1-\gamma-\epsilon)}\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}}.$$

Specifically, the estimate is consistent in global time when $d=2$, that is $\|u(t)-u_{h}(t)\|_{L^{2}}\sim h^{2}$. In this case, graded temporal meshes can effectively resolve the impaction of weak singularity. The convergence analysis is prposed in subsection 2.4.

2.2 Time semi-discrete scheme

Here, the L1 formula (2.2), one of the most classical discrete method, is used to approximate the Caputo derivative at nonuniform meshes $\{t_{n}|0=t_{0}<t_{1}<\dots<t_{N}=T\}$. Denote $\tau_{n}=t_{n}-t_{n-1}$, $n\geq 1$.

$$\bar{\partial}_{t}^{\nu}v(t_{n}): =\sum_{k=1}^{n}\int_{t_{k-1}}^{t_{k}}\omega_{1-\nu}(t_{n}-s)\frac{v(t_{k})-v(t_{k-1})}{\tau_{k}}ds$$

$$=\sum_{k=1}^{n}A_{n-k}^{(n)}\nabla_{\tau}v(t_{k}),$$ (2.2)

where $\nabla_{\tau}v(t_{k})=v(t_{k})-v(t_{k-1})$, $A_{n-k}^{(n)}:=\int_{t_{k-1}}^{t_{k}}\frac{\omega_{1-\nu}(t_{n}-s)}{\tau_{k}}ds$. The corresponding numerical scheme (2.3) is as following:

$$\begin{cases}\bar{\partial}_{t}^{\nu}V^{n}-\Delta U^{n}=a_{1}(x)\omega_{2-\alpha}(t_{n}),&1\leq n\leq N,\\ V^{n}=\bar{\partial}_{t}^{\nu}U^{n},&1\leq n\leq N,\\ U(x,0)=V(x,0)=0,&x\in\Omega,\\ U(x,t)=V(x,t)=0,&(x,t)\in\partial\Omega\times(0,T),\end{cases}$$ (2.3)

where $U$ and $V$ are numerical solutions corresponding to $u$ and $v$ in (2.1). Some important lemmas are introduced.

Lemma 2.2.

([9])For $V^{n}$, $0\leq n\leq N$, one has

$$(\bar{\partial}_{t}^{\nu}V^{n},V^{n})\geq\frac{1}{2}\bar{\partial}_{t}^{\nu}\|V^{n}\|_{L^{2}}^{2}.$$

Lemma 2.3.

([12])Let $(g^{n})_{n=1}^{N}$ and $(\lambda_{l})_{l=0}^{N-1}$ be given nonnegative sequences. Assume that there exists a constant $\Lambda$ such that $\Lambda\geq\sum_{l=0}^{N-1}\lambda_{l}$, and that the maximum step satisfies

$$\max_{1\leq n\leq N}\tau_{n}\leq\frac{1}{\sqrt[\nu]{\Gamma{(2-\nu)}\Lambda}}.$$

Then, for any nonnegative sequence $(v^{k})_{k=0}^{N}$ and $(w^{k})_{k=0}^{N}$ satisfying

$$\sum_{k=1}^{n}A_{n-k}^{(n)}\nabla_{\tau}\big{[}(v^{k})^{2}+(w^{k})^{2}\big{]}\leq\sum_{k=1}^{n}\lambda_{n-k}\big{(}v^{k}+w^{k}\big{)}^{2}+(v^{n}+w^{n})g^{n},~{}~{}1\leq n\leq N,$$

it holds that

$$v^{n}+w^{n}\leq 4E_{\nu}(4\Lambda t_{n}^{\nu})\bigg{(}v^{0}+w^{0}+\max_{1\leq k\leq n}\sum_{j=1}^{k}P_{k-j}^{(k)}g^{j}\bigg{)},~{}~{}1\leq n\leq N,$$

where $E_{\nu}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma{(1+k\nu)}}$ is the Mittag-Leffler function.

Lemma 2.4.

For the sequence $(P_{n-j}^{(n)})_{j=1}^{n}$, some properties are given in [9].

$$\sum_{j=k}^{n}P_{n-j}^{(n)}A_{j-k}^{(j)}\equiv 1,~{}~{}1\leq k\leq n,$$

$$0\leq P_{n-j}^{(n)}\leq\Gamma{(2-\nu)}\tau_{j}^{\nu},~{}~{}\sum_{j=1}^{n}P_{n-j}^{(n)}\omega_{1-\nu}(t_{j})\leq C,~{}~{}1\leq j\leq n\leq N.$$

2.3 Stability

Theorem 2.

If $a_{0}=0$, $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{d}{4},1]$ and $0<\epsilon\ll 1$, then

$$\|V^{n}\|_{L^{2}}+\|\nabla U^{n}\|_{L^{2}}\leq C\bigg{(}\|V^{0}\|_{L^{2}}+\|\nabla U^{0}\|_{L^{2}}+2\|a_{1}\|_{L^{2}}\max_{1\leq k\leq n}\sum_{j=1}^{k}P_{k-j}^{(k)}\omega_{2-\alpha}(t_{j})\bigg{)},~{}~{}1\leq n\leq N.$$

Proof.

Taking the inner product with $V^{n}$ and $\Delta U^{n}$ for the first two equations of (2.1), respectively. It gives that

$$(\bar{\partial}_{t}^{\nu}V^{n},V^{n})-(\Delta U^{n},V^{n}) =\omega_{2-\alpha}(t_{n})(a_{1},V^{n}),$$

$$(V^{n},\Delta U^{n}) =(\bar{\partial}_{t}^{\nu}U^{n},\Delta U^{n}).$$

Adding above equations, it comes that

$$(\bar{\partial}_{t}^{\nu}V^{n},V^{n})+(\bar{\partial}_{t}^{\nu}\nabla U^{n},\nabla U^{n})=\omega_{2-\alpha}(t_{n})(a_{1},V^{n}).$$

By Lemma 2.2, one has

$$\bar{\partial}_{t}^{\nu}(\|V^{n}\|_{L^{2}}^{2}+\|\nabla U^{n}\|_{L^{2}}^{2})\leq 2\omega_{2-\alpha}(t_{n})\|a_{1}\|_{L^{2}}\|V^{n}\|_{L^{2}}$$

$$\leq 2\omega_{2-\alpha}(t_{n})\|a_{1}\|_{L^{2}}(\|V^{n}\|_{L^{2}}+\|\nabla U^{n}\|_{L^{2}}).$$

The desired result follows from Lemma 2.3. ∎

Remark 2.1.

If $d=1$ in Theorem 2, then the $L^{\infty}$ stability is derived by the embedding inequality

$$\|U^{n}\|_{L^{\infty}}\leq C\bigg{(}\|V^{0}\|_{L^{2}}+\|\nabla U^{0}\|_{L^{2}}+2\|a_{1}\|_{L^{2}}\max_{1\leq k\leq n}\sum_{j=1}^{k}P_{k-j}^{(k)}\omega_{2-\alpha}(t_{j})\bigg{)},~{}~{}1\leq n\leq N.$$

Theorem 3.

If $a_{0}=0$ and $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{1}{2},1]$ and $0<\epsilon\ll 1$, then

$$\|\nabla V^{n}\|_{L^{2}}+\|\Delta U^{n}\|_{L^{2}}\leq C\bigg{(}\|\nabla V^{0}\|_{L^{2}}+\|\Delta U^{0}\|_{L^{2}}+2\|\nabla a_{1}\|_{L^{2}}\max_{1\leq k\leq n}\sum_{j=1}^{k}P_{k-j}^{(k)}\omega_{2-\alpha}(t_{j})\bigg{)}$$

is valid for any $1\leq n\leq N$. Here the constant $C>0$ is independent of $n$.

Proof.

Taking the inner product with $-\Delta V^{n}$ and $-\Delta^{2}U^{n}$ for the first two equations of (2.1), respectively. It gives that

$$-(\bar{\partial}_{t}^{\nu}V^{n},\Delta V^{n})+(\Delta U^{n},\Delta V^{n}) =-\omega_{2-\alpha}(t_{n})(a_{1},\Delta V^{n}),$$

$$-(V^{n},\Delta^{2}U^{n}) =-(\bar{\partial}_{t}^{\nu}U^{n},\Delta^{2}U^{n}).$$

Adding above equations, it comes that

$$(\bar{\partial}_{t}^{\nu}\nabla V^{n},\nabla V^{n})+(\bar{\partial}_{t}^{\nu}\Delta U^{n},\Delta U^{n})=\omega_{2-\alpha}(t_{n})(\nabla a_{1},\nabla V^{n}).$$

By Lemma 2.2, one has

$$\bar{\partial}_{t}^{\nu}(\|\nabla V^{n}\|_{L^{2}}^{2}+\|\Delta U^{n}\|_{L^{2}}^{2}) \leq 2\omega_{2-\alpha}(t_{n})\|\nabla a_{1}\|_{L^{2}}\|\nabla V^{n}\|_{L^{2}}$$

$$\leq 2\omega_{2-\alpha}(t_{n})\|\nabla a_{1}\|_{L^{2}}(\|\nabla V^{n}\|_{L^{2}}+\|\Delta U^{n}\|_{L^{2}}).$$

The desired result follows from Lemma 2.3. ∎

Remark 2.2.

When $d=2$, $a_{0}=0$, $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{d}{4},1]$ and $0<\epsilon\ll 1$, from Theorem 3 and the embedding inequality, one gets

$$\|U^{n}\|_{L^{\infty}}\leq C\bigg{(}\|\nabla V^{0}\|_{L^{2}}+\|\Delta U^{0}\|_{L^{2}}+2\|\nabla a_{1}\|_{L^{2}}\max_{1\leq k\leq n}\sum_{j=1}^{k}P_{k-j}^{(k)}\omega_{2-\alpha}(t_{j})\bigg{)},~{}~{}1\leq n\leq N.$$

2.4 Convergence of the case in $\mathbb{R}^{2}$

In this part, we propose a convergence estimate for the numerical scheme (2.3). Our main focus is on the bounded domain $\Omega\subset\mathbb{R}^{2}$. The convergence analysis is considered under the following conditions: $a_{0}=0$ and $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{1}{2},1]$, $0<\epsilon\ll 1$. This is consistent with the feasibility conditions for the SFOR framework. From Theorem 1, we know that the spatial error is consistent in global time for the two-dimensional situation. Graded temporal meshes can be used directly to resolve the initial singularity. However, the case in $\mathbb{R}^{1}$ is more complicated. The choice of time step size may affect the spatial error as $t\rightarrow 0$ if $a_{0}=0$ and $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{1}{4},1]$, $0<\epsilon\ll 1$. Therefore, convergence is discussed in two subsections 2.4-2.5.

Eq. (2.1) becomes the following form at $t_{n}=T(n/N)^{r}$, $r\geq 1$, $0\leq n\leq N$. It is easy to verify that $\tau_{n}\leq CN^{-1}t_{n}^{1-\frac{1}{r}}$, $n\geq 2$. Denote $r_{1}^{n}:=\bar{\partial}_{t}^{\nu}v^{n}-\partial_{t}^{\nu}v^{n}$ and $r_{2}^{n}:=\partial_{t}^{\nu}u^{n}-\bar{\partial}_{t}^{\nu}u^{n}$.

$$\begin{cases}\bar{\partial}_{t}^{\nu}v^{n}-\Delta u^{n}=a_{1}(x)\omega_{2-\alpha}(t_{n})+r_{1}^{n},&1\leq n\leq N,\\ v^{n}=\bar{\partial}_{t}^{\nu}u^{n}+r_{2}^{n},&1\leq n\leq N,\\ u(x,0)=v(x,0)=0,&x\in\Omega,\\ u(x,t)=v(x,t)=0,&(x,t)\in\partial\Omega\times(0,T).\end{cases}$$

Denote $\bar{u}:=u^{n}-U^{n}$ and $\bar{v}:=v^{n}-V^{n}$. One has the error system (2.4):

$$\begin{cases}\bar{\partial}_{t}^{\nu}\bar{v}^{n}-\Delta\bar{u}^{n}=r_{1}^{n},&1\leq n\leq N,\\ \bar{v}^{n}=\bar{\partial}_{t}^{\nu}\bar{u}^{n}+r_{2}^{n},&1\leq n\leq N,\\ \bar{u}(x,0)=\bar{v}(x,0)=0,&x\in\Omega,\\ \bar{u}(x,t)=\bar{v}(x,t)=0,&(x,t)\in\partial\Omega\times(0,T).\end{cases}$$ (2.4)

The estimate of $r_{1}^{n}$ is based on the regularity result of Theorem 14.

Lemma 2.5.

For $r_{1}^{n}$ in (2.4), it holds that

$$\|r_{1}^{n}\|\leq Ct_{n}^{-\nu}N^{-q_{1}}\|a_{1}\|_{L^{2}},~{}~{}q_{1}=\min\{r(1-\nu),2-\nu\}.$$

Proof.

For $n=1$, it holds that

$$\|r_{1}^{1}\|_{L^{2}} \leq C\bigg{\|}\int_{0}^{\tau_{1}}(\tau_{1}-s)^{-\nu}\bigg{(}v^{\prime}(s)-\frac{v^{1}-v^{0}}{\tau_{1}}\bigg{)}ds\bigg{\|}_{L^{2}}$$

$$\leq C\bigg{\|}\int_{0}^{\tau_{1}}(\tau_{1}-s)^{-\nu}\bigg{(}v^{\prime}(s)-\frac{1}{\tau_{1}}\int_{0}^{\tau_{1}}v^{\prime}(y)dy\bigg{)}ds\bigg{\|}_{L^{2}}$$

$$\leq C\tau_{1}^{-1}\bigg{\|}\int_{0}^{\tau_{1}}(\tau_{1}-s)^{-\nu}\int_{0}^{\tau_{1}}\big{(}v^{\prime}(s)-v^{\prime}(y)\big{)}dyds\bigg{\|}_{L^{2}}$$

$$\leq C\tau_{1}^{-1}\int_{0}^{\tau_{1}}(\tau_{1}-s)^{-\nu}\int_{0}^{\tau_{1}}\big{(}\|v^{\prime}(s)\|_{L^{2}}+\|v^{\prime}(y)\|_{L^{2}}\big{)}dyds$$

$$\leq C\tau_{1}^{-1}\int_{0}^{\tau_{1}}(\tau_{1}-s)^{-\nu}\int_{0}^{\tau_{1}}\big{(}\|v^{\prime}(s)\|_{L^{2}}+\|v^{\prime}(y)\|_{L^{2}}\big{)}dyds$$

$$\leq C\tau_{1}^{-1}\int_{0}^{\tau_{1}}(\tau_{1}-s)^{-\nu}\int_{0}^{\tau_{1}}\big{(}s^{-\nu}+y^{-\nu}\big{)}dyds$$

$$\leq C\tau_{1}^{1-2\nu}.$$

For the case $n\geq 2$, one has

$$\|r_{1}^{n}\|_{L^{2}} \leq C\bigg{\|}\sum_{k=0}^{n-1}\int_{t_{k}}^{t_{k+1}}(t_{n}-s)^{-\nu}\bigg{(}v^{\prime}(s)-\frac{v^{k+1}-v^{k}}{\tau_{k+1}}\bigg{)}ds\bigg{\|}_{L^{2}}$$

$$\leq C\sum_{k=0}^{n-1}\bigg{\|}\int_{t_{k}}^{t_{k+1}}(t_{n}-s)^{-\nu}\bigg{(}v^{\prime}(s)-\frac{v^{k+1}-v^{k}}{\tau_{k+1}}\bigg{)}ds\bigg{\|}_{L^{2}}$$

$$:=C\sum_{k=0}^{n-1}r_{1}^{n,k}.$$

Then, we consider the estimate of $r_{1}^{n}$ term by term. For the case $k=0$, it gives that

$$r_{1}^{n,0} \leq C\bigg{\|}\int_{0}^{t_{1}}(t_{n}-s)^{-\nu}\bigg{(}v^{\prime}(s)-\frac{v^{1}-v^{0}}{\tau_{1}}\bigg{)}ds\bigg{\|}_{L^{2}}$$

$$\leq C\int_{0}^{t_{1}}(t_{n}-s)^{-\nu}\|v^{\prime}(s)\|_{L^{2}}ds+C\tau_{1}^{-1}\int_{0}^{t_{1}}(t_{n}-s)^{-\nu}\int_{0}^{t_{1}}\|v^{\prime}(y)\|_{L^{2}}dyds$$

$$\leq C(t_{n}-t_{1})^{-\nu}\int_{0}^{t_{1}}s^{-\nu}ds+C\tau_{1}^{-1}\int_{0}^{t_{1}}(t_{n}-s)^{-\nu}\int_{0}^{t_{1}}y^{-\nu}dyds$$

$$\leq C(t_{n}-t_{1})^{-\nu}\tau_{1}^{1-\nu}+C\tau_{1}^{-\nu}\int_{0}^{t_{1}}(t_{n}-s)^{-\nu}ds$$

$$\leq C(t_{n}-t_{1})^{-\nu}\tau_{1}^{1-\nu}$$

$$\leq Ct_{n}^{-\nu}\tau_{1}^{1-\nu},$$

where $(t_{n}-t_{1})^{-\nu}=t_{n}^{-\nu}(1-t_{1}/t_{n})^{-\nu}\leq Ct_{n}^{-\nu}$.

For the case $k=n-1$, following identity is used

$$v^{\prime}(s)-\frac{v^{n}-v^{n-1}}{\tau_{n}}=\frac{1}{\tau_{n}}\int_{t_{n-1}}^{t_{n}}v^{\prime}(s)-v^{\prime}(y)dy=\frac{1}{\tau_{n}}\int_{t_{n-1}}^{t_{n}}\int_{y}^{s}v^{\prime\prime}(z)dzdy.$$

From Theorem 14, it holds that

$$\bigg{\|}v^{\prime}(s)-\frac{v^{n}-v^{n-1}}{\tau_{n}}\bigg{\|}_{L^{2}}\leq\frac{1}{\tau_{n}}\int_{t_{n-1}}^{t_{n}}\int_{\min\{s,y\}}^{\max\{s,y\}}\|v^{\prime\prime}(z)\|_{L^{2}}dzdy\leq C\tau_{n}\sup_{s\in(t_{n-1},t_{n})}s^{-\nu-1}.$$

And, one has

$$r_{1}^{n,n-1} \leq C\bigg{\|}\int_{t_{n-1}}^{t_{n}}(t_{n}-s)^{-\nu}\bigg{(}v^{\prime}(s)-\frac{v^{n}-v^{n-1}}{\tau_{n}}\bigg{)}ds\bigg{\|}_{L^{2}}$$

$$\leq C\tau_{n}\sup_{s\in(t_{n-1},t_{n})}s^{-\nu-1}\int_{t_{n-1}}^{t_{n}}(t_{n}-s)^{-\nu}ds$$

$$\leq C\tau_{n}^{2-\nu}\sup_{s\in(t_{n-1},t_{n})}s^{-\nu-1}$$

$$\leq C\bigg{(}\frac{\tau_{n}}{t_{n}}\bigg{)}^{2-\nu}t_{n}^{1-2\nu}$$

$$\leq Ct_{n}^{-\nu}N^{-q}t_{n}^{-q/r}t_{n}^{1-\nu}$$

$$\leq Ct_{n}^{-\nu}N^{-q},$$

the inequality holds using $\tau_{n}\leq CN^{-1}t_{n}^{1-\frac{1}{r}}$. For the case $1\leq k\leq n-2$, we have

$$r_{1}^{n,k} \leq C\bigg{\|}\int_{t_{k}}^{t_{k+1}}(t_{n}-s)^{-\nu}\bigg{(}v^{\prime}(s)-\frac{v^{k+1}-v^{k}}{\tau_{k+1}}\bigg{)}ds\bigg{\|}_{L^{2}}$$

$$\leq C\int_{t_{k}}^{t_{k+1}}(t_{n}-s)^{-\nu-1}\bigg{\|}v(s)-\frac{(s-t_{k})v^{k+1}+(t_{k+1}-s)v^{k}}{\tau_{k+1}}\bigg{\|}_{L^{2}}ds$$

$$\leq C\tau_{k+1}^{2}\sup_{s\in(t_{k},t_{k+1})}s^{-\nu-1}\int_{t_{k}}^{t_{k+1}}(t_{n}-s)^{-\nu-1}ds$$

$$\leq C\bigg{(}\tau_{k+1}^{2-\nu}t_{k+1}^{\nu}\sup_{s\in(t_{k},t_{k+1})}s^{-\nu-1}\bigg{)}\bigg{(}\tau_{k+1}^{\nu}t_{k+1}^{-\nu}\int_{t_{k}}^{t_{k+1}}(t_{n}-s)^{-\nu-1}ds\bigg{)}$$

$$\leq C\bigg{(}\tau_{k+1}^{2-\nu}t_{k+1}^{\nu}\sup_{s\in(t_{k},t_{k+1})}s^{-\nu-1}\bigg{)}\bigg{(}\tau_{k+1}^{\nu}\int_{t_{k}}^{t_{k+1}}s^{-\nu}(t_{n}-s)^{-\nu-1}ds\bigg{)}.$$

Denote $\Phi_{k}:=\tau_{k+1}^{2-\nu}t_{k+1}^{\nu}\sup_{s\in(t_{k},t_{k+1})}s^{-\nu-1}\sim N^{-q}$. Hence, we obtain

$$\sum_{k=1}^{n-2}r_{1}^{n,k} \leq C\sum_{k=1}^{n-2}\Phi_{k}\tau_{k+1}^{\nu}\int_{t_{k}}^{t_{k+1}}s^{-\nu}(t_{n}-s)^{-\nu-1}ds$$

$$\leq C\max_{k}\Phi_{k}\bigg{(}\tau_{n+1}^{\nu}\int_{t_{1}}^{t_{n-1}}s^{-\nu}(t_{n}-s)^{-\nu-1}ds\bigg{)}$$

$$\leq Ct_{n}^{-\nu}N^{-q},$$

the last inequality holds by

$$\int_{t_{1}}^{t_{n-1}}s^{-\nu}(t_{n}-s)^{-\nu-1}ds \leq\int_{t_{1}}^{t_{n}/2}s^{-\nu}(t_{n}-s)^{-\nu-1}ds+\int_{t_{n}/2}^{t_{n-1}}s^{-\nu}(t_{n}-s)^{-\nu-1}ds$$

$$\leq C(t_{n}/2)^{1-\nu}(t_{n}-t_{n}/2)^{-\nu-1}+C(t_{n}/2)^{-\nu}\tau_{n}^{-\nu}$$

$$\leq Ct_{n}^{-\nu}\tau_{n}^{-\nu}+Ct_{n}^{-\nu}\tau_{n}^{-\nu}$$

$$\leq Ct_{n}^{-\nu}\tau_{n+1}^{-\nu}.$$

The proof completes based on above estimates. ∎

Following the idea in Lemma 2.5, if $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{1}{2},1]$, $0<\epsilon\ll 1$, from the regularity of $u$ in Theorem 12, we get the estimate for $r_{2}^{n}$ in Lemma 2.6.

Lemma 2.6.

For $r_{2}^{n}$ in (2.4), it holds that

$$\|\nabla r_{2}^{n}\|\leq Ct_{n}^{-\nu}(\tau_{1}+N^{-q_{2}})\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}},~{}~{}q_{2}=\min\{r(1+\epsilon\alpha),2-\nu\}.$$

Theorem 4.

If $a_{0}=0$, $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{1}{2},1]$ and $0<\epsilon\ll 1$, then

$$\|\bar{v}^{n}\|_{L^{2}}+\|\nabla\bar{u}^{n}\|_{L^{2}}\leq C(\tau_{1}+N^{-q_{1}}+N^{-q_{2}})\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}},~{}~{}1\leq n\leq N.$$

Proof.

Taking the inner product with $\bar{v}^{n}$ and $\Delta\bar{u}^{n}$ for the first two equations of (2.4), respectively. It gives that

$$(\bar{\partial}_{t}^{\nu}\bar{v}^{n},\bar{v}^{n})-(\Delta\bar{u}^{n},\bar{v}^{n}) =(r_{1}^{n},\bar{v}^{n}),$$

$$(\bar{v}^{n},\Delta\bar{u}^{n}) =(\bar{\partial}_{t}^{\nu}\bar{u}^{n},\Delta\bar{u}^{n})+(r_{2}^{n},\Delta\bar{u}^{n}).$$

Adding above equations, it comes that

$$(\bar{\partial}_{t}^{\nu}\bar{v}^{n},\bar{v}^{n})+(\bar{\partial}_{t}^{\nu}\nabla\bar{u}^{n},\nabla\bar{u}^{n})=(r_{1}^{n},\bar{v}^{n})+(\nabla r_{2}^{n},\nabla\bar{u}^{n}).$$

By Lemma 2.2, one has

$$\frac{1}{2}\bar{\partial}_{t}^{\nu}(\|\bar{v}^{n}\|_{L^{2}}^{2}+\|\nabla\bar{u}^{n}\|_{L^{2}}^{2}) \leq\|r_{1}^{n}\|_{L^{2}}\|\bar{v}^{n}\|_{L^{2}}+\|\nabla r_{2}^{n}\|_{L^{2}}\|\nabla\bar{u}^{n}\|_{L^{2}}$$

$$\leq(\|r_{1}^{n}\|_{L^{2}}+\|\nabla r_{2}^{n}\|_{L^{2}})(\|\bar{v}^{n}\|_{L^{2}}+\|\nabla\bar{u}^{n}\|_{L^{2}}).$$

The desired result follows from Lemmas 2.3 and 2.4. ∎

Remark 2.3.

In [7], the authors presented the analysis framework for sub-diffusion equations. The convergence of the fully discrete scheme of diffusion-wave equations follows analogously. From the results of Theorems 1 and 4, it implies that the fully consistent $L^{2}$ norm error reaches optimal $O(h^{2}+N^{-(2-\nu)})$ in $\Omega\times(0,T]$, when the time grid parameter $r=\frac{2-\nu}{1-\nu}$.

2.5 Convergence of the case in $\mathbb{R}^{1}$

The framework of convergence is the same as that for Theorem 4. The truncation error $\|r_{1}^{n}\|$, $1\leq n\leq N$ is available from Lemma 2.5. The point is on the estimate $\|\nabla r_{2}^{n}\|$, $1\leq n\leq N$. The regularity of $\|\partial_{t}^{m}(-\Delta)^{\frac{1}{2}}u\|_{L^{2}}$ is derived from Theorem 13. Then, one has the following results.

Lemma 2.7.

For $r_{1}^{n}$, $r_{2}^{n}$ in (2.4) with $\Omega\subset\mathbb{R}^{1}$, it holds that

$$\|r_{1}^{n}\|\leq Ct_{n}^{-\nu}N^{-q_{1}}\|a_{1}\|_{L^{2}},~{}~{}q_{1}=\min\{r(1-\nu),2-\nu\},$$

$$\|\nabla r_{2}^{n}\|\leq Ct_{n}^{-\nu}N^{-q_{3}}\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}},~{}~{}q_{3}=\min\{r(1-\nu/2),2-\nu\}.$$

The convergence estimate is obtained by Theorem 4 and Lemma 2.7.

Theorem 5.

If $a_{0}=0$, $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{1}{4},1]$ and $0<\epsilon\ll 1$, then

$$\|\bar{v}^{n}\|_{L^{2}}+\|\nabla\bar{u}^{n}\|_{L^{2}}\leq C(N^{-q_{1}}+N^{-q_{3}})\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}},~{}~{}1\leq n\leq N.$$

Combining the results of Theorems 1 and 5, we obtain the convergence of the fully discrete solution $u_{h}^{n}$ at $t_{n}$, $n\geq 1$. The conclusion is different from the one in Remark 2.3. Specifically, the spatial discrete error is $O(h^{2})$, $\alpha\in(1,\frac{4}{3}]$, for $\forall t>0$ in Theorem 6.

Theorem 6.

If $a_{0}=0$, $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{1}{4},1)$ and $0<\epsilon\ll 1$, then

$$\|u(t_{n})-u_{h}^{n}\|_{L^{2}}\leq C(h^{2}+N^{-q_{1}}+N^{-q_{3}}),~{}~{}1<\alpha\leq\frac{1}{1-(\gamma+\epsilon)},$$

$$\|u(t_{n})-u_{h}^{n}\|_{L^{2}}\leq C(t_{n}^{1-\alpha(1-\gamma-\epsilon)}h^{2}+N^{-q_{1}}+N^{-q_{3}}),~{}~{}\frac{1}{1-(\gamma+\epsilon)}<\alpha<2.$$

When $\frac{1}{1-(\gamma+\epsilon)}<\alpha<2$, $t_{1}^{1-\alpha(1-\gamma-\epsilon)}h^{2}$ may blow up as the time step parameter $r$ is large. The situation becomes better when $t_{n}$ away from 0.

Remark 2.4.

When $t_{n}>t^{*}$, $t^{*}$ is a fixed positive constant, one has

$$\|u(t_{n})-u_{h}^{n}\|_{L^{2}}\leq C(h^{2}+N^{-q_{1}}+N^{-q_{3}}).$$

Combining the above results, we have an interesting finding. First, let us go back to the problem of initial value reconstruction. The stability of the backward problem is as follows. There are some suggestions of the choices of observation time $T$ for $\alpha\in(1,\frac{4}{3}]$ or $(\frac{4}{3},2)$.

Theorem 7.

([2]) Letting $\alpha\in(1,2)$, $\beta,\gamma\in[0,1]$ be such that $\gamma+\beta\neq 0$, we suppose that the pair $(u,a_{1})$ in the space $L^{2}\left(0,T;H_{0}^{1}(\Omega)\cap H^{2}(\Omega)\right)\times D((-\Delta)^{\beta})$ is a solution of our backward problem, which corresponds to the measurement data $u(\cdot,T)$. If $\|(-\Delta)^{\beta}a_{1}\|_{L^{2}(\Omega)}\leq M$ for some positive constant $M$, then the estimate

$$\|a_{1}\|_{D((-\Delta)^{-\gamma})}\leq C_{T}M^{\frac{1-\gamma}{1+\beta}}\left\|u(\cdot,T)\right\|_{L^{2}(\Omega)}^{\frac{\beta+\gamma}{\beta+1}}.$$

is valid provided that one of the following conditions hold:

1. $\alpha\in(1,\frac{4}{3}]$ and $T>0$;

2. $\alpha\in(\frac{4}{3},2)$ and $T^{\alpha}\not\in\cup_{k=1}^{\infty}\{\frac{t_{1}}{\lambda_{k}},\cdots,\frac{t_{P}}{\lambda_{k}}\}$.

Here $C_{T}>0$ is a constant that is independent of $a_{1}$ and $u$, but may depend on $\alpha,\beta,\gamma,T,M$, and $\Omega$.

Then, the choice of $T$ also plays an important role in the stochastic convergence analysis of the regularization method based on the measurements of the scattering points. Naturally, the question arises of whether the choice of $T$ has an effect on the numerical inversion. Remark 2.5 gives a positive conclusion on this. Specifically, the convergence of the numerical method is guaranteed and optimally achieved under a suitable time grid. Therefore, we call it the preserving numerical framework.

Remark 2.5.

(Preserving numerical framework) If $a_{0}=0$, $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{1}{4},1)$, $0<\epsilon\ll 1$ and $r=1+\frac{2}{2-\alpha}$, then

1. $\|u(t_{n})-u_{h}^{n}\|_{L^{2}}\leq C(h^{2}+N^{-(2-\frac{\alpha}{2})})\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}}$, $\alpha\in(1,\frac{4}{3}]$ and $t_{n}>0$;

2. $\|u(t_{n})-u_{h}^{n}\|_{L^{2}}\leq C(h^{2}+N^{-(2-\frac{\alpha}{2})})\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}}$, $\alpha\in(\frac{4}{3},2)$ and $t_{n}>t^{*}$, $t^{*}$ is a fixed positive constant.

3 Scattered point measurement-based regularization

For stating the Tikhonov regularization method based on scattered point measurement, we collect a set of scattered points $\{x_{i}\}_{i=1}^{n}$ which are such that $x_{i}\neq x_{j}$ for $i\neq j$ and are quasi-uniformly distributed in $\Omega$, that is, there exists a positive constant $B$ such that $d_{\max}\leq Bd_{\min}$ , where $d_{\max}>0$ and $d_{\min}>0$ are defined by

$$d_{\max}=\sup_{x\in\Omega}\inf_{1\leq i\leq n}|x-x_{i}|,\quad d_{\min}=\inf_{1\leq i\neq j\leq n}|x_{i}-x_{j}|.$$

Furthermore, for any $u,v\in C(\overline{\Omega})$ and $y\in\mathbb{R}^{n}$, we define

$$(y,v)_{n}:=\frac{1}{n}\sum_{i=1}^{n}y_{i}v(x_{i}),\quad(u,v)_{n}:=\frac{1}{n}\sum_{i=1}^{n}u(x_{i})v(x_{i}),$$

and the discrete semi-norm

$$\|u\|_{n}:=\left(\sum_{i=1}^{n}\frac{u^{2}(x_{i})}{n}\right)^{\frac{1}{2}},\quad u\in C(\overline{\Omega}).$$

In view of the stability results from the above sections, it follows that the forward operator $S$ is bounded and one-to-one from $X:=D((-\Delta)^{\beta})$ to $H^{2}(\Omega)$. Moreover, let $a^{*}\in X$ be the unknown initial value of the problem (1.1). We assume that the measurement data contains noise and is presented in the following form:

$$m_{i}=(Sa^{*})(x_{i})+e_{i},\quad i=1,2,...,n,$$ (3.1)

where $\{e_{i}\}_{i=1}^{n}$ denote a sequence of random variables that are independent and identically distributed over the probability space. The expectation is such that $\mathbb{E}[e_{i}]=0$, and the variances are bounded by $\sigma^{2}$, that is, $\mathbb{E}[e_{i}^{2}]\leq\sigma^{2}$.

We denote $\mathbf{x}:=(x_{1},x_{2},\cdots,x_{n})$, $\mathbf{m}=(m_{1},m_{2},...,m_{n})^{T}$, and $\mathbf{e}:=(e_{1},e_{2},\cdots,e_{n})$. Then the above term (3.1) can be rephrased as the vector form:

$$\mathbf{m}=(Sa^{*})(\mathbf{x})+\mathbf{e}.$$

Let $S_{\tau,h}$ be the fully discrete approximation of the operator $S$. The forward problem is solved by the numerical scheme (2.3) in finite element space $X_{h}$. We seek a numerical solution, denoted as $a_{n}^{*}$, for the unknown initial value $a^{*}$, utilizing the Tikhonov regularization form:

$$\arg\min_{a\in X_{h}}\|(S_{\tau,h}a)(\mathbf{x})-\mathbf{m}\|_{n}^{2}+\rho_{n}\|a\|^{2}_{X},$$ (3.2)

where $X_{h}\subset D((-\Delta)^{\beta})$ with $\beta\in\mathbb{R}$ and $\rho_{n}>0$ is called a regularization parameter.

Here we present our main theorem on the stochastic convergence of the SFOR framework. Estimates $\|S_{\tau,h}a^{*}-Sa^{*}\|_{L^{2}}^{2}\leq e_{S}:=O(h^{4}+N^{-(4-\alpha)})$ have been discussed in Remarks 2.3 and 2.5. The feasibility of the SFOR method suggests $\beta\in(\frac{d}{4},1)$. The observation time $T$ satisfies the requirements in Theorem 7.

Theorem 8.

Let $a_{n}^{*}\in X_{h}$ with dimensions $N_{h}$ be the unique solution of our Tikhonov regularization form (3.2). Then there exist constants $\rho_{0}>0$ and $C>0$ such that the following estimates

$$\mathbb{E}\bigl{[}\|S_{\tau,h}a_{n}^{*}-Sa^{*}\|^{2}_{n}\bigr{]}\leq C(\rho_{n}+e_{S})\|a^{*}\|^{2}_{X}+C\left(1+\frac{e_{S}}{\rho_{n}}+\frac{N_{h}e_{S}}{\rho_{n}^{1-\frac{d}{4}/(1+\beta)}}\right)\frac{\sigma^{2}}{n\rho_{n}^{\frac{d}{4}/(1+\beta)}},$$

and

$$\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2}_{X}\bigr{]}\leq C\frac{\rho_{n}+e_{S}}{\rho_{n}}\|a^{*}\|^{2}_{X}+C\left(1+\frac{e_{S}}{\rho_{n}}+\frac{N_{h}e_{S}}{\rho_{n}^{1-\frac{d}{4}/(1+\beta)}}\right)\frac{\sigma^{2}}{n\rho_{n}^{1+\frac{d}{4}/(1+\beta)}}.$$

are valid for any $0<\rho_{n}\leq\rho_{0}$. Here, the constant $C$ is independent of $n,a^{*},a_{n}^{*}$.

Proof.

When the parameter $r$ of graded mesh is taken $1+\frac{2}{2-\alpha}$, it holds that

$$\|S_{\tau,h}a^{*}-Sa^{*}\|_{L^{2}}\leq C(h^{2}+N^{-(2-\frac{\alpha}{2})})\|a^{*}\|_{X},$$

details can be found in Remarks 2.3 and 2.5. Applying the Lemma 3.10 in [3] and Lemma 2.4 in [2], we have

$$\|S_{\tau,h}a^{*}-Sa^{*}\|_{n} \leq C\|S_{\tau,h}a^{*}-Sa^{*}\|_{L^{2}}+Ch^{2}\|Sa^{*}\|_{H^{2}}$$

$$\leq C\|S_{\tau,h}a^{*}-Sa^{*}\|_{L^{2}}+Ch^{2}\|a^{*}\|_{L^{2}}$$

$$\leq C\|S_{\tau,h}a^{*}-Sa^{*}\|_{L^{2}}+Ch^{2}\|(-\Delta)^{\gamma}a^{*}\|_{L^{2}}$$

$$\leq C(h^{2}+N^{-(2-\frac{\alpha}{2})})\|a^{*}\|_{X}.$$

Replace $e(h)$ therein by $e_{S}$ in Assumption 2.3, the proof is complete following the idea in Theorem 2.10 [3]. ∎

Remark 3.1.

For Theorem 8, if $e_{S}\leq C\rho_{n}$ and $N_{h}e_{S}\leq C\rho_{n}^{1-\frac{d}{4}/(1+\beta)}$, we have

$$\mathbb{E}\bigl{[}\|S_{\tau,h}a_{n}^{*}-Sa^{*}\|^{2}_{n}\bigr{]}\leq C\rho_{n}\|a^{*}\|^{2}_{X}+\frac{C\sigma^{2}}{n\rho_{n}^{\frac{d}{4}/(1+\beta)}},$$

and

$$\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2}_{X}\bigr{]}\leq C\|a^{*}\|^{2}_{X}+\frac{C\sigma^{2}}{n\rho_{n}^{1+\frac{d}{4}/(1+\beta)}}.$$

Theorem 9 shows that the error estimate is based on key parameters such as the noise level, the regularization parameter, and the number of observation points. Then, the optimal regularization parameter is given in Remark 3.2.

Lemma 3.1.

[18, Theorems 3.3 and 3.4] There exists a constant $C>0$ such that for all $u\in H^{k}(\Omega)$ with $k>\frac{d}{2}$ , the following estimates are valid:

$$\|u\|^{2}_{L^{2}(\Omega)}\leq C\left(\|u\|^{2}_{n}+n^{-\frac{2k}{d}}\|u\|^{2}_{H^{k}(\Omega)}\right),$$ (3.3)

$$\|u\|^{2}_{n}\leq C\left(\|u\|^{2}_{L^{2}(\Omega)}+n^{-\frac{2k}{d}}\|u\|^{2}_{H^{k}(\Omega)}\right).$$

Theorem 9.

Suppose that the unknown initial value $a^{*}\in X=D((-\Delta)^{\beta})$ with $\beta\in(\frac{d}{4},1]$, and for the regularization parameter $\rho_{n}$, if $e_{S}\leq C\rho_{n}$ and $N_{h}e_{S}\leq C\rho_{n}^{1-\frac{d}{4}/(1+\beta)}$, then there exists a constant $C=C(\beta,\rho_{n},T,\Omega)$ such that

$$\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2+\frac{2}{\beta}}_{L^{2}(\Omega)}\bigr{]}\leq C\left[\|a^{*}\|^{2}_{X}+\frac{\sigma^{2}}{n\rho_{n}^{1+\frac{d}{4}/(1+\beta)}}\right]^{\frac{2}{\beta}}\frac{1}{n^{4(1+\beta)/d}}$$

$$+C\rho_{n}\left[\|a^{*}\|^{2}_{X}+\frac{\sigma^{2}}{n\rho_{n}^{1+\frac{d}{4}/(1+\beta)}}\right]^{\frac{1}{\beta}+1}$$

$$+C\left[\|a^{*}\|^{2}_{X}+\frac{\sigma^{2}}{n\rho_{n}^{1+\frac{d}{4}/(1+\beta)}}\right]^{\frac{1}{\beta}}\left(1+n^{-\frac{4}{d}}h^{-2}\right)e_{S}\|a_{n}^{*}\|^{2}_{X}.$$

Proof.

In view of the estimate (3.3) in Lemma 3.1, we see that

$$\|S_{\tau,h}a_{n}^{*}-Sa^{*}\|^{2}_{n}+n^{-\frac{4}{d}}\|S_{\tau,h}a_{n}^{*}-Sa^{*}\|^{2}_{H^{2}(\Omega)}\geq C\|S_{\tau,h}a_{n}^{*}-Sa^{*}\|^{2}_{L^{2}(\Omega)}.$$

For $\|S_{\tau,h}a_{n}^{*}-Sa^{*}\|^{2}_{L^{2}(\Omega)}$, one has

$$\|S_{\tau,h}a_{n}^{*}-Sa^{*}\|^{2}_{L^{2}(\Omega)} \geq C\|Sa_{n}^{*}-Sa^{*}\|^{2}_{L^{2}(\Omega)}-C\|S_{\tau,h}a_{n}^{*}-Sa_{n}^{*}\|^{2}_{L^{2}(\Omega)}$$

$$\geq C\|Sa_{n}^{*}-Sa^{*}\|^{2}_{L^{2}(\Omega)}-Ce_{S}\|a_{n}^{*}\|^{2}_{X}.$$

Moreover, from the second estimate in Remark 3.1 with $X=D((-\Delta)^{\beta})$, $\beta>0$, combined with Theorem 7 with $\gamma=0$ further implies that

$$\|S_{\tau,h}a_{n}^{*}-Sa^{*}\|^{2}_{n}+n^{-\frac{4}{d}}\|S_{\tau,h}a_{n}^{*}-Sa^{*}\|^{2}_{H^{2}(\Omega)}+Ce_{S}\|a_{n}^{*}\|^{2}_{X}\geq CM_{n}^{-\frac{2}{\beta}}\|a_{n}^{*}-a^{*}\|^{2+\frac{2}{\beta}}_{L^{2}(\Omega)}.$$

Here $M_{n}:=\|(-\Delta)^{\beta}(a_{n}^{*}-a^{*})\|_{L^{2}(\Omega)}$. On the other hand, Lemma 2.4 in [2] implies that

$$\|S_{\tau,h}a_{n}^{*}-Sa^{*}\|^{2}_{H^{2}(\Omega)} \leq C\|S_{\tau,h}a_{n}^{*}-Sa_{n}^{*}\|^{2}_{H^{2}(\Omega)}+\|Sa_{n}^{*}-Sa^{*}\|^{2}_{H^{2}(\Omega)}$$

$$\leq Ch^{-2}e_{S}\|a_{n}^{*}\|_{X}^{2}+C\|a_{n}^{*}-a^{*}\|^{2}_{L^{2}(\Omega)}.$$

Collecting all the above estimates and using the first inequality in Remark 3.1 with $X=D((-\Delta)^{\beta})$, we see that

$$\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2+\frac{2}{\beta}}_{L^{2}(\Omega)}\bigr{]}\leq C\mathbb{E}\bigl{[}M_{n}^{\frac{2}{\beta}}\bigr{]}\left[n^{-\frac{4}{d}}\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2}_{L^{2}(\Omega)}\bigr{]}+\rho_{n}\|a^{*}\|^{2}_{X}+\frac{\sigma^{2}}{n\rho_{n}^{\frac{d}{4}/(1+\beta)}}\right]$$

$$+C\mathbb{E}\bigl{[}M_{n}^{\frac{2}{\beta}}\bigr{]}\left(1+n^{-\frac{4}{d}}h^{-2}\right)e_{S}\|a_{n}^{*}\|^{2}_{X}.$$

Moreover, for $\varepsilon>0$, by the use of the Young inequality $|\xi\zeta|\leq\frac{C}{\varepsilon}|\xi|^{p}+\varepsilon|\zeta|^{q}$ with $p=\beta+1$ and $q=\frac{\beta+1}{\beta}$, we obtain

$$n^{-\frac{4}{d}}\|a_{n}-a^{*}\|^{2}_{L^{2}(\Omega)}\leq\frac{C}{\varepsilon}n^{-\frac{4(1+\beta)}{d}}+\varepsilon\|a_{n}-a^{*}\|^{2+\frac{2}{\beta}}_{L^{2}(\Omega)}.$$

Consequently, we see that

$$\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2+\frac{2}{\beta}}_{L^{2}(\Omega)}\bigr{]}$$

$$\leq C\mathbb{E}\bigl{[}M_{n}^{\frac{2}{\beta}}\bigr{]}\left[\frac{1}{\varepsilon}n^{-\frac{4(1+\beta)}{d}}+\varepsilon\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2+\frac{2}{\beta}}_{L^{2}(\Omega)}\bigr{]}+\rho_{n}\|a^{*}\|^{2}_{X}+\frac{\sigma^{2}}{n\rho_{n}^{\frac{d}{4}/(1+\beta)}}\right]$$

$$+C\mathbb{E}\bigl{[}M_{n}^{\frac{2}{\beta}}\bigr{]}\left(1+n^{-\frac{4}{d}}h^{-2}\right)e_{S}\|a_{n}^{*}\|^{2}_{X}.$$

By letting $\varepsilon=\frac{1}{2C\mathbb{E}\bigl{[}M_{n}^{\frac{2}{\beta}}\bigr{]}}$, we can see that the term $\varepsilon\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2+\frac{2}{\beta}}_{L^{2}(\Omega)}\bigr{]}$ on the right hand side of the above inequality can be absorbed, and then we get

$$\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2+\frac{2}{\beta}}_{L^{2}(\Omega)}\bigr{]}\leq C\mathbb{E}\bigl{[}M_{n}^{\frac{2}{\beta}}\bigr{]}\left[2C\mathbb{E}\bigl{[}M_{n}^{\frac{2}{\beta}}\bigr{]}n^{-\frac{4(1+\beta)}{d}}+\rho_{n}\|a^{*}\|^{2}_{X}+\frac{\sigma^{2}}{n\rho_{n}^{\frac{d}{4}/(1+\beta)}}\right]$$

$$+C\mathbb{E}\bigl{[}M_{n}^{\frac{2}{\beta}}\bigr{]}\left(1+n^{-\frac{4}{d}}h^{-2}\right)e_{S}\|a_{n}^{*}\|^{2}_{X}.$$

Now by using the second estimate in Remark 3.1 with $X=D((-\Delta)^{\beta})$, and noting the definition $M_{n}:=\|(-\Delta)^{\beta}(a_{n}^{*}-a^{*})\|_{L^{2}(\Omega)}$, it follows that

$$\mathbb{E}\bigl{[}M_{n}^{\frac{2}{\beta}}\bigr{]}\leq C\left[\|a^{*}\|^{2}_{X}+\frac{\sigma^{2}}{n\rho_{n}^{1+\frac{d}{4}/(1+\beta)}}\right]^{\frac{1}{\beta}},$$

from which we further know that

$$\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2+\frac{2}{\beta}}_{L^{2}(\Omega)}\bigr{]}\leq C\mathbb{E}\bigl{[}M^{\frac{2}{\beta}}\bigr{]}\left[\frac{\mathbb{E}\bigl{[}M^{\frac{2}{\beta}}\bigr{]}}{n^{4(1+\beta)/d}}+\rho_{n}\|a^{*}\|^{2}_{X}+\frac{\sigma^{2}}{n\rho_{n}^{\frac{d}{4}/(1+\beta)}}\right]$$

$$+C\mathbb{E}\bigl{[}M_{n}^{\frac{2}{\beta}}\bigr{]}\left(1+n^{-\frac{4}{d}}h^{-2}\right)e_{S}\|a_{n}^{*}\|^{2}_{X}$$

$$\leq C\left[\|a^{*}\|^{2}_{X}+\frac{\sigma^{2}}{n\rho_{n}^{1+\frac{d}{4}/(1+\beta)}}\right]^{\frac{2}{\beta}}\frac{1}{n^{4(1+\beta)/d}}$$

$$+C\rho_{n}\left[\|a^{*}\|^{2}_{X}+\frac{\sigma^{2}}{n\rho_{n}^{1+\frac{d}{4}/(1+\beta)}}\right]^{\frac{1}{\beta}+1}$$

$$+C\left[\|a^{*}\|^{2}_{X}+\frac{\sigma^{2}}{n\rho_{n}^{1+\frac{d}{4}/(1+\beta)}}\right]^{\frac{1}{\beta}}\left(1+n^{-\frac{4}{d}}h^{-2}\right)e_{S}\|a_{n}^{*}\|^{2}_{X}.$$

We can complete the proof of the theorem. ∎

Remark 3.2.

Remark 3.1 and Theorem 9 indicate that the optimal regularization parameter has the form: $\rho_{n}^{\frac{1}{2}+\frac{d}{8}/(1+\beta)}=O(\sigma n^{-\frac{1}{2}}\|a^{*}\|^{-1}_{X})$ for $\beta\in(\frac{d}{4},1]$, and the optimal error will be

$$\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2+\frac{2}{\beta}}_{L^{2}(\Omega)}\bigr{]}\leq C\left(\rho_{n}+\frac{1}{n^{4(1+\beta)/d}}\right)+C\left(1+n^{-\frac{4}{d}}h^{-2}\right)e_{S}\|a_{n}^{*}\|^{2}_{X}.$$

If $\left(1+n^{-\frac{4}{d}}h^{-2}\right)e_{S}\|a_{n}^{*}\|^{2}_{X}<\rho_{n}+\frac{1}{n^{4(1+\beta)/d}}$, we have

$$\mathbb{E}\bigl{[}\|a_{n}^{*}-a^{*}\|^{2+\frac{2}{\beta}}_{L^{2}(\Omega)}\bigr{]}\leq C\left(\rho_{n}+\frac{1}{n^{4(1+\beta)/d}}\right).$$

4 Regularity theory for model (1.1)

In this section, we recall the well-posedness result of the initial-boundary value problem (1.1). For this, we make several settings. Let $L^{2}(\Omega)$ be the square-integrable function space with inner product $(\cdot,\cdot)_{L^{2}(\Omega)}$ (or $(\cdot,\cdot)$ for short) and let $H^{1}(\Omega)$, $H^{2}(\Omega)$ etc. be the usual Sobolev spaces.

The set $\{\lambda_{k},\varphi_{k}\}_{k=1}^{\infty}$ constitutes the Dirichlet eigensystem of the elliptic operator $-\Delta:H^{2}(\Omega)\cap H_{0}^{1}(\Omega)\to L^{2}(\Omega)$, specifically,

$$\begin{cases}-\Delta\varphi_{k}=\lambda_{k}\varphi_{k}&\text{in }\Omega,\\ \varphi_{k}=0&\text{on }\partial\Omega,\end{cases}$$

where $\lambda_{k}$ is the eigenvalue of the operator $-\Delta$ and satisfies $0<\lambda_{1}\leq\lambda_{2}\leq\ldots,\lambda_{k}\rightarrow\infty$ as $k\rightarrow\infty$, and $\varphi_{k}$ is the eigenfunctions corresponding to the value $\lambda_{k}$ and $\{\varphi_{k}\}_{k=1}^{\infty}$ forms an orthonormal basis in $L^{2}(\Omega)$. We have the asymptotic behavior of the eigenvalue $\lambda_{k}\sim k^{2/d}$ as $k\to\infty$. Then for $\gamma\in\mathbb{R}$, fractional power $(-\Delta)^{\gamma}$ can be defined

$$(-\Delta)^{\gamma}\psi:=\sum_{k=1}^{\infty}\lambda_{k}^{\gamma}(\psi,\varphi_{k})\varphi_{k},\quad\psi\in D((-\Delta)^{\gamma}),$$

where

$$\mathcal{D}((-\Delta)^{\gamma}):=\left\{\psi\in L^{2}(\Omega);\sum_{k=1}^{\infty}\left|\lambda_{k}^{\gamma}(\psi,\varphi_{k})\right|^{2}<\infty\right\}.$$

The space $D((-\Delta)^{\gamma})$ is a Hilbert space equipped with the inner product

$$(\psi,\phi)_{D((-\Delta)^{\gamma})}=\left((-\Delta)^{\gamma}\psi,(-\Delta)^{\gamma}\phi\right)_{L^{2}(\Omega)}.$$

Moreover, we define the norm

$$\left\|\psi\right\|_{\mathcal{D}((-\Delta)^{\gamma})} =\left((-\Delta)^{\gamma}\psi,(-\Delta)^{\gamma}\psi\right)_{L^{2}(\Omega)}^{\frac{1}{2}}=\left(\sum_{n=1}^{\infty}\left|\lambda_{n}^{\gamma}(\psi,\varphi_{n})\right|^{2}\right)^{\frac{1}{2}}.$$

For short, we also denote the inner product $(\cdot,\cdot)_{\mathcal{D}((-\Delta)^{\gamma})}$ and the norm $\|\cdot\|_{\mathcal{D}((-\Delta)^{\gamma})}$ as $(\cdot,\cdot)_{\gamma}$ and $\|\cdot\|_{\gamma}$ if no conflict occurs. Furthermore, it satisfies $\mathcal{D}((-\Delta)^{\gamma})\subset{H^{2\gamma}(\Omega)}$ for $\gamma>0$. In particular, we have $\mathcal{D}((-\Delta)^{\frac{1}{2}})=H_{0}^{1}(\Omega)$, $\mathcal{D}((-\Delta)^{-\frac{1}{2}})=H^{-1}(\Omega)$ and the norm equivalence $\left\|\cdot\right\|_{\mathcal{D}((-\Delta)^{\gamma})}\sim\left\|\cdot\right\|_{H^{2\gamma}(\Omega)}$ with $\gamma=\pm\frac{1}{2}$.

The regularity of the solution is based on the boundedness of the Mittag-Leffler functions in Lemma 4.1.

Lemma 4.1.

([14]) If $0<\alpha<2$, $\beta$ is an arbitrary complex number and $\mu$ is an arbitrary real number such that

$$\frac{\pi\alpha}{2}<\mu<\min\{\pi,\pi\alpha\},$$

then

$$|E_{\alpha,\beta}(z)|\leq\frac{C}{1+|z|},\quad\mu\leq|\arg z|\leq\pi,$$

where $E_{\alpha,\beta}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma(\alpha k+\beta)},~{}z\in\mathbb{C}$.

Theorem 10.

If $a_{0}=0$ and $a_{1}\in D((-\Delta)^{\gamma})$, $\gamma\in(\frac{d}{4},1]$, then

$$|\partial_{t}u|\leq C\|(-\Delta)^{\gamma}a_{1}\|_{L^{2}}.$$

Proof.

The solution to the problem (1.1) can be expressed as:

$$u(x,t)=\sum_{n=1}^{\infty}tE_{\alpha,2}(-\lambda_{n}t^{\alpha})(a_{1},\varphi_{n})\varphi_{n}(x).$$ (4.1)

From the definition of the Mittag-Leffler function, Lemma 4.1 and (4.1), one has

$$|\partial_{t}u(x,t)| =\bigg{|}\sum_{n=1}^{\infty}E_{\alpha,1}(-\lambda_{n}t^{\alpha})(a_{1},\varphi_{n})\varphi_{n}(x)\bigg{|}$$

$$\leq C\sum_{n=1}^{\infty}\big{|}E_{\alpha,1}(-\lambda_{n}t^{\alpha})\big{|}\big{|}(a_{1},\varphi_{n})\big{|}$$

$$\leq C\sum_{n=1}^{\infty}\lambda_{n}^{-\gamma}\big{|}E_{\alpha,1}(-\lambda_{n}t^{\alpha})\big{|}\lambda_{n}^{\gamma}\big{|}(a_{1},\varphi_{n})\big{|}$$

$$\leq C\sqrt{\sum_{n=1}^{\infty}\lambda_{n}^{-2\gamma}\big{|}E_{\alpha,1}(-\lambda_{n}t^{\alpha})\big{|}^{2}}\sqrt{\sum_{n=1}^{\infty}\lambda_{n}^{2\gamma}\big{|}(a_{1},\varphi_{n})\big{|}^{2}}$$

$$\leq C\|(-\Delta)^{\gamma}a_{1}\|_{L^{2}},$$

where the last inequality follows from the inequality

$$\sum_{n=1}^{\infty}\lambda_{n}^{-2\gamma}\big{|}E_{\alpha,1}(-\lambda_{n}t^{\alpha})\big{|}^{2}\leq\sum_{n=1}^{\infty}\frac{\lambda_{n}^{-2\gamma}}{(1+\lambda_{n}t^{\alpha})^{2}}\leq\sum_{n=1}^{\infty}\lambda_{n}^{-2\gamma}\leq C,$$

where $\lambda_{n}\sim n^{\frac{2}{d}}$ and $\gamma\in(\frac{d}{4},1]$. This finishes the proof of the theorem. ∎

Theorem 11.

If $a_{0}=0$ and $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{d}{4},1]$ and $0<\epsilon\ll 1$, then

$$|\partial_{tt}u|\leq Ct^{\epsilon\alpha-1}\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}}.$$

Proof.

Similar to Theorem 10, one has

$$|\partial_{tt}u(x,t)| \leq C\sum_{n=1}^{\infty}\big{|}\lambda_{n}^{1-\epsilon}t^{\alpha-1}E_{\alpha,\alpha}(-\lambda_{n}t^{\alpha})\big{|}\big{|}\lambda_{n}^{\epsilon}(a_{1},\varphi_{n})\big{|}$$

$$\leq Ct^{\epsilon\alpha-1}\sum_{n=1}^{\infty}\big{|}(\lambda_{n}t^{\alpha})^{1-\epsilon}E_{\alpha,\alpha}(-\lambda_{n}t^{\alpha})\big{|}\big{|}\lambda_{n}^{\epsilon}(a_{1},\varphi_{n})\big{|}$$

$$\leq Ct^{\epsilon\alpha-1}\sup_{n}\frac{(\lambda_{n}t^{\alpha})^{1-\epsilon}}{1+\lambda_{n}t^{\alpha}}\sum_{n=1}^{\infty}\big{|}\lambda_{n}^{\epsilon}(a_{1},\varphi_{n})\big{|},$$

where the last inequality follows from the inequality $\frac{(\lambda_{n}t^{\alpha})^{1-\epsilon}}{1+\lambda_{n}t^{\alpha}}<C$, $n\geq 1$. Then, it gives that

$$|\partial_{tt}u(x,t)| \leq Ct^{\epsilon\alpha-1}\sum_{n=1}^{\infty}\big{|}\lambda_{n}^{\epsilon}(a_{1},\varphi_{n})\big{|}$$

$$\leq Ct^{\epsilon\alpha-1}\sum_{n=1}^{\infty}\lambda_{n}^{-\gamma}\lambda_{n}^{\epsilon+\gamma}\big{|}(a_{1},\varphi_{n})\big{|}$$

$$\leq Ct^{\epsilon\alpha-1}\sqrt{\sum_{n=1}^{\infty}\lambda_{n}^{-2\gamma}}\sqrt{\sum_{n=1}^{\infty}\lambda_{n}^{2(\epsilon+\gamma)}\big{|}(a_{1},\varphi_{n})\big{|}^{2}}$$

$$\leq Ct^{\epsilon\alpha-1}\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}},$$

where $\sqrt{\sum_{n=1}^{\infty}\lambda_{n}^{-2\gamma}}<C$, $\lambda_{n}\sim n^{\frac{2}{d}}$ and $\gamma\in(\frac{d}{4},1]$. We complete the proof of the theorem. ∎

Lemma 4.2.

If $a_{0}=0$ and $a_{1}\in L^{2}(\Omega)$, then

$$\|\partial_{t}^{m}u\|_{L^{2}}\leq Ct^{1-m}\|a_{1}\|_{L^{2}},~{}~{}m=0,1.$$

Proof.

From (4.1), indicates that

$$\|u\|_{L^{2}}^{2} =\bigg{\|}\sum_{n=1}^{\infty}tE_{\alpha,2}(-\lambda_{n}t^{\alpha})(a_{1},\varphi_{n})\varphi_{n}(x)\bigg{\|}_{L^{2}}^{2}$$

$$=\sum_{n=1}^{\infty}t^{2}E_{\alpha,2}(-\lambda_{n}t^{\alpha})^{2}|(a_{1},\varphi_{n})|^{2}$$

$$\leq\sup_{n}t^{2}E_{\alpha,2}(-\lambda_{n}t^{\alpha})^{2}\sum_{n=1}^{\infty}|(a_{1},\varphi_{n})|^{2}$$

$$\leq\sup_{n}\frac{Ct^{2}}{(1+\lambda_{n}t^{\alpha})^{2}}\|a_{1}\|_{L^{2}}^{2}\leq Ct^{2}\|a_{1}\|_{L^{2}}^{2}.$$

Similarly, we get

$$\|\partial_{t}u\|_{L^{2}}^{2} =\bigg{\|}\sum_{n=1}^{\infty}E_{\alpha,1}(-\lambda_{n}t^{\alpha})(a_{1},\varphi_{n})\varphi_{n}(x)\bigg{\|}_{L^{2}}^{2}$$

$$=\sum_{n=1}^{\infty}E_{\alpha,1}(-\lambda_{n}t^{\alpha})^{2}|(a_{1},\varphi_{n})|^{2}$$

$$\leq\sup_{n}E_{\alpha,1}(-\lambda_{n}t^{\alpha})^{2}\sum_{n=1}^{\infty}|(a_{1},\varphi_{n})|^{2}$$

$$\leq\sup_{n}\frac{C}{(1+\lambda_{n}t^{\alpha})^{2}}\|a_{1}\|_{L^{2}}^{2}\leq C\|a_{1}\|_{L^{2}}^{2}.$$

The proof of the lemma is complete. ∎

Lemma 4.3.

If $a_{0}=0$ and $a_{1}\in D((-\Delta)^{\gamma})$, $\gamma\in(0,1]$, then

$$\|\partial_{t}^{m}u\|_{L^{2}}\leq Ct^{\gamma\alpha-m+1}\|(-\Delta)^{\gamma}a_{1}\|_{L^{2}},~{}~{}m=2,3.$$

Proof.

From (4.1), indicates that

$$\|\partial_{tt}u\|_{L^{2}}^{2} =\bigg{\|}\sum_{n=1}^{\infty}\lambda_{n}t^{\alpha-1}E_{\alpha,\alpha}(-\lambda_{n}t^{\alpha})(a_{1},\varphi_{n})\varphi_{n}(x)\bigg{\|}_{L^{2}}^{2}$$

$$=\sum_{n=1}^{\infty}\lambda_{n}^{2}t^{2\alpha-2}E_{\alpha,\alpha}(-\lambda_{n}t^{\alpha})^{2}|(a_{1},\varphi_{n})|^{2}$$

$$=\sum_{n=1}^{\infty}(\lambda_{n}t^{\alpha})^{2-2\gamma}t^{2\gamma\alpha-2}E_{\alpha,\alpha}(-\lambda_{n}t^{\alpha})^{2}\lambda_{n}^{2\gamma}|(a_{1},\varphi_{n})|^{2}$$

$$\leq Ct^{2\gamma\alpha-2}\sup_{n}(\lambda_{n}t^{\alpha})^{2-2\gamma}E_{\alpha,\alpha}(-\lambda_{n}t^{\alpha})^{2}\sum_{n=1}^{\infty}\lambda_{n}^{2\gamma}|(a_{1},\varphi_{n})|^{2}$$

$$\leq Ct^{2\gamma\alpha-2}\sup_{n}\frac{(\lambda_{n}t^{\alpha})^{2-2\gamma}}{(1+\lambda_{n}t^{\alpha})^{2}}\sum_{n=1}^{\infty}\lambda_{n}^{2\gamma}|(a_{1},\varphi_{n})|^{2}$$

$$\leq Ct^{2\gamma\alpha-2}\|(-\Delta)^{\gamma}a_{1}\|_{L^{2}}^{2},$$

where the last inequality holds since $\sup_{n}\frac{(\lambda_{n}t^{\alpha})^{2-2\gamma}}{(1+\lambda_{n}t^{\alpha})^{2}}<C$. Similarly, we get

$$\|\partial_{ttt}u\|_{L^{2}}^{2} =\bigg{\|}\sum_{n=1}^{\infty}\lambda_{n}t^{\alpha-2}E_{\alpha,\alpha-1}(-\lambda_{n}t^{\alpha})(a_{1},\varphi_{n})\varphi_{n}(x)\bigg{\|}_{L^{2}}^{2}$$

$$=\sum_{n=1}^{\infty}(\lambda_{n}t^{\alpha})^{2-2\gamma}t^{2\gamma\alpha-4}E_{\alpha,\alpha-1}(-\lambda_{n}t^{\alpha})^{2}\lambda_{n}^{2\gamma}|(a_{1},\varphi_{n})|^{2}$$

$$\leq Ct^{2\gamma\alpha-4}\sup_{n}(\lambda_{n}t^{\alpha})^{2-2\gamma}E_{\alpha,\alpha-1}(-\lambda_{n}t^{\alpha})^{2}\sum_{n=1}^{\infty}\lambda_{n}^{2\gamma}|(a_{1},\varphi_{n})|^{2}$$

$$\leq Ct^{2\gamma\alpha-4}\|(-\Delta)^{\gamma}a_{1}\|_{L^{2}}^{2}.$$

Collecting all the above estimates, we finish the proof of the lemma. ∎

Theorem 12.

If $a_{0}=0$ and $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in\mathbb{R}$ and $\epsilon\in(0,1]$, then

$$\|\partial_{t}^{m}(-\Delta)^{\gamma+\epsilon}u\|_{L^{2}}\leq Ct^{1-m}\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}},~{}~{}m=0,1,$$

$$\|\partial_{t}^{m}(-\Delta)^{\gamma}u\|_{L^{2}}\leq Ct^{\epsilon\alpha-m+1}\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}},~{}~{}m=2,3.$$

Proof.

The desired results are directly obtained following the idea in lemmas 4.2 and 4.3. ∎

Theorem 13.

If $a_{0}=0$ and $a_{1}\in D((-\Delta)^{\gamma+\epsilon})$, $\gamma\in(\frac{1}{4},1]$ and $0<\epsilon\ll 1$, then

$$\|\partial_{t}^{m}(-\Delta)^{\frac{1}{2}}u\|_{L^{2}}\leq Ct^{1-m-\frac{\alpha}{4}}\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}},~{}~{}m=0,1,2,3.$$

Proof.

For $\xi\in(0,1]$, it gives that

$$\|(-\Delta)^{\frac{1}{2}}u\|_{L^{2}}^{2} =\bigg{\|}\sum_{n=1}^{\infty}\lambda_{n}^{\frac{1}{2}}tE_{\alpha,2}(-\lambda_{n}t^{\alpha})(a_{1},\varphi_{n})\varphi_{n}(x)\bigg{\|}_{L^{2}}^{2}$$

$$=\sum_{n=1}^{\infty}\lambda_{n}t^{2}E_{\alpha,2}(-\lambda_{n}t^{\alpha})^{2}|(a_{1},\varphi_{n})|^{2}$$

$$=\sum_{n=1}^{\infty}t^{\alpha(2\xi-2)+2}(\lambda_{n}t^{\alpha})^{2-2\xi}E_{\alpha,2}(-\lambda_{n}t^{\alpha})^{2}\lambda_{n}^{2\xi-1}|(a_{1},\varphi_{n})|^{2}$$

$$\leq\sup_{n}\frac{(\lambda_{n}t^{\alpha})^{2-2\xi}}{(1+\lambda_{n}t^{\alpha})^{2}}t^{\alpha(2\xi-2)+2}\sum_{n=1}^{\infty}\lambda_{n}^{2\xi-1}|(a_{1},\varphi_{n})|^{2}$$

$$\leq Ct^{\alpha(2\xi-2)+2}\|(-\Delta)^{\xi-\frac{1}{2}}a_{1}\|_{L^{2}}^{2},$$

where the last inequality holds since $\sup_{n}\frac{(\lambda_{n}t^{\alpha})^{2-2\xi}}{(1+\lambda_{n}t^{\alpha})^{2}}<C$. We take $\xi=\frac{3}{4}\leq\gamma+\epsilon+\frac{1}{2}$, it arrives that

$$\|(-\Delta)^{\frac{1}{2}}u\|_{L^{2}} \leq Ct^{1-\frac{\alpha}{4}}\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}}.$$

Similarly, one has $\|\partial_{t}(-\Delta)^{\frac{1}{2}}u\|_{L^{2}}\leq Ct^{-\frac{\alpha}{4}}\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}}$. The proof is complete. ∎

Some properties of $E_{\alpha,\beta}(-\lambda_{n}t^{\alpha})$, $\alpha>0$, $\beta\in\mathbb{R}$, are need to deduce the estimates of $\|\partial_{t}^{m}(-\Delta)^{\frac{1}{2}}u\|_{L^{2}}$, $m=2,3$,

$$\partial_{t}[tE_{\alpha,1}(-\lambda_{n}t^{\alpha})] =\partial_{t}\bigg{[}t\sum_{k=0}^{\infty}\frac{(-\lambda_{n}t^{\alpha})^{k}}{\Gamma{(k\alpha+1)}}\bigg{]}$$

$$=\partial_{t}\bigg{[}t+\sum_{k=1}^{\infty}\frac{(-\lambda_{n})^{k}t^{k\alpha+1}}{\Gamma{(k\alpha+1)}}\bigg{]}$$

$$=1+\sum_{k=1}^{\infty}\frac{(-\lambda_{n}t^{\alpha})^{k}}{\Gamma{(k\alpha+1)}}+\sum_{k=1}^{\infty}\frac{(-\lambda_{n}t^{\alpha})^{k}}{\Gamma{(k\alpha)}}$$

$$:=E_{\alpha,1}(-\lambda_{n}t^{\alpha})+E_{\alpha,0}(-\lambda_{n}t^{\alpha}),$$

and, we get

$$\partial_{t}[E_{\alpha,1}(-\lambda_{n}t^{\alpha})] =\partial_{t}[t^{-1}tE_{\alpha,1}(-\lambda_{n}t^{\alpha})]$$

$$=-t^{-1}E_{\alpha,1}(-\lambda_{n}t^{\alpha})+t^{-1}\partial_{t}[tE_{\alpha,1}(-\lambda_{n}t^{\alpha})]$$

$$=-t^{-1}E_{\alpha,1}(-\lambda_{n}t^{\alpha})+t^{-1}\big{(}E_{\alpha,1}(-\lambda_{n}t^{\alpha})+E_{\alpha,0}(-\lambda_{n}t^{\alpha})\big{)}$$

$$=t^{-1}E_{\alpha,0}(-\lambda_{n}t^{\alpha}).$$

Based on above results, it holds that

$$\|\partial_{t}^{2}(-\Delta)^{\frac{1}{2}}u\|_{L^{2}}^{2} =\bigg{\|}\sum_{n=1}^{\infty}\lambda_{n}^{\frac{1}{2}}t^{-1}E_{\alpha,0}(-\lambda_{n}t^{\alpha})(a_{1},\varphi_{n})\varphi_{n}(x)\bigg{\|}_{L^{2}}^{2}$$

$$=\sum_{n=1}^{\infty}\lambda_{n}t^{-2}E_{\alpha,0}(-\lambda_{n}t^{\alpha})^{2}|(a_{1},\varphi_{n})|^{2}$$

$$=\sum_{n=1}^{\infty}t^{\alpha(2\xi-2)-2}(\lambda_{n}t^{\alpha})^{2-2\xi}E_{\alpha,0}(-\lambda_{n}t^{\alpha})^{2}\lambda_{n}^{2\xi-1}|(a_{1},\varphi_{n})|^{2}$$

$$\leq\sup_{n}\frac{(\lambda_{n}t^{\alpha})^{2-2\xi}}{(1+\lambda_{n}t^{\alpha})^{2}}t^{\alpha(2\xi-2)-2}\sum_{n=1}^{\infty}\lambda_{n}^{2\xi-1}|(a_{1},\varphi_{n})|^{2}$$

$$\leq Ct^{\alpha(2\xi-2)-2}\|(-\Delta)^{\xi-\frac{1}{2}}a_{1}\|_{L^{2}}^{2},$$

where $\xi\in(0,1]$, taking $\xi=\frac{3}{4}\leq\gamma+\epsilon+\frac{1}{2}$, i. e.

$$\|\partial_{t}^{2}(-\Delta)^{\frac{1}{2}}u\|_{L^{2}} \leq Ct^{-\frac{\alpha}{4}-1}\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}}^{2}.$$

Similarly, we get

$$\|\partial_{t}^{3}(-\Delta)^{\frac{1}{2}}u\|_{L^{2}} \leq Ct^{-\frac{\alpha}{4}-2}\|(-\Delta)^{\gamma+\epsilon}a_{1}\|_{L^{2}}^{2}.$$

Theorem 14.

If $a_{0}=0$ and $a_{1}\in L^{2}(\Omega)$, then

$$\|\partial_{t}^{m}v\|_{L^{2}}\leq Ct^{1-\nu-m}\|a_{1}\|_{L^{2}},~{}~{}m=0,1,2,3.$$

Proof.

From (4.1), indicates that

$$\|v\|_{L^{2}}^{2} =\bigg{\|}\frac{1}{\Gamma{(1-\nu)}}\int_{0}^{t}(t-s)^{-\nu}\partial_{s}u(s)ds\bigg{\|}_{L^{2}}^{2}$$

$$=\bigg{\|}\frac{1}{\Gamma{(1-\nu)}}\sum_{n=1}^{\infty}\int_{0}^{t}(t-s)^{-\nu}E_{\alpha,1}(-\lambda_{n}s^{\alpha})ds(a_{1},\varphi_{n})\varphi_{n}(x)\bigg{\|}_{L^{2}}^{2}$$

$$=\sum_{n=1}^{\infty}\big{(}t^{1-\nu}E_{\alpha,2-\nu}(-\lambda_{n}t^{\alpha})\big{)}^{2}|(a_{1},\varphi_{n})|^{2}$$

$$\leq t^{2-2\nu}\sup_{n}E_{\alpha,2-\nu}(-\lambda_{n}t^{\alpha})^{2}\sum_{n=1}^{\infty}|(a_{1},\varphi_{n})|^{2}$$

$$\leq Ct^{2-2\nu}\|a_{1}\|_{L^{2}}^{2}.$$

Similarly, we get

$$\|\partial_{t}^{m}v\|_{L^{2}}^{2} =\sum_{n=1}^{\infty}\big{[}\partial_{t}^{m}\big{(}t^{1-\nu}E_{\alpha,2-\nu}(-\lambda_{n}t^{\alpha})\big{)}\big{]}^{2}|(a_{1},\varphi_{n})|^{2}$$

$$=t^{2-2\nu-2m}\sum_{n=1}^{\infty}\big{(}E_{\alpha,2-\nu-m}(-\lambda_{n}t^{\alpha})\big{)}^{2}|(a_{1},\varphi_{n})|^{2}$$

$$\leq Ct^{2-2\nu-2m}\|a_{1}\|_{L^{2}}^{2}.$$

We complete the proof of the theorem. ∎

5 Numerical experiment

In this section, we carry out some numerical experiments to check the theoretical results. In Example 1, we consider testing the convergence statements for one dimensional cases. Let us return to the literature where the SFOR method is proposed [12]. Remark 2.2 therein explains why auxiliary variables were introduced to extract the singular term $a_{1}(x)\omega_{2-\alpha}(t)$. Based on this, they adopted the following numerical framework:

$$\begin{cases}\partial_{t}^{\nu}\textbf{v}-\Delta\textbf{u}=t\Delta a_{1}(x),&(x,t)\in\Omega\times(0,T),\\ \textbf{v}=\partial_{t}^{\nu}\textbf{u},&(x,t)\in\Omega\times(0,T),\\ \textbf{u}(x,0)=\textbf{v}(x,0)=0,&x\in\Omega,\\ \textbf{u}(x,t)=\textbf{v}(x,t)=0,&(x,t)\in\partial\Omega\times(0,T),\end{cases}$$ (5.1)

where $u=\textbf{u}+ta_{1}(x)$. The limitation of (5.1) is $(\Delta a_{1}(x),\phi(x))$, $x\in\Omega$ exists, where $\phi(x)$ is basis function from finite element space. Our numerical framework (2.1) relaxes this requirement. Furthermore, we find that the optimal convergence is reached, despite the presence of $a_{1}(x)\omega_{2-\alpha}(t)$, when the mesh parameter $r=\frac{4-\alpha}{2-\alpha}$, $1<\alpha<2-\epsilon$, $\epsilon$ is a fixed positive constant.

Theoretically, it could not work for the case $\alpha\rightarrow 2^{-}$ since $r=\frac{4-\alpha}{2-\alpha}\rightarrow+\infty$. For $\Delta a_{1}\in L^{2}(\Omega)$, scheme (5.1) is helpful. The reason is presented in Remark 5.1. For more general $a_{1}(x)$, we can use scheme (2.1) with bounded $r$. It gives detailed instructions on how to apply the SFOR method in different application cases.

Remark 5.1.

The point is v become more regular. Following the idea of Theorem 14, one has $\|\partial_{t}^{m}\textbf{v}\|_{L^{2}}\leq Ct^{1+\nu-m}\|a_{1}\|_{L^{2}}$, $m=0,1,2,3$. Combining the analysis in Theorems 4 and 5, it implies that the mesh parameter $r=2$ is enough.

Example 1.

(One dimensional) Forward problems (1.1) with $\Omega=(0,\pi)$, $T=0.1$,

$$(a)~{}~{}a_{1}(x) =\sin(x),$$

$$(b)~{}~{}a_{1}(x) =x,~{}~{}x\in(0,\pi/2],~{}~{}a_{1}(x)=\pi-x,~{}~{}x\in(\pi/2,\pi).$$

The size of the space grids $h=\frac{\pi}{200}$, $N$ is the number of partitions in the time grids. $e_{L^{2}}=\max_{1\leq n\leq N}\|u_{h}^{n}-U_{h}^{n}\|_{L^{2}}$, where $u_{h}^{n}$ and $U_{h}^{n}$ are the reference solution ($h=\frac{\pi}{200},N=2048$) and the numerical solution, respectively. Furthermore, to test the convergence rate, let $Order=\log_{2}(e_{L^{2}}(N/2)/e_{L^{2}}(N))$.

Example 2.

(One dimensional) Backward problems (1.1) with $\Omega=(0,\pi)$, $T=0.1$,

$$(a)~{}~{}a_{1}(x) =\sin(x),~{}~{}\alpha=1.5,$$

$$(b)~{}~{}a_{1}(x) =x,~{}~{}x\in(0,\pi/2],~{}~{}a_{1}(x)=\pi-x,~{}~{}x\in(\pi/2,\pi),~{}~{}\alpha=1.5.$$

We use the numerical framework (2.3) with the mesh parameter $r=\frac{4-\alpha}{2-\alpha}$, where the size of the space grids $h=\frac{\pi}{20}$, $N=2048$ is the number of partitions in the time grids. Reconstruction $a_{1}^{rec}$ is found by gradient descent method with initial guess $a_{1}^{rec}=0$.

Table 3: Scheme (2.1) for Example 1 (a) with $\alpha=1.25$.

$N$	$r=1$		$r_{opt}=\frac{4-\alpha}{2-\alpha}$		$r=\frac{4-\alpha}{\alpha}$
	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order
16	9.7511e-03	-	3.2933e-03	-	3.5919e-03	-
32	5.8969e-03	0.726	1.3099e-03	1.330	1.5200e-03	1.241
64	3.4701e-03	0.765	5.0947e-04	1.362	6.2187e-04	1.289
128	1.9768e-03	0.812	1.9505e-04	1.385	2.4760e-04	1.329
Optimal order	1.375

In example 2, we consider recovering the initial function $a_{1}$. For the case $(a)$, let $\sigma=0.05$, then the noise lever $\frac{\sigma}{\|Sa_{1}\|_{L^{\infty}}}\approx 50\%$. The number of observation points $n$ are taken $11$, $49$ and $199$ in numerical tests. The optimal regularization parameter of $H^{1}$ regularization $\rho_{n}=O((\sigma n^{-\frac{1}{2}}\|a^{*}\|^{-1}_{H^{1}})^{12/7})=O(5.1\times 10^{-4})$, $O(1.4\times 10^{-4})$ and $O(4.3\times 10^{-5})$. In Figure 3, 11 observation points are presented. Then, $L^{2}$ error $\frac{\|a_{1}^{rec}-a_{1}\|_{L^{2}}}{\|a_{1}\|_{L^{2}}}$ is given when $\rho_{11}=10^{-k}$. Furthermore, we compare $a_{1}$ with the optimal $a_{1}^{rec}$ when $\rho_{11}=10^{-3.25}\approx 5.6\times 10^{-4}$. It shows that there is a big gap between reconstruction $a_{1}^{rec}$ and $a_{1}$. To obtain a better $a_{1}^{rec}$, increasing observational data is a viable strategy. The related results in Figures 6-9 support the conclusion of Remark 3.2.

For the case $(b)$, let $\sigma=0.015$, then the noise lever $\frac{\sigma}{\|Sa_{1}\|_{L^{\infty}}}\approx 10\%$. The number of observation points $n$ is taken $199$ in numerical tests. The optimal regularization parameter $\rho_{199}=O((\sigma n^{-\frac{1}{2}}\|a^{*}\|^{-1}_{H^{1}})^{12/7})=O(3.0\times 10^{-6})$. $a_{1}^{rec}$ found by $H^{1}$ regularization with $\rho_{199}=3\times 10^{-6}$ is presented in Figure 11. It verifies our method also work for $a_{1}$ with singularity point.

Example 3.

(Two dimensional) Forward problems (1.1) with $\Omega=(0,1)\times(0,1)$, $T=0.1$,

$$a_{1}(x,y) =\ln(1+10x)(x-1)\sin^{\frac{3}{4}}(\pi y).$$

The size of the space grids $h=\frac{1}{30}$, $N$ is the number of partitions in the time grids. Let the numerical solution on fine mesh be the reference solution ($h=\frac{1}{30}$,$N=1280$).

In Example 3, when we take $r=\frac{4-\alpha}{2-\alpha}$, the convergence rate reaches optimal $2-\alpha/2$. Under the same accuracy requirements, the use of optimal mesh parameters can save the cost of calculation. Furthermore, it improves the efficiency of numerical inversion.

Example 4.

(Two dimensional) Backward problems (1.1) with $\Omega=(0,1)\times(0,1)$, $T=0.1$,

$$a_{1}(x,y) =\ln(1+10x)(x-1)\sin^{\frac{3}{4}}(\pi y),~{}~{}\alpha=1.25.$$

We use the numerical framework (2.3) with the mesh parameter $r=\frac{4-\alpha}{2-\alpha}$, where the size of the space grids $h=\frac{1}{30}$, $N=160$ is the number of partitions in the time grids. Reconstruction $a_{1}^{rec}$ is found by gradient descent method with initial guess $a_{1}^{rec}=0$.

In Example 4, let $\sigma=0.01$, then the noise lever $\frac{\sigma}{\|Sa_{1}\|_{L^{\infty}}}\approx 12\%$. We take the points of space grids as measurement points. The optimal regularization parameter of $H^{1}$ regularization $\rho_{841}=O((\sigma n^{-\frac{1}{2}}\|a^{*}\|^{-1}_{H^{1}})^{3/2})=O(1.5\times 10^{-6})$. The reconstruction results are presented with different regularization parameters. $L^{2}$ errors $\frac{\|a_{1}^{rec}-a_{1}\|_{L^{2}}}{\|a_{1}\|_{L^{2}}}=15.0\%$, $11.0\%$ and $24.5\%$ when $\rho_{841}=10^{-5}$, $2\times 10^{-6}$ and $10^{-7}$, respectively. $a_{1}^{rec}$ with $\rho_{841}=2\times 10^{-6}$ close to the exact one, see Figures 13 and 15. However, the peak in Figure 13 is not consistent with that of the exact one. And in Figure 15, the reconstruction becomes blurry when $\rho_{841}$ is too small.

Table 4: Scheme (2.1) for Example 1 (a) with $\alpha=1.5$.

$N$	$r=1$		$r_{opt}=\frac{4-\alpha}{2-\alpha}$		$r=\frac{4-\alpha}{\alpha}$
	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order
16	2.4091e-02	-	8.1512e-03	-	1.4993e-02	-
32	1.6884e-02	0.513	3.5616e-03	1.195	8.8464e-03	0.761
64	1.1451e-02	0.560	1.5057e-03	1.242	5.0422e-03	0.811
128	7.4624e-03	0.618	6.2368e-04	1.272	2.7729e-03	0.863
Optimal order	1.25

Table 5: Scheme (2.1) for Example 1 (a) with $\alpha=1.75$.

$N$	$r=1$		$r_{opt}=\frac{4-\alpha}{2-\alpha}$		$r=\frac{4-\alpha}{\alpha}$
	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order
16	4.2309e-02	-	2.3254e-02	-	4.1819e-02	-
32	3.3563e-02	0.334	1.1796e-02	0.979	3.2089e-02	0.382
64	2.5699e-02	0.385	5.7768e-03	1.030	2.3768e-02	0.433
128	1.8823e-02	0.449	2.7613e-03	1.065	1.6855e-02	0.496
Optimal order	1.125

Table 6: Scheme (5.1) for Example 1 (a).

$N$	$\alpha=1.01$	$r_{opt}=2$	$\alpha=1.01$	$r=\frac{4-\alpha}{\alpha}$	$\alpha=1.99$	$r_{opt}=2$
	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order
16	1.2195e-04	-	1.8560e-04	-	5.5727e-05	-
32	4.5553e-05	1.421	7.0649e-05	1.394	2.6807e-05	1.056
64	1.6669e-05	1.450	2.6175e-05	1.433	1.2964e-05	1.048
128	5.9895e-06	1.477	9.4821e-06	1.465	6.2090e-06	1.062
Optimal order	1.495				1.005

Table 7: Scheme (2.1) for Example 1 (b) with $\alpha=1.25$.

$N$	$r=1$		$r_{opt}=\frac{4-\alpha}{2-\alpha}$		$r=\frac{4-\alpha}{\alpha}$
	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order
16	1.2467e-02	-	4.2304e-03	-	4.6017e-03	-
32	7.5376e-03	0.726	1.6830e-03	1.330	1.9470e-03	1.241
64	4.4349e-03	0.765	6.5472e-04	1.362	7.9643e-04	1.290
128	2.5261e-03	0.812	2.5069e-04	1.385	3.1707e-04	1.328
Optimal order	1.375

Table 8: Scheme (2.1) for Example 1 (b) with $\alpha=1.5$.

$N$	$r=1$		$r_{opt}=\frac{4-\alpha}{2-\alpha}$		$r=\frac{4-\alpha}{\alpha}$
	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order
16	3.0851e-02	-	1.0472e-02	-	1.9206e-02	-
32	2.1619e-02	0.513	4.5766e-03	1.194	1.1330e-02	0.761
64	1.4660e-02	0.560	1.9351e-03	1.242	6.4568e-03	0.811
128	9.5536e-03	0.618	8.0164e-04	1.271	3.5505e-03	0.863
Optimal order	1.25

Table 9: Scheme (2.1) for Example 1 (b) with $\alpha=1.75$.

$N$	$r=1$		$r_{opt}=\frac{4-\alpha}{2-\alpha}$		$r=\frac{4-\alpha}{\alpha}$
	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order
16	5.4227e-02	-	2.9853e-02	-	5.3599e-02	-
32	4.3015e-02	0.334	1.5145e-02	0.979	4.1127e-02	0.382
64	3.2936e-02	0.385	7.4165e-03	1.030	3.0461e-02	0.433
128	2.4123e-02	0.449	3.5449e-03	1.065	2.1601e-02	0.496
Optimal order	1.125

Table 10: Scheme (2.1) for Example 1 (b) with $\alpha=1.99$.

$N$	$r=1$		$r=5$		$r=10$
	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order
16	7.0532e-03	-	3.0649e-02	-	5.3655e-02	-
32	6.1265e-03	0.203	2.6342e-02	0.218	4.6214e-02	0.215
64	5.1452e-03	0.252	2.1830e-02	0.271	3.8240e-02	0.273
128	4.1324e-03	0.316	1.7280e-02	0.337	3.0184e-02	0.341
Optimal order	1.005

Figure 1: Observation.

Figure 2: $L^{2}$ error.

Figure 3: $\rho_{11}=10^{-3.25}$.

Figure 4: Observation.

Figure 5: $L^{2}$ error.

Figure 6: $\rho_{49}=10^{-3.5}$.

Figure 7: Observation.

Figure 8: $L^{2}$ error.

Figure 9: $\rho_{199}=10^{-3.75}$.

Figure 10: Observation.

Figure 11: $\rho_{199}=3\times 10^{-6}$.

Table 11: Scheme (2.1) for Example 3 with $\alpha=1.25$.

$N$	$r=1$		$r_{opt}=\frac{4-\alpha}{2-\alpha}$		$r=\frac{4-\alpha}{\alpha}$
	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order
20	2.2836e-03	-	1.0643e-03	-	9.0033e-04	-
40	1.3385e-03	0.771	4.2493e-04	1.325	3.6861e-04	1.288
80	7.6092e-04	0.815	1.6494e-04	1.365	1.4657e-04	1.331
160	4.1144e-04	0.887	6.1987e-05	1.412	5.6218e-05	1.383
Optimal order	1.375

Table 12: Scheme (2.1) for Example 3 with $\alpha=1.75$.

$N$	$r=1$		$r_{opt}=\frac{4-\alpha}{2-\alpha}$		$r=\frac{4-\alpha}{\alpha}$
	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order	$e_{L^{2}}(N)$	Order
20	1.4366e-02	-	8.7758e-03	-	1.4216e-02	-
40	1.1009e-02	0.384	4.3576e-03	1.010	1.0541e-02	0.431
80	8.0645e-03	0.449	2.0575e-03	1.083	7.4764e-03	0.496
160	5.5304e-03	0.544	9.4409e-04	1.124	4.9702e-03	0.589
Optimal order	1.125

Figure 12: Exact $a_{1}$.

Figure 13: $\rho_{841}=10^{-5}$.

Figure 14: $\rho_{841}=2\times 10^{-6}$.

Figure 15: $\rho_{841}=10^{-7}$.

6 Concluding remarks

We develop a numerical framework for fractional wave equations under the lower regularity assumptions. The optimal convergence of it is guaranteed by choosing suitable time grids parameter for smooth and nonsmooth solutions. In the numerical framework, a scattered point measurement-based regularization method is used to solve backward problems with uncertain data. The optimal error estimates of stochastic convergence not only balance discretization errors, the noise, and the number of observation points, but also propose an a priori choice of regularization parameters. Despite the presence of large observation errors, we can still obtain more precise inversion results by increasing the number of observation data, even for initial functions with singularity points.

Declarations

On behalf of all authors, the corresponding author states that there is no conflict of interest. No datasets were generated or analyzed during the current study.