A Generalized Iterated Tikhonov Method in the Fourier Domain for Determining the Unknown Source of the Time-Fractional Diffusion Equation

Abstract

In this paper, an inverse problem of determining a source in a time-fractional diffusion equation is investigated. A Fourier extension scheme is used to approximate the solution to avoid the impact on smoothness caused by directly using singular system eigenfunctions for approximation. A modified implicit iteration method is proposed as a regularization technique to stabilize the solution process. The convergence rates are derived when a discrepancy principle serves as the principle for choosing the regularization parameters. Numerical tests are provided to further verify the efficacy of the proposed method.

1. Introduction

Fractional derivatives and integrals capture memory and genetic properties by considering historical dependencies and accumulated past influences, offering a more comprehensive mathematical framework for modeling complex systems in physics, biology, materials science, and genetics. Benefiting from the ability of fractional derivatives and integrals to describe the memory and genetic properties of different substances, replacing the integer derivatives in classical diffusion equations with fractional derivatives can more accurately describe anomalous diffusion phenomena [1,2]. For instance, fractional derivatives enable more accurate modeling of non-Gaussian diffusion which is observed in biological systems [3]; fractional diffusion equations enhance accuracy in modeling environmental transport like pollutant dispersion in groundwater and air [4] and the use of fractional calculus in modeling financial market dynamics improve the accuracy of market models [5]. Therefore, the study of fractional-order diffusion equations has attracted widespread attention in recent years [6,7,8,9,10,11,12,13,14]. Symmetry are intricately linked to the time fractional diffusion equation and can help in solving the equation analytically, determining conservation laws, and developing numerical methods.

Many authors have carried out extensive mathematical studies of direct problems, namely initial value problems and initial boundary value problems for time-fractional differential equations [15,16,17,18,19]. However, in some practical applications, due to condition limitations, we must determine part of the boundary data, initial data, diffusion coefficients, or source terms by additional measurement data, which will yield some inverse problems of fractional-order diffusion. Accurately determining boundary data, initial data, diffusion coefficients, and source terms is crucial across various fields due to its profound impact on modeling accuracy and predictive capabilities, e.g., Accurate diffusion coefficients enable reliable modeling of transport phenomena; precise determination of boundary conditions and initial data is vital for predicting drug distribution within tissues and organs [20]; accurate characterization of boundary conditions, diffusion coefficients, and source terms is essential for modeling fluid flow in porous media, predicting contaminant transport, and assessing natural resource reserves [4]. Therefore, research on such problems has received increased attention in recent years, and its research results have also provided effective solutions to problems in many fields, such as industry, medical care, geophysics, and environmental protection.

In this paper, we consider the following time-fractional diffusion equation ($0<\alpha <1$):

$$\left\{\begin{array}{cc}{}_0{\partial }_t^{\alpha }u-u_{xx}=f\left(x\right)q\left(t\right),\hfill & 0<x<1,0<t<T,\hfill \\ u(0,t)=u(1,t)=0,\hfill & 0\le t\le T,\hfill \\ u(x,0)=0,\hfill & 0\le x\le 1,\hfill \end{array}\right.$$

(1)

${}_0{\partial }_t^{\alpha }u$ here is the usual left-sided Caputo fractional derivative [1], namely

$${}_0{\partial }_t^{\alpha }u=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{\partial u(x,s)}{\partial s}\frac{ds}{{(t-s)}^{\alpha }},0<\alpha <1,$$

where $\Gamma (·)$ is the Gamma function. Problem (1) is forward if the source term $f\left(x\right)h\left(t\right)$ is given. Now, we assume that the time source term $q\in C[0,T]$ is known, and the inverse problem here is to identify the source term $f\left(x\right)$ from a final additional data

$$u(x,T)=g\left(x\right),0\le x\le 1.$$

(2)

We assume the measured data $g^{\delta }\in L^2(0,1)$, which satisfies

$$\parallel g^{\delta }{-g\parallel }_{L^2(0,1)}\le \delta ,$$

(3)

where the constant $\delta >0$ represents noise level.

The increasing interest and research efforts dedicated to solving inverse source problems for fractional diffusion equations stem from the recognition of their significance in various scientific and engineering disciplines. Researchers are increasingly focusing on developing efficient and reliable techniques for solving these inverse source problems [21]. In general, the inverse source problems are ill-posed, i.e., solutions to these problems do not always exist, and even if they do, they do not continuously depend on the given data. Several methods have been developed for solving inverse source problems for fractional diffusion equations, including the truncation method [22], the boundary element method [23], the spectral method [24], the iterative method [25], the Tikhonov regularization method [26,27]. In these methods, the eigenfunctions of the corresponding differential operators are often used to construct approximations of the source terms. This brings about a problem where the eigenfunctions function needs to satisfy specific boundary conditions, which results in the solution being generally less smooth in the sense of regularization theory [28]. To overcome this shortcoming, this paper introduces the Fourier extension method to construct the approximation of the source term. Fourier extension method is devised to ameliorate the Gibbs phenomenon. It extends the function, which is non-periodic, to a function that is periodic at a larger interval. Recently, methods based on Fourier extension have demonstrated highly accurate approximations for non-periodic functions without being restricted by any boundary conditions. The regularization method in the Hilbert scale is introduced to stabilize the solution process. This approach has been successfully used to solve several other ill-posed problems [28,29]. In contrast to the standard Tikhonov method used in previous literature, this paper adopts a generalized iterated Tikhonov method. The implicit iterated Tikhonov regularization method has been well developed for solving the ill-posed problem [30,31,32]. In this paper, we design a specific penalty term in the Fourier domain to overcome the overflow problem of the original method at higher smoothness. We will give error estimates and further verify the efficacy of the method through numerical examples.

This paper is organized as follows: The modified implicit iteration regularization method based on Fourier extension is presented in Section 2. In Section 3, we present the error estimates for the regularization solution. Section 4 provides some numerical tests to confirm the effectiveness of the method. Finally, a simple conclusion is presented in Section 5.

2. The Generalized Iterated Tikhonov Method with Fourier Extension Approximation

Let ${\Omega }_1=(0,1)$, ${\Omega }_2=(0,2)$ and let $L^2\left({\Omega }_i\right)$, $H^s\left({\Omega }_i\right)$$(i=1,2)$ be the usual Lebesgue and Sobolev spaces, respectively. ${\langle ·,·\rangle }_i$ and ${\parallel ·\parallel }_{{\Omega }_i}$ denote the inner products and norm in $L^2\left({\Omega }_i\right)$, ${\parallel ·\parallel }_{s,{\Omega }_i}$ denotes the norm in $H^s\left({\Omega }_i\right)$. Let $H_p^s\left({\Omega }_2\right)$ be the subspace of $H^2\left({\Omega }_2\right)$ containing all functions with period 2. First, we recap some basic definitions and lemmas.

Definition 1

([1]). For $\alpha >0$ and $\beta \in \mathcal{R}$, the generalized Mittag-Leffler function is defined by

$$E_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+\beta )}.$$

(4)

Lemma 1

([33]). For $\lambda >0$ and $0<\alpha <1$, we have

$$\frac{d}{dt}{E}_{\alpha ,1}\left(-\lambda t^{\alpha }\right)=-\lambda t^{\alpha -1}{E}_{\alpha ,\alpha }\left(-\lambda t^{\alpha }\right),t>0.$$

Lemma 2

([22]). If $0<\alpha <1$ and $\eta \ge 0$, we have $E_{\alpha ,1}(-\eta )$ is a monotone decreasing function and

$$1=E_{\alpha ,1}\left(0\right)>E_{\alpha ,1}(-\eta )>0,\eta >0.$$

(5)

The solution of (1) can be obtained by the separation of variables as [25]:

$$u(x,t)=\sum _{\ell =1}^{\infty }\left({\int }_0^{t}q\left(\tau \right){(t-\tau )}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-{\lambda }_{\ell }{(t-\tau )}^{\alpha }\right)d\tau \right){\langle f,X_{\ell }\rangle }_1{X}_{\ell },$$

(6)

where

$${\lambda }_{\ell }={\ell }^2{\pi }^2,X_{\ell }=\sqrt{2}sin\ell \pi x,\ell =1,2,\dots .$$

For ease of notation, denote

$$h_{\ell }\left(T\right)={\int }_0^{T}q\left(\tau \right){(T-\tau )}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-{\lambda }_{\ell }{(T-\tau )}^{\alpha }\right)d\tau ,$$

and substitute the condition (2) into (6), we have

$$g\left(x\right)=\sum _{\ell =1}^{\infty }{h}_{\ell }\left(T\right){\langle f,X_{\ell }\rangle }_1{X}_{\ell },$$

(7)

consequently, we obtain

$$f\left(x\right)=\sum _{\ell =1}^{\infty }\frac{{\langle g,X_{\ell }\rangle }_1{X}_{\ell }}{h_{\ell }\left(T\right)}.$$

Now, if we define an operator $K:L^2\left({\Omega }_1\right)\to L^2\left({\Omega }_1\right)$ as

$$\begin{array}{cc}\hfill Kf\left(x\right)& ={\int }_{{\Omega }_1}k(x,\omega )f\left(\omega \right)d\omega ,\hfill \end{array}$$

(8)

with

$$k(x,\omega )=\sum _{\ell =1}^{\infty }{h}_{\ell }\left(T\right)X_{\ell }\left(x\right)X_{\ell }\left(\omega \right).$$

(9)

Then, from (7), the problem of finding the source term $f\left(x\right)$ from noisy data $g^{\delta }$ can be formulated as the operator equation:

$$Kf\left(x\right)=g^{\delta }\left(x\right).$$

(10)

The ill-posedness of (10) can be obtained by the following Lemma.

Lemma 3.

Let $q\left(t\right)\in C[0,T]$ satisfy

$$0<q_1\le q\left(t\right)\le q_2,\forall t\in [0,T],$$

(11)

then we have

$$\frac{q_1\left(1-E_{\alpha ,1}\left(-{\lambda }_1{T}^{\alpha }\right)\right)}{{\lambda }_{\ell }}\le h_{\ell }\left(T\right)\le \frac{q_2}{{\lambda }_{\ell }}.$$

(12)

Proof. 

Using Lemma 1,

$$\begin{array}{cc}\hfill {\int }_0^T{(T-\tau )}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-{\lambda }_{\ell }{(T-\tau )}^{\alpha }\right)d\tau & =-\frac{1}{{\lambda }_{\ell }}{\int }_0^T-{\lambda }_{\ell }{(T-\tau )}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-{\lambda }_{\ell }{(T-\tau )}^{\alpha }\right)d\tau \hfill \\ \hfill & =-\frac{1}{{\lambda }_{\ell }}{\int }_0^T-\frac{d}{d\tau }{E}_{\alpha ,1}\left(-{\lambda }_{\ell }{(T-\tau )}^{\alpha }\right)d\tau \hfill \\ \hfill & =\frac{1}{{\lambda }_{\ell }}\left(1-E_{\alpha ,1}(-{\lambda }_{\ell }{T}^{\alpha })\right).\hfill \end{array}$$

(13)

Therefore, (12) can be obtained by (5), (11) and (13).    □

Remark 1.

Obviously, if we use $K^{-1}{g}^{\delta }$ to approximate f, the component of error with respect to $X_{\ell }$ will be amplified by $\frac{1}{h_{\ell }\left(T\right)}$. From (12)

$$\frac{1}{h_{\ell }\left(T\right)}\ge \frac{{\lambda }_{\ell }}{q_2}.$$

Note that ${lim}_{\ell \to \infty }\lambda \to \infty$, so the problem (1) is ill-posed.

The Fourier transformation of a function $f\in L^2\left({\Omega }_2\right)$ is

$$\mathcal{F}\left[f\left(x\right)\right]=\sum _{k=-\infty }^{\infty }{\mathbf{f}}_k{e}^{ik\pi x},$$

(14)

with the Fourier coefficients

$${\mathbf{f}}_k=\frac{1}{2}{\int }_{{\Omega }_2}f\left(x\right)e^{-ik\pi x}dx,k=0,\pm 1,\pm 2,\dots .$$

(15)

Accordingly, we can define the operator $\mathcal{F}:{\ell }^2\mathit{\to }{L}^2\left({\Omega }_2\right)$ for any vector $\overrightarrow{\mathbf{v}}={\left({\mathbf{v}}_k\right)}_{k=-\infty }^{\infty }\in {\ell }^2$ as

$$\left(\mathcal{F}\overrightarrow{\mathbf{v}}\right)\left(x\right)=\sum _{k=-\infty }^{\infty }{\mathbf{v}}_k{e}^{ik\pi x}.$$

(16)

Then we can transform Equation (10) into

$$\left(K\mathcal{F}\overrightarrow{\mathbf{v}}\right)\left(x\right)=g^{\delta }\left(x\right),x\in {\Omega }_1.$$

(17)

Let $A=K\mathcal{F}$, then above equation can be rewritten as

$$A\overrightarrow{\mathbf{v}}=g^{\delta }.$$

(18)

The problem to solve (18) can be reformulated as finding the minimum of the functional

$$\overrightarrow{\mathbf{v}}\mapsto {\parallel A\overrightarrow{\mathbf{v}}-g^{\delta }\parallel }^2.$$

(19)

This minimization problem again does not depend continuously on the data, and the Tikhonov regularization in Hilbert adds a penalty term to the functional to restore stability:

$$J_{\alpha }\left(x\right)=\parallel A\overrightarrow{\mathbf{v}}-g^{\delta }{\parallel }^2+\alpha {\parallel B^p\overrightarrow{\mathbf{v}}\parallel }_{{\ell }^2}^2,$$

(20)

where $\alpha >0$ is the regularization parameter, B is an unbounded densely defined self-adjoint strictly positive definite operator, and p is some nonnegative real number. If we suppose that $f\in H^s\left({\Omega }_1\right),s>0$, then we can take $B:{\ell }^2\mathit{\to }{\ell }^2$ as $B\overrightarrow{\mathbf{v}}={\left({(1+{\left|k\right|}^2)}^{\frac{1}{2}}{\mathbf{v}}_k\right)}_{k=-\infty }^{\infty }$. The minimizer of (20) can be obtained by solving the operator equation

$$(A^{*}A+\alpha B^{2p}){\overrightarrow{\mathbf{v}}}_{\alpha }^{\delta }=A^{*}{g}^{\delta }.$$

(21)

If we choose the regularization parameter using the discrepancy principle:

$$\parallel A{\overrightarrow{\mathbf{v}}}_{\alpha }^{\delta }-g^{\delta }{\parallel }_{{\Omega }_1}=C\delta$$

(22)

with $C>1$, then the corresponding convergence rates can be derived using the theories developed in [34]. There are two issues with the above process. First, in the actual process, s is usually unknown, and when s is large, a larger p must be selected to obtain a higher convergence rate, which will cause overflow. In addition, when the error is small, the regularization parameter selected by (22) is very small, which will also cause numerical instability.

To circumvent these two issues, this paper introduces a new operator R and uses a more stable iterated Tikhonov regularization method to solve (18). Specifically, the regularized solution ${\overrightarrow{\mathbf{v}}}_n^{\delta }$ of (18) is defined by

$$\begin{array}{cc}\hfill {\overrightarrow{\mathbf{v}}}_0^{\delta }& =0,\hfill \\ \hfill {\overrightarrow{\mathbf{v}}}_m^{\delta }& ={\overrightarrow{\mathbf{v}}}_{m-1}^{\delta }-{(A^{*}A+{\beta }_m{R}^2)}^{-1}{A}^{*}(A{\overrightarrow{\mathbf{v}}}_{m-1}^{\delta }-g^{\delta }),m=1,2,\cdots ,n,\hfill \end{array}$$

(23)

where the operator $R:{\ell }^2\mathit{\to }{\ell }^2$ is defined by $R\overrightarrow{\mathbf{v}}={\left(e^{\left|k\right|}{\mathbf{v}}_k\right)}_{k=-\infty }^{\infty }$ and ${\beta }_m>0$ are properly chosen real numbers. The number of iterations n will be chosen as the solution of equation:

$$\parallel A{\overrightarrow{\mathbf{v}}}_n^{\delta }-g^{\delta }{\parallel }_{{\Omega }_1}=C\delta$$

(24)

with $C\ge 1$.

Remark 2.

The above formula is inconvenient in the implementation process. Later in the specific numerical process, similar to [31,32], we adopt the following form of discrepancy principle: first take ${\beta }_1=1$, ${\beta }_k=q^{k-1}{\beta }_1$ with some $q<1$ and choose n as the first integer for which

$$\parallel A{\overrightarrow{\mathbf{v}}}_n^{\delta }-g^{\delta }{\parallel }_{{\Omega }_1}\le C\delta <{\parallel A{\overrightarrow{\mathbf{v}}}_k^{\delta }-g^{\delta }\parallel }_{{\Omega }_1},0\le k<n,$$

(25)

with some $C>1$.

Lemma 4

([32]). If we let $T=A{R}^{-1}$,

$${\sigma }_n:=\sum _{m=1}^n\frac{1}{{\beta }_m},d_n\left(\lambda \right):=\frac{1}{\lambda }\left(1-\prod _{m=1}^n\frac{{\beta }_m}{\lambda +{\beta }_m}\right)\mathrm{and}{r}_n\left(\lambda \right):=1-\lambda d_n\left(\lambda \right),$$

(26)

then the regularized solution ${\overrightarrow{\mathbf{v}}}_m^{\delta }$ possesses the representation

$${\overrightarrow{\mathbf{v}}}_m^{\delta }-{\overrightarrow{\mathbf{v}}}_0^{\delta }=R^{-1}{d}_n\left(T^{*}T\right)T^{*}{g}^{\delta }$$

(27)

and we have

$$\begin{array}{ccc}\hfill d_n\left(\lambda \right)\le {\sigma }_n,& \hfill & \hfill \lambda d_n\left(\lambda \right)\le 1,\\ \hfill \lambda r_n\left(\lambda \right)\le {\sigma }_n^{-1},& \hfill & \hfill r_n\left(\lambda \right)\le 1.\end{array}$$

(28)

Lemma 5

([35]). Suppose that $d_n\left(\lambda \right)$ and $r_n\left(\lambda \right)$ is defined by (26), then we have

$$\sqrt{\lambda }{d}_n\left(\lambda \right)\le \sqrt{{\sigma }_n},\sqrt{\lambda }{r}_n\left(\lambda \right)\le \sqrt{{\sigma }_n^{-1}}.$$

(29)

Next, we consider the discretization form of (23). Suppose that the data $g^{\delta }$ are given at nodes ${\left\{x_i\right\}}_{i=1}^M$ and

$$V_K=\mathrm{span}\{e^{ikx}\left|\right|k|\le K\}.$$

Then the discrete form of Equation (18) can be given as

$$\mathsf{A}\mathsf{v}={\mathrm{g}}^{\delta },$$

(30)

where the vector

$$\mathsf{v}={(v_{-K},v_{-K+1},\dots ,v_{K-1},v_K)}^T,{\mathrm{g}}^{\delta }={(g^{\delta }\left(x_1\right),g^{\delta }\left(x_2\right),\dots ,g^{\delta }\left(x_M\right))}^T$$

and the matrix ${\mathsf{A}}_{M\times (2K+1)}$ contains all $A{e}^{ik{x}_i}$ (obtained by the direct solver) as its columns. Moreover, if we let

$$\mathsf{R}=\left(\begin{array}{ccccc}{e}^K& & & & \\ & e^{K-1}& & & \\ & & \ddots & & \\ & & & e^{K-1}& \\ & & & & e^K\end{array}\right),$$

then the discrete form of (23) can be given as:

$$\begin{array}{cc}\hfill {\mathsf{v}}_0^{\delta }& =0,\hfill \\ \hfill {\mathsf{v}}_m^{\delta }& ={\mathsf{v}}_{m-1}^{\delta }-{({\mathsf{A}}^{*}\mathsf{A}+{\beta }_m{\mathsf{R}}^2)}^{-1}{\mathsf{A}}^{*}(\mathsf{A}{\mathsf{v}}_{m-1}^{\delta }-{\mathrm{g}}^{\delta }),m=1,2,\cdots ,n.\hfill \end{array}$$

(31)

We can give a step-by-step version of the algorithm as Algorithm 1:

Algorithm 1: Implicit iterative Tikhonov regularization algorithm

1. Start with initial data ${\mathrm{g}}^{\delta }$, $\delta$, M, K, $C=1.1$, ${\mathsf{v}}_0^{\delta }=0$, ${\beta }_1=1$, $q=\frac{1}{10}$. 2. Compute the matrix $\mathsf{A}$ by the direct solver. 3. if $\parallel \mathsf{A}{\mathsf{v}}_0^{\delta }-{\mathrm{g}}^{\delta }\parallel >C_1\delta$ then 4:    Compute ${\mathsf{v}}^{\delta }:={\mathsf{v}}_0^{\delta }-{({\mathsf{A}}^{*}\mathsf{A}+{\beta }_1{\mathsf{R}}^2)}^{-1}{\mathsf{A}}^{*}(\mathsf{A}{\mathsf{v}}_0^{\delta }-{\mathrm{g}}^{\delta })$ and set $m=1$. 5:   while $\parallel \mathsf{A}{\mathsf{v}}^{\delta }-{\mathrm{g}}^{\delta }\parallel >C\delta$ do 6:       ${\mathsf{v}}_0^{\delta }:={\mathsf{v}}^{\delta }$, ${\beta }_1={\beta }_{1}q$, ${\mathsf{v}}^{\delta }:={\mathsf{v}}_0^{\delta }-{({\mathsf{A}}^{*}\mathsf{A}+{\beta }_1{\mathsf{R}}^2)}^{-1}{\mathsf{A}}^{*}(\mathsf{A}{\mathsf{v}}_0^{\delta }-{\mathrm{g}}^{\delta })$ and set $m=m+1$. 7:   end while

3. Source Conditions and Convergence Rates

We denote the operator

$$\begin{array}{cc}\hfill P_N\overrightarrow{\mathbf{v}}& :=(\dots ,0,0,{\mathbf{v}}_{-N},{\mathbf{v}}_{-N+1},\dots ,{\mathbf{v}}_0,\dots ,{\mathbf{v}}_{N-1},{\mathbf{v}}_N,0,0,\dots ),\hfill \\ \hfill D^s\overrightarrow{\mathbf{v}}& :={\left({(1+k^2)}^{\frac{s}{2}}{\mathbf{v}}_k\right)}_{k=-\infty }^{\infty }.\hfill \end{array}$$

(32)

To obtain the meaningful convergence rates of the solution, some a priori bounded conditions for f are usually needed. Here, we assume that $f\in H^s\left({\Omega }_1\right)$. So, by the extension theorem in Sobolev spaces and Parseval’s formula, we can assume the exact solution $\overrightarrow{\mathbf{f}}$ of (18) satisfies

$$\parallel D^s\overrightarrow{\mathbf{f}}{\parallel }_{{\ell }^2}\le E.$$

(33)

If we let

$${\overrightarrow{\mathbf{f}}}_N=P_N\overrightarrow{\mathbf{f}},g_N=A{\overrightarrow{\mathbf{f}}}_N\mathrm{and}{\overrightarrow{\mathbf{v}}}_{n,N}=R^{-1}{d}_n\left(T^{*}T\right)T^{*}{g}_N,$$

(34)

then it can be deduced that

$$\begin{array}{cccc}\hfill A({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{v}}}_{n,N})& =\hfill & \hfill & T{d}_n\left(T^{*}T\right)T^{*}(g^{\delta }-g_N),\hfill \\ \hfill A(\overrightarrow{\mathbf{f}}-{\overrightarrow{\mathbf{v}}}_{n,N})& =\hfill & \hfill & T{r}_n\left(T^{*}T\right)R{\overrightarrow{\mathbf{f}}}_N,\hfill \\ \hfill g^{\delta }-A{\overrightarrow{\mathbf{v}}}_n^{\delta }& =\hfill & \hfill & r_n\left(T^{*}T\right)g^{\delta },\hfill \\ \hfill R({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{v}}}_{n,N})& =\hfill & \hfill & d_n\left(T^{*}T\right)T^{*}(g^{\delta }-g_N),\hfill \\ \hfill R(\overrightarrow{\mathbf{f}}-{\overrightarrow{\mathbf{v}}}_{n,N})& =\hfill & \hfill & r_n\left(T^{*}T\right)R{\overrightarrow{\mathbf{f}}}_N.\hfill \end{array}$$

(35)

In further studies, the following results are needed.

Lemma 6.

If the condition (33) holds, then we have

$$\parallel \overrightarrow{\mathbf{f}}-P_N\overrightarrow{\mathbf{f}}{\parallel }_{{\ell }^2}\le N^{-s}E,\mathrm{and}{\parallel R{\overrightarrow{\mathbf{f}}}_N\parallel }_{{\ell }^2}\le C_{N}E,$$

(36)

where

$$C_N=max\left(1,\frac{e^N}{N^s}\right).$$

(37)

Proof. 

It is easy to obtain

$$\begin{array}{cc}\hfill \parallel \overrightarrow{\mathbf{f}}-P_N\overrightarrow{\mathbf{f}}{\parallel }_{{\ell }^2}^2& =\sum _{k>N}{\mathbf{f}}_k^2\le N^{-2s}\sum _{k>N}{k}^{2s}{\mathbf{f}}_k^2\le N^{-2s}{E}^2,\hfill \end{array}$$

(38)

and

$$\begin{array}{cc}\hfill \parallel R{\overrightarrow{\mathbf{f}}}_N{\parallel }_{{\ell }^2}^2& =\sum _{k=-N}^N{e}^{2\left|k\right|}{\mathbf{f}}_k^2=\sum _{k=-N}^N\frac{e^{2\left|k\right|}}{{(1+|k|}^2{)}^s}{(1+|k|}^2{)}^s{\mathbf{f}}_k^2\le max\left(1,\frac{e^{2N}}{N^{2s}}\right)E^2.\hfill \end{array}$$

(39)

□

If we denote the operator $\widehat{K}$ as

$$\widehat{K}f=\sum _{\ell =1}^{\infty }\frac{1}{{\lambda }_{\ell }}(f,X_{\ell })X_{\ell },$$

(40)

then from Lemma 3,

$$\frac{1}{q_2}\parallel Kf\parallel \le \parallel \widehat{K}f\parallel \le \frac{1}{q_1(1-E_{\alpha ,1}(-{\pi }^2{T}^{\alpha }))}\parallel Kf\parallel .$$

(41)

And the following result holds.

Lemma 7

([28]). If the condition (33) holds, then there exists a constant M such that

$$\parallel \widehat{K}\mathcal{F}(I-P_N)\overrightarrow{\mathbf{f}}{\parallel }_{{\Omega }_1}\le M{N}^{-\widehat{s}-2}E,$$

(42)

where $\widehat{s}=s-\frac{1}{2}$.

Lemma 8.

If the conditions (3) and (33) hold, we have

$$\parallel A({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{f}}}_N){\parallel }_{{\Omega }_1}\le (C+1)\delta +q_{2}M{N}^{-\widehat{s}-2}E,$$

(43)

$$\parallel R({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{f}}}_N){\parallel }_{{\ell }^2}\le \sqrt{{\sigma }_n}(\delta +q_{2}M{N}^{-\widehat{s}-2}E)+C_{N}E,$$

(44)

and

$$\parallel A{\overrightarrow{\mathbf{v}}}_n^{\delta }-g^{\delta }{\parallel }_{{\Omega }_1}\le \delta +q_{2}M{N}^{-\widehat{s}-2}E+\sqrt{{\sigma }_n^{-1}}{C}_{N}E.$$

(45)

Proof. 

Using the triangle inequality, together with (3), (24), (41) and (42), we can obtain

$$\parallel A({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{f}}}_N){\parallel }_{{\Omega }_1}\le \parallel A{\overrightarrow{\mathbf{v}}}_n^{\delta }-g^{\delta }{\parallel }_{{\Omega }_1}+\parallel g^{\delta }{-g\parallel }_{{\Omega }_1}+{\parallel A(I-P_N)\overrightarrow{\mathbf{f}}\parallel }_{{\Omega }_1}\le (C+1)\delta +q_{2}M{N}^{-\widehat{s}-2}E.$$

(46)

And using the triangle inequality, (3), (29), (35), (41) and (42), we can obtain

$$\begin{array}{cc}\parallel R({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{f}}}_N){\parallel }_{{\ell }^2}\hfill & \le \parallel R({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{v}}}_{n,N}){\parallel }_{{\ell }^2}+{\parallel R({\overrightarrow{\mathbf{v}}}_{n,N}-{\overrightarrow{\mathbf{f}}}_N)\parallel }_{{\ell }^2}\hfill \\ \hfill & =\parallel d_n\left(T^{*}T\right)T^{*}(g^{\delta }-g_N){\parallel }_{{\ell }^2}+{\parallel r_n\left(T^{*}T\right)R{\overrightarrow{\mathbf{f}}}_N\parallel }_{{\ell }^2}\hfill \\ \hfill & \le \sqrt{{\sigma }_n}\parallel (g^{\delta }-g_N){\parallel }_{{\Omega }_1}+{\parallel R{\overrightarrow{\mathbf{f}}}_N\parallel }_{{\ell }^2}\hfill \\ \hfill & \le \sqrt{{\sigma }_n}(\delta +q_{2}M{N}^{-\widehat{s}-2}E)+C_{N}E.\hfill \end{array}$$

(47)

Moreover, in terms of the triangle inequality, (28), (29), (35), (36), (41) and (42), we have

$$\begin{array}{cc}\parallel A{\overrightarrow{\mathbf{v}}}_n^{\delta }-g^{\delta }{\parallel }_{{\Omega }_1}\hfill & =\parallel r_n\left(\mathcal{T}{\mathcal{T}}^{*}\right)g^{\delta }{\parallel }_{{\Omega }_1}\hfill \\ \hfill & \le \parallel r_n\left(\mathcal{T}{\mathcal{T}}^{*}\right)(g^{\delta }-g){\parallel }_{{\Omega }_1}+\parallel r_n\left(\mathcal{T}{\mathcal{T}}^{*}\right)(g-g_N){\parallel }_{{\Omega }_1}+{\parallel r_n\left(\mathcal{T}{\mathcal{T}}^{*}\right)g_N\parallel }_{{\Omega }_1}\hfill \\ \hfill & \le \delta +\parallel g-g_N{\parallel }_{{\Omega }_1}+\parallel r_n\left(\mathcal{T}{\mathcal{T}}^{*}\right)\mathcal{T}\parallel ·\parallel \mathcal{R}{\overrightarrow{\mathbf{f}}}_N\parallel \hfill \\ \hfill & \le \delta +q_{2}M{N}^{-\widehat{s}-2}E+\sqrt{{\sigma }_n^{-1}}{C}_{N}E.\hfill \end{array}$$

(48)

□

Lemma 9

([28]). If the vector sequences ${\overrightarrow{\mathbf{v}}}^{\delta }$ satisfy

$$\parallel \widehat{K}\mathcal{F}{\overrightarrow{\mathbf{v}}}^{\delta }{\parallel }_{{\Omega }_1}\le c_1\delta ,{\parallel R{\overrightarrow{\mathbf{v}}}^{\delta }\parallel }_{{\ell }^2}\le c_2{e}^{c_3{\delta }^{-\frac{1}{\widehat{s}+2}}}{\delta }^{\frac{s}{\widehat{s}+2}},\mathrm{as}\delta \mathit{\to }0,$$

(49)

where $c_1,c_2,c_3$ are some nonnegative real numbers. Then

$$\parallel \mathcal{F}\left(D^{\widehat{s}}{\overrightarrow{\mathbf{v}}}^{\delta }\right){\parallel }_{{\Omega }_1}\le C_0$$

(50)

holds with a constant $C_0$.

Lemma 10

([36]). Let $\Omega =(a,b)\in \mathbb{R}$ and $f\in H^s(\Omega )$, then there exists a constant K depending on ${ϵ}_0$ and j, s, such that for any $0<ϵ\le {ϵ}_0$ and $0\le j\le s$

$${\parallel f\parallel }_{j,\Omega }\le {K(ϵ\parallel f\parallel }_{s,\Omega }+{ϵ}^{-j/(s-j)}{\parallel f\parallel }_{0,\Omega }).$$

(51)

Now, we can derive the following convergence results.

Theorem 1.

If the conditions (3) and (33) hold and ${\overrightarrow{\mathbf{v}}}_n^{\delta }$ is defined by (23) and (24). Let $f_n^{\delta }=\mathcal{F}{\overrightarrow{\mathbf{v}}}_n^{\delta }$, then we have

$$\parallel f_n^{\delta }{-f\parallel }_{{\Omega }_1}=O\left({\delta }^{\frac{\widehat{s}}{\widehat{s}+2}}\right).$$

(52)

Proof. 

If we choose N such that

$$q_{2}M{N}^{-\widehat{s}-2}E=\frac{C-1}{2}\delta ,$$

(53)

that is

$$N={\left(\frac{2q_{2}E}{(C-1)\delta }\right)}^{\frac{1}{\widehat{s}+2}}.$$

(54)

Then from (43) and (45), we can obtain

$$\parallel A({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{f}}}_N){\parallel }_{{\Omega }_1}\le \frac{3C+1}{2}\delta ,$$

(55)

and

$$\frac{C-1}{2}\sqrt{{\sigma }_N}\delta \le C_{N}E.$$

(56)

Therefore, together with (41) and (44), there exist constants $c_1,c_2,c_3$ such that

$$\parallel \widehat{K}\mathcal{F}({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{f}}}_N){\parallel }_{{\Omega }_1}\le c_1\delta ,$$

(57)

and

$$\parallel R({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{f}}}_N){\parallel }_{{\ell }^2}\le c_2{e}^{c_3{\delta }^{-\frac{1}{\widehat{s}+2}}}{\delta }^{\frac{s}{\widehat{s}+2}}.$$

(58)

Thus, from the assertion of Lemma 9, there exist a constant $C_0$

$$\parallel \mathcal{F}\left[D^{\widehat{s}}({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{f}}}_N)\right]{\parallel }_{{\Omega }_1}\le C_0.$$

(59)

Then

$$\begin{array}{cc}{∥\mathcal{F}\left[D^{\widehat{s}}({\overrightarrow{\mathbf{v}}}_n^{\delta }-\overrightarrow{\mathbf{f}})\right]∥}_{{\Omega }_1}\hfill & \le {∥\mathcal{F}\left[D^{\widehat{s}}({\overrightarrow{\mathbf{v}}}_n^{\delta }-{\overrightarrow{\mathbf{f}}}_N)\right]∥}_{{\Omega }_1}+{∥\mathcal{F}\left[D^{\widehat{s}}(\overrightarrow{\mathbf{f}}-{\overrightarrow{\mathbf{f}}}_N)\right]∥}_{{\Omega }_1}\hfill \\ \hfill & \le C_0+{∥\mathcal{F}\left[D^{\widehat{s}}(\overrightarrow{\mathbf{f}}-{\overrightarrow{\mathbf{f}}}_N)\right]∥}_{{\Omega }_2}\hfill \\ \hfill & \le C_0+cE.\hfill \end{array}$$

(60)

Now if we let $h=\widehat{K}\mathcal{F}({\overrightarrow{\mathbf{v}}}_n^{\delta }-\overrightarrow{\mathbf{f}})$, then we have

$${\parallel h\parallel }_{{\Omega }_1}\le \frac{1}{q_1(1-E_{\alpha ,1}(-{\pi }^2{T}^{\alpha }))}\parallel A({\overrightarrow{\mathbf{v}}}_n^{\delta }-\overrightarrow{\mathbf{f}})\parallel \le \frac{1}{q_1(1-E_{\alpha ,1}(-{\pi }^2{T}^{\alpha }))}(C+1)\delta .$$

(61)

Noting that

$${\parallel h\parallel }_{\widehat{s}+2,{\Omega }_1}={∥\mathcal{F}\left[D^{\widehat{s}}({\overrightarrow{\mathbf{v}}}_n^{\delta }-\overrightarrow{\mathbf{f}})\right]∥}_{{\Omega }_1},$$

(62)

so by Lemma 10, this together with

$$\parallel f_n^{\delta }{-f\parallel }_{{\Omega }_1}=\parallel \mathcal{F}({\overrightarrow{\mathbf{v}}}_n^{\delta }-\overrightarrow{\mathbf{f}}){\parallel }_{{\Omega }_1}={\parallel h^{\prime \prime }\parallel }_{{\Omega }_1}$$

(63)

implies the assertion. □

Remark 3.

Regarding the above results and proofs, there are two points to note:

• The parameter N is introduced only for the convenience of the theoretical proof process, and it does not appear in the specific implementation of the method.
• Although the convergence results (52) look similar to those in [25], they are based on different source conditions. The source conditions in [25] are much stricter than condition (33). Compared with condition (33), the source conditions for obtaining the convergence rate in the literature [25] requires that the coefficients of the function f with respect to the eigenfunctions $\left\{X_{\ell }\right\}$ have decay properties, which is only true when the function f satisfies specific boundary conditions. Condition (33) is the smoothness in the sense of general Sobolev space, and there is no boundary constraint on the function f.
• If we use operator $B^p$ instead of R, then from [34] we can obtain a convergence result similar to (52) when $p>(s-2)/2$, but this is numerically difficult to achieve when s is large.

4. Numerical Results and Discussion

We provide some numerical tests to demonstrate the efficacy of the proposed method. The analytic solution of Equation (1) is usually difficult to obtain. Therefore, we need to acquire the data of g through numerical methods. Furthermore, we also need to give the discrete form of Equation (18). We use the spectral method in [37] to achieve these purposes in this paper. Let $x_i=i\Delta x$, $i=0,1,2,\dots ,M,$ then we obtain $g\left(x_i\right)$ by the method in [37] and noise data $g^{\delta }\left(x_i\right)$ are generated as

$$g^{\delta }\left(x_i\right)=g\left(x_i\right)+{\delta }_i,$$

(64)

where ${\left\{{\delta }_i\right\}}_{i=1}^M$ are generated by MATLAB function $\mathrm{randn}(M,1)\ast {\delta }_1$, ${\delta }_1$ is a variable constant that reflects the error level. In our numerical tests, we take $M=129$, $T=1$, $C=1.1$ in (25). To evaluate the accuracy of the method, we compute the relative error of the numerical solution by

$$e_r\left(f\right)={\left(\frac{{\sum }_{i=0}^M{\left(f_n^{\delta }\left(x_i\right)-f\left(x_i\right)\right)}^2}{{\sum }_{i=0}^{M}f{\left(x_i\right)}^2}\right)}^{\frac{1}{2}}.$$

We will also compare the proposed method (M1) with the method (M2) in [25].

Example 1.

We consider the following cases:

(1)

Case 1: Take $q\left(t\right)=e^{-t}$ and

$$f\left(x\right)=e^x.$$

(2)

Case 2: Take $q\left(t\right)=e^{-t}$ and

$$f\left(x\right)=\left\{\begin{array}{cc}3x+1,\hfill & 0\le x<\frac{1}{3},\hfill \\ -3x+3,\hfill & \frac{1}{3}\le x<\frac{2}{3},\hfill \\ 3x-1,\hfill & \frac{2}{3}<x\le 1.\hfill \end{array}\right.$$

(3)

Case 3: Take $q\left(t\right)=t^2+1$ and

$$f\left(x\right)=e^x.$$

(4)

Case 4: Take $q\left(t\right)=t^2+1$ and

$$f\left(x\right)=\left\{\begin{array}{cc}3x+1,\hfill & 0\le x<\frac{1}{3},\hfill \\ -3x+3,\hfill & \frac{1}{3}\le x<\frac{2}{3},\hfill \\ 3x-1,\hfill & \frac{2}{3}<x\le 1.\hfill \end{array}\right.$$

Figure 1 shows how the relative error of the numerical solution changes with the parameter $\alpha$. It can be seen that as $\alpha$ increases, the relative error becomes larger. The change is more significant when $\alpha <0.5$. We also exhibit the changes in relative errors with ${\delta }_1$ in Figure 2. It can be seen that due to the high smoothness of the solutions in Cases 1 and 3, their convergence rates are significantly faster than those in Case 2 and Case 4.

Table 1. $e_r\left(f\right)$ for various ${\delta }_1$ and $\alpha$ (Case 1).

${\mathit{\delta }}_1$	$\mathit{\alpha }=0.1$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.9$
	M1	M2	M1	M2	M1	M2
1 × ${10}^{-1}$	1.02 × ${10}^{-1}$	4.10 × ${10}^{-1}$	1.28 × ${10}^{-1}$	4.82 × ${10}^{-1}$	1.33 × ${10}^{-1}$	4.94 × ${10}^{-1}$
1 × ${10}^{-2}$	1.45 × ${10}^{-2}$	2.92 × ${10}^{-1}$	2.78 × ${10}^{-2}$	3.05 × ${10}^{-1}$	2.94 × ${10}^{-2}$	3.22 × ${10}^{-1}$
1 × ${10}^{-3}$	4.30 × ${10}^{-3}$	2.25× ${10}^{-1}$	7.07 × ${10}^{-3}$	2.53 × ${10}^{-1}$	8.13 × ${10}^{-3}$	2.71 × ${10}^{-1}$
1 × ${10}^{-4}$	6.54 × ${10}^{-4}$	1.98× ${10}^{-1}$	1.59× ${10}^{-3}$	2.13 × ${10}^{-1}$	2.11× ${10}^{-3}$	2.20 × ${10}^{-1}$

,

Table 2. $e_r\left(f\right)$ for various ${\delta }_1$ and $\alpha$ (Case 2).

${\mathit{\delta }}_1$	$\mathit{\alpha }=0.1$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.9$
	M1	M2	M1	M2	M1	M2
1 × ${10}^{-1}$	1.20 × ${10}^{-1}$	3.25 × ${10}^{-1}$	1.84 × ${10}^{-1}$	3.67 × ${10}^{-1}$	1.73 × ${10}^{-1}$	3.83 × ${10}^{-1}$
1 × ${10}^{-2}$	2.06 × ${10}^{-2}$	2.41 × ${10}^{-1}$	3.18 × ${10}^{-2}$	2.81 × ${10}^{-1}$	3.30 × ${10}^{-2}$	3.08 × ${10}^{-1}$
1 × ${10}^{-3}$	1.35 × ${10}^{-2}$	1.87 × ${10}^{-1}$	1.35 × ${10}^{-2}$	2.19 × ${10}^{-1}$	2.13 × ${10}^{-2}$	2.77 × ${10}^{-1}$
1 × ${10}^{-4}$	9.04 × ${10}^{-3}$	1.45 × ${10}^{-1}$	1.09 × ${10}^{-2}$	1.74 × ${10}^{-1}$	1.41 × ${10}^{-2}$	2.01 × ${10}^{-1}$

,

Table 3. $e_r\left(f\right)$ for various ${\delta }_1$ and $\alpha$ (Case 3).

${\mathit{\delta }}_1$	$\mathit{\alpha }=0.1$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.9$
	M1	M2	M1	M2	M1	M2
1 × ${10}^{-1}$	1.13 × ${10}^{-1}$	3.81 × ${10}^{-1}$	1.13 × ${10}^{-1}$	4.01 × ${10}^{-1}$	1.13 × ${10}^{-1}$	4.21 × ${10}^{-1}$
1 × ${10}^{-2}$	1.53 × ${10}^{-2}$	2.79 × ${10}^{-1}$	1.83 × ${10}^{-2}$	2.82 × ${10}^{-1}$	1.84 × ${10}^{-2}$	3.05 × ${10}^{-1}$
1 × ${10}^{-3}$	2.77 × ${10}^{-3}$	2.05 × ${10}^{-1}$	3.82 × ${10}^{-3}$	2.14 × ${10}^{-1}$	6.92 × ${10}^{-3}$	2.37 × ${10}^{-1}$
1 × ${10}^{-4}$	1.95 × ${10}^{-3}$	1.79 × ${10}^{-1}$	2.45 × ${10}^{-3}$	1.86 × ${10}^{-1}$	2.51 × ${10}^{-3}$	1.97 × ${10}^{-1}$

and

Table 4. $e_r\left(f\right)$ for various ${\delta }_1$ and $\alpha$ (Case 4).

${\mathit{\delta }}_1$	$\mathit{\alpha }=0.1$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.9$
	M1	M2	M1	M2	M1	M2
1 × ${10}^{-1}$	1.92 × ${10}^{-1}$	3.05 × ${10}^{-1}$	2.12 × ${10}^{-1}$	3.12 × ${10}^{-1}$	2.32 × ${10}^{-1}$	3.25 × ${10}^{-1}$
1 × ${10}^{-2}$	2.21 × ${10}^{-2}$	2.67 × ${10}^{-1}$	2.30 × ${10}^{-2}$	2.75 × ${10}^{-1}$	2.32 × ${10}^{-2}$	2.88 × ${10}^{-1}$
1 × ${10}^{-3}$	1.84 × ${10}^{-2}$	2.02 × ${10}^{-1}$	1.81 × ${10}^{-2}$	2.23 × ${10}^{-1}$	2.17 × ${10}^{-2}$	2.53 × ${10}^{-1}$
1 × ${10}^{-4}$	4.35 × ${10}^{-3}$	1.81 × ${10}^{-1}$	7.38 × ${10}^{-3}$	1.99 × ${10}^{-1}$	8.23 × ${10}^{-3}$	2.34 × ${10}^{-1}$

lists the relative errors of M1 and M2 for various ${\delta }_1$ and $\alpha$. It can be seen that different q has no substantial impact on the error, and the parameter $\alpha$ has no significant effect on the convergence rate. Moreover, since the solutions do not meet the boundary conditions required by the eigenfunctions, the results of M1 are significantly better than M2.

5. Conclusions

In this paper, we present a reconstruction method for the source term of a time-fractional diffusion equation. The Fourier extension approximation is introduced to overcome the shortcomings of directly using singular system functions. We present a generalized iterated Tikhonov scheme to stabilize the solution process, and the error estimates are obtained with more natural source conditions. Numerical results for several different situations further verify the effectiveness of the method. Moreover, we point out that the framework of the proposed method can be applied to multi-dimensional inverse source problems and other ill-posed problems, which will be discussed in a future paper.