Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Galerkin spectral method for a multi-term time-fractional diffusion equation and an application to inverse source problem

  • In this paper, we employ the Galerkin spectral method to handle a multi-term time-fractional diffusion equation, and investigate the numerical stability and convergence of the proposed method. In addition, we find an interesting application of the Galerkin spectral method to solving an inverse source problem efficiently from the noisy final data in a general bounded domain, and the uniqueness and the ill-posedness for the inverse problem are proved based on expression of the solution. Furthermore, we compare the calculation results of spectral method and finite difference method without any regularization method, and get a norm estimate of the coefficient matrix of a spectral method discretized. And for that we conclude that the spectral method itself can act as a regularization method for some inverse problem (such as inverse source problem). Finally, several numerical examples are used to illustrate the effectiveness and accuracy of the method.

    Citation: L.L. Sun, M.L. Chang. Galerkin spectral method for a multi-term time-fractional diffusion equation and an application to inverse source problem[J]. Networks and Heterogeneous Media, 2023, 18(1): 212-243. doi: 10.3934/nhm.2023008

    Related Papers:

    [1] Kexin Li, Hu Chen, Shusen Xie . Error estimate of L1-ADI scheme for two-dimensional multi-term time fractional diffusion equation. Networks and Heterogeneous Media, 2023, 18(4): 1454-1470. doi: 10.3934/nhm.2023064
    [2] Farman Ali Shah, Kamran, Dania Santina, Nabil Mlaiki, Salma Aljawi . Application of a hybrid pseudospectral method to a new two-dimensional multi-term mixed sub-diffusion and wave-diffusion equation of fractional order. Networks and Heterogeneous Media, 2024, 19(1): 44-85. doi: 10.3934/nhm.2024003
    [3] Feiyang Peng, Yanbin Tang . Inverse problem of determining diffusion matrix between different structures for time fractional diffusion equation. Networks and Heterogeneous Media, 2024, 19(1): 291-304. doi: 10.3934/nhm.2024013
    [4] Diandian Huang, Xin Huang, Tingting Qin, Yongtao Zhou . A transformed $ L1 $ Legendre-Galerkin spectral method for time fractional Fokker-Planck equations. Networks and Heterogeneous Media, 2023, 18(2): 799-812. doi: 10.3934/nhm.2023034
    [5] Patrick Henning, Mario Ohlberger . The heterogeneous multiscale finite element method for advection-diffusion problems with rapidly oscillating coefficients and large expected drift. Networks and Heterogeneous Media, 2010, 5(4): 711-744. doi: 10.3934/nhm.2010.5.711
    [6] Leqiang Zou, Yanzi Zhang . Efficient numerical schemes for variable-order mobile-immobile advection-dispersion equation. Networks and Heterogeneous Media, 2025, 20(2): 387-405. doi: 10.3934/nhm.2025018
    [7] Yinlin Ye, Hongtao Fan, Yajing Li, Ao Huang, Weiheng He . An artificial neural network approach for a class of time-fractional diffusion and diffusion-wave equations. Networks and Heterogeneous Media, 2023, 18(3): 1083-1104. doi: 10.3934/nhm.2023047
    [8] Caihong Gu, Yanbin Tang . Global solution to the Cauchy problem of fractional drift diffusion system with power-law nonlinearity. Networks and Heterogeneous Media, 2023, 18(1): 109-139. doi: 10.3934/nhm.2023005
    [9] Fermín S. V. Bazán, Luciano Bedin, Mansur I. Ismailov, Leonardo S. Borges . Inverse time-dependent source problem for the heat equation with a nonlocal Wentzell-Neumann boundary condition. Networks and Heterogeneous Media, 2023, 18(4): 1747-1771. doi: 10.3934/nhm.2023076
    [10] Panpan Xu, Yongbin Ge, Lin Zhang . High-order finite difference approximation of the Keller-Segel model with additional self- and cross-diffusion terms and a logistic source. Networks and Heterogeneous Media, 2023, 18(4): 1471-1492. doi: 10.3934/nhm.2023065
  • In this paper, we employ the Galerkin spectral method to handle a multi-term time-fractional diffusion equation, and investigate the numerical stability and convergence of the proposed method. In addition, we find an interesting application of the Galerkin spectral method to solving an inverse source problem efficiently from the noisy final data in a general bounded domain, and the uniqueness and the ill-posedness for the inverse problem are proved based on expression of the solution. Furthermore, we compare the calculation results of spectral method and finite difference method without any regularization method, and get a norm estimate of the coefficient matrix of a spectral method discretized. And for that we conclude that the spectral method itself can act as a regularization method for some inverse problem (such as inverse source problem). Finally, several numerical examples are used to illustrate the effectiveness and accuracy of the method.



    Let Ω be a bounded domain in Rd with sufficient smooth boundary Ω. We consider the following initial boundary value problem (IBVP) for the anomalous diffusion model with multi-term time-fractional derivatives

    {sj=1qjαj0+u(x,t)(Lu)(x,t)=f(x)p(t),(x,t)Ω×I,u(x,0)=0,xΩ,u(x,t)=0,(x,t)Ω×I, (1.1)

    where I=(0,T) with T>0 be fixed and αj0+ is the Caputo fractional derivative defined by

    αj0+u(x,t)=1Γ(1αj)t0u(x,s)sds(ts)αj,0<αj<1.

    For a fixed positive integer s, the orders α=(α1,...,αs) and the coefficients q=(q1,...,qs) are restricted in the admissible sets

    B:={(α1,...,αs)Rs;¯αα1>α2>>αsα_}, (1.2)
    Q:={(q1,...,qs)Rs;q1=1,qj[q_,¯q],(j=2,...,s)}, (1.3)

    with fixed 0<α_<¯α<1 and 0<q_<¯q.

    L is a symmetric uniformly elliptic operator defined on D(L)=H2(Ω)H10(Ω) given by

    Lu(x,t)=di=1xi(dj=1aij(x)xju(x,t))+c(x)u(x,t),xΩ,

    in which the coefficients satisfy

    aij=aji,1i,jd,aijC1(ˉΩ),μdi=1ξ2idi,j=1aij(x)ξiξj,xˉΩ,ξ=(ξ1,,ξd)Rd,foraconstantμ>0,c(x)0,xˉΩ,c(x)C(ˉΩ).

    If the number of fractional derivatives s, the orders α and its coefficients q, elliptic operator L and source functions p(t),f(x) are given appropriately, then IBVP Eq (1.1) is a direct problem. Here the spatial source term f(x) is unknown, so the inverse problem in the paper is to determine the spatial source term f(x) based on problem Eq (1.1) and an additional terminal data

    u(x,T)=g(x),xΩ. (1.4)

    Since the measurement is noise-contaminated inevitably, we denote the noisy measurement of g as gδ(x) satisfying

    gδ(x)g(x)δ. (1.5)

    It is well known that the time-fractional diffusion equations (TFDEs) have a wide range of applications in physics, chemistry and other aspects [1,2,3,4,5]. The most representative is continuous time random walk problem in general non-Markovian processes. However, with the increasing demand on the accuracy of the problem, the single-time fractional diffusion equation gradually failed to meet the needs of the problem, so Schumer et al. [6] proposed using the multi-term time-fractional diffusion equation (MTFDE) to increase the accuracy of the model. MTFDE is not only a useful tool for describing the behavior of anomalous diffusion phenomena in highly heterogeneous aquifers and complex viscoelastic materials [7], but can also be applied indirectly in the numerical solution of distributed-order fractional differential equations [8].

    There are many studies on the direct problem for IBVP Eq (1.1) in recent years such as some uniqueness and existence results [9], the maximum principle [10], and analytic solutions [11,12]. In terms of numerical work, there are also some papers on the numerical solutions by finite difference methods [13,14,15] and by finite element methods [16,17,18,19].

    Due to the spectral method needs fewer grid points to obtain highly accurate solution. Therefore, the spectral method is more suitable for discrete MTFDE because of its huge computations. The study of the spectral method can be divided into two categories in MTFDE. On the one hand, the time fractional derivative is discreted by finite difference method and spectral method is applied for the space variable. Guo et al. [20] studied the numerical approximation of distributed-order time-space fractional reaction-diffusion equation. Jiang et al. [21] considered a Legendre spectral method on graded meshes for handling MTFDE with non-smooth solutions in two-dimensional case. Liu et al. [22] applied an alternating direction implicit Legendre spectral method to handle the multi-term time fractional Oldroyd-B fluid type diffusion equation in the two-dimensional case. On the other hand, the spectral method is applied both in space and in time. Zheng et al. [23] considered a high order numerical method for handling MTFDE. Zaky [24] applied a Legendre spectral tau method to handle MTFDE. However, we are interested not only in obtaining a numerical solution with high precision by using spectral method, but also in determining the spatial source term in an inverse problem of the MTFDE.

    As far as we know, the theory as well as numerical methods of inverse source problems for single term (i.e, . s=1) case in Eq (1.1) are relatively abundant. Zhang et al. [25] testified a uniqueness result to inverse the spatial source term in one-dimensional case by using one-point Cauchy data and proposed an efficient numerical method. Wei et al. [26] studied to identify a spatial source term in a multi-dimensional time-fractional diffusion equation from boundary measured data. The uniqueness for the inverse source problem is proved by the Laplace transformation method. Yan et al. [27] studied to identify a spatial source term in a multi-dimensional time-fractional diffusion-wave equation from a part of noisy boundary data. The uniqueness of inverse spatial source term problem is proved by the Titchmarsh convolution theorem and the Duhamel principle. Sun et al. [28] devoted to recovering simultaneously the fractional order and the space-dependent source term from partial Cauchy's boundary data in a multi-dimensional time-fractional diffusion equation. Recently, Yeganeh et al. [29] came up with an interesting idea. They used a method based on a finite difference scheme in time and a local discontinuous Galerkin method in space to determine a spatial source term in a time-fractional diffusion equation. This has not been seen in previous studies.

    For a multi-term case, however, research results on inverse source problem are relatively few at present. Jiang et al. [30] built a weak unique continuation property for time-fractional diffusion-advection equations, and they considered an inverse problem on determining identifying the spatial source term by interior measurements. Li et al. [31] considered an inverse problem of recovering time-dependent source term from the Cauchy data in a MTFDE, and they applied the conjugate gradient method to identify the approximate source term. Recently, Sun et al. [32] investigated an inverse the spatial source term in MTFDE with nohomogeneous boundary condition from partially disturbed boundary data. They proposed the Levenberg-Marquardt regularization method to compute an inverse source problem. In addition, simultaneous inversion of source term and other terms has been studied. For instance, Malik et al. [33] studied an inverse problem of identifying a time-dependent source term along with diffusion/temperature concentration from a non-local over-specified condition for a space-time fractional diffusion equation. Sun et al. [34] considered a nonlinear inverse problem for simultaneously recovering the potential function and the fractional orders in a MTFDE from the noisy boundary Cauchy data in the one-dimensional case. For other inverse source problems, we can refer to [35,36,37]. Nevertheless, to the best knowledge of the authors, no one has used the spectral method to determine a spatial source term in an inverse problem of MTFDE.

    In this paper, we will focus on two aspect research on our proposed model Eq (1.1). One is the study on the direct problem. The theoretical analysis and the numerical scheme of the Galerkin spectral method are proposed to solve the IBVP Eq (1.1). On the other hand, we use the Galerkin spectral method to investigate an inverse spatial source term problem in a MTFDE from a noisy final data in a general bounded domain. In the propose method, multi-term Caputo fractional derivatives are discretized by the L1-formula and the Galerkin spectral method is applied for the space variable. At the same time, the comparison between the Galerkin spectral method and the finite difference method is added in this paper. Finally, we not only prove the validity of the Galerkin spectral method in the application of MTFDE, but also verify the superiority of the Galerkin spectral method in forward and inverse problems by comparing the numerical results with the finite difference method.

    The remainder of this paper is organized as follows. Some preliminaries are presented in Section 2. The detailed convergence analysis of the presented method is shown in Section 3. Uniqueness and ill-posedness for the inverse problem are showed in Section 4. We present the Galerkin spectral method and the finite difference method algorithm in Section 5. Numerical results for four examples are investigated in Section 6. Finally, we give a conclusion in Section 7.

    We firstly introduce some preliminaries as follows in this section.

    Definition 2.1. The multinomial Mittag-Leffler function is defined by (see [9,38])

    E(θ1,,θs),θ0(z1,,zs):=k=0k1++ks=k(k;k1,,ks)sj=1zkjjΓ(θ0+sj=1θjkj),

    where θ0,θjR, and zjC (j=1,,s), and (k;k1,,ks) denotes the multinomial coefficient

    (k;k1,,ks):=k!k1!ks!withk=sj=1kj,

    where kj(j=1,,s) are non-negative integers.

    For the convenience of the later, if the orders α=(α1,...,αs) and its coefficients q=(q1,...,qs) satisfy Eq (1.2) and Eq (1.3), then we adopt the abbreviation

    E(n)α,β(t):=E(α1,α1α2,,α1αs),β(λntα1,q2tα1α2,,qstα1αs),t>0,n=1,2,,

    where α=(α1,α1α2,,α1αs) and λn denotes the n-th eigenvalues of elliptic operator L with the homogeneous Dirichlet boundary condition.

    Lemma 2.2. ([9]) Let 0<αs<αs1<<α1<1. Then,

    ddt{tα1E(n)α,1+α1(t)}=tα11E(n)α,α1(t),t>0.

    Lemma 2.3. ([9,39]) Let 0<αs<αs1<<α1<1. Then the function tα11E(n)α,α1(t) is positive for t>0.

    Lemma 2.4. ([40]) Let 0<αs<αs1<<α1<1. Then

    |1λntα1E(n)α,α1+1(t)|sj=2M(1+qjtα1αj)1+λntα1,t>0,n=1,2,,

    where M ia a positive constant.

    Proposition 2.5. Let 0<αs<αs1<<α1<1. Then we have E(n)α,1+α1(t)>0.

    Proof. By Lemma 2.2 and 2.3, we know that

    ddt{tα1E(n)α,1+α1(t)}=tα11E(n)α,α1(t)>0.

    Hence it is obvious E(n)α,1+α1(t)>0.

    Proposition 2.6. For λn>0 and 0<αs<αs1<<α1<1, we have 0<1λntα1E(n)α,1+α1(t)<1 for t>0. Moreover, 1λntα1E(n)α,1+α1(t) is a strictly decreasing function on t>0.

    Proof. By Lemma 2.2 and 2.3, we have

    ddt{1λntα1E(n)α,1+α1(t)}=λntα11E(n)α,α1(t)<0.

    We notice that 1λntα1E(n)α,1+α1(t) is a continuous function on t. Hence, we have limt0(1λntα1E(n)α,1+α1(t))=1. By Lemma 2.4 we know limt(1λntα1E(n)α,1+α1(t))=0. The proof is completed.

    Lemma 2.7. For λn>0 and 0<αs<αs1<<α1<1, such that

    0E(n)α,1+α1(T)1Tα1λn.

    Proof. By Proposition 2.6, we know 0<E(n)α,1+α1(T)<1λnTα1.

    Lemma 2.8. Let k be a positive integer, Δt>0 and limkkΔt=T. We have the following properties from the sequence defined by

    ωk=sj=1qjΔtαjΓ(2αj)(k1αj(k1)1αj), (2.1)

    where αj and qj are defined in equation Eq (1.1). Here we denote ω1k equivalent to 1ωk.

    (1) ωk is a decreasing sequence with respect to k;

    (2) ω1>>ωk>0 for each k;

    (3) ω1k is an increasing sequence with respect to k;

    (4) limkkα1ω1k=Δtα1q1(1α1)Γ(2α1)+sj=2qjTα1αj(1αj)Γ(2αj).

    Proof. First, let Φ(x)=ωx be a continuous function with respect to x, According to dΦ(x)dx<0, we know Φ(x) is decreasing function, then we have ωk is decreasing the sequence. And we can direct to check that ω1>>ωk>0 for each k. In addition, we can easy to derive ω1k is increasing function. Further, we have

    limkkα1ω1k=limkkα1q1Δtα1Γ(2α1)(k1α1(k1)1α1)+sj=2qjΔtαjΓ(2αj)(k1αj(k1)1αj)=:limk1Ψ1k+Ψ2k,

    where

    Ψ1k=q1Δtα1Γ(2α1)(k1α1(k1)1α1)kα1,Ψ2k=sj=2qjΔtαjΓ(2αj)(k1αj(k1)1αj)kα1.

    First we have

    limkΨ1k=limkq1Δtα1Γ(2α1)(1(11k)1α1)k1α1kα1=limkq1Δtα1Γ(2α1)limk(11+(1α1)1k+o(1k))k=q1(1α1)Γ(2α1)Δtα1,

    where we use the Taylor's formula of (11/x)1α1 in the second equation. Similarly, we have

    limkΨ2k=limksj=2qjΔtαjΓ(2αj)(k1αj(k1)1αj)kα1=limksj=2qjΔtαjΓ(2αj)(k1αj(k1)1αj)kαjkαjΔtα1Δtα1kα1=limksj=2qj1Γ(2αj)(1(11k)1αj)k(kΔt)αj(kΔt)α1Δtα1=sj=2qjTα1αj(1αj)Γ(2αj)Δtα1.

    For all these reasons, we arrive at

    limkkα1ω1k=Δtα1q1(1α1)Γ(2α1)+sj=2qjTα1αj(1αj)Γ(2αj). (2.2)

    This proof is complete.

    In this section, we discuss semi-discretize and full discretize schemes of the problem Eq (1.1) in the one-dimensional case. First, let tn:=nΔt,n=0,1,,N, where Δt:=TN (at leat satisfying Δt<1) is the time step. We can denote un an approximation to u(x,tn), and in order to convenience we write p(tn),f(x),a(x),c(x) is equal to pn,f,a,c, where f(x),p(tn) are the source terms of Eq (1.1), and a(x),c(x) come from the definition of the symmetric uniformly elliptic operator L. Then we introduce the time-fractional derivative is approximated by the L1 formula [41].

    sj=1qjαj0+u(x,tn)=nk=1ωk(u(x,tnk+1)u(x,tnk))+rnΔt,

    where ωk=sj=1qjΔtαjΓ(2αj)(k1αj(k1)1αj) and rnΔt is the truncation error with the estimate

    rnΔtcuΔt2α1++cuΔt2αsc1Δt2α1.

    Here c1 is a constant depending on u,aj and T.

    In this paper, we use H1-norm defined by

    νH1(Ω)=(ν2L2(Ω)+1ω1dνdx2L2(Ω))12.

    where ω1 is defined in Lemma 2.8.

    In order to establish the complete semi-discrete problem, we define the error term rn by

    rn:=1ω1[nj=1ωj(u(x,tnj+1)u(x,tnj))sj=1qjαj0+u(x,tn)]. (3.1)

    Then we have

    |rn|=1ω1|rnΔt|c2Δt2α1, (3.2)

    where c2 is a constant depending on u,aj and T.

    Using the first term on right hand side of Eq (3.1) as an approximation of sj=1qjαj0+u(x,tn) leads to the following finite difference scheme to Eq (1.1)

    nk=1ωk(unk+1unk)(aunx)xcun=fpn. (3.3)

    Then we have

    un1ω1(aunx)x1ω1cun=n1k=1ωkωk+1ω1unk+ωnω1u0+1ω1fpn. (3.4)

    Multiply both sides by νH10(Ω) and integrate over Ω and by Green formula. We can find unH10(Ω), such that

    (un,ν)+1ω1(aunx,νx)1ω1(cun,ν)=n1k=1ωkωk+1ω1(unk,ν)+ωnω1(u0,ν)+1ω1(fpn,ν). (3.5)

    For convenience and without loss of generality, we consider the case aij(x)1,c(x)=0. Stability of the scheme Eq (3.5) is given in the following result.

    Theorem 3.1. The semi-discretized problem Eq (3.5) is unconditionally stable in sense that for all Δt>0, holds

    unH1(Ω)u0L2(Ω)+1ω1nk=1pkfL2(Ω).

    Proof. we prove this theorem using induction. When n = 1, scheme Eq (3.5) has

    (u1,ν)+1ω1(u1x,νx)=(u0,ν)+1ω1(1k=1p1f,ν).

    Taking ν=u1 and Schwarz inequality, we have

    u1H1(Ω)u0L2(Ω)+1ω11k=1pkfL2(Ω). (3.6)

    Now for the following inequality to hold

    ujH1(Ω)u0L2(Ω)+1ω1jk=1pkfL2(Ω),j=1,2,,n1. (3.7)

    We must prove that unH1(Ω)u0L2(Ω)+1ω1nk=1pkfL2(Ω). Taking ν=un in Eq (3.5) and by using Eq (3.7), we have

    unH1(Ω)n1k=1ωkωk+1ω1unkL2(Ω)+ωnω1u0L2(Ω)+1ω1fL2(Ω)pn[n1k=1ωkωk+1ω1+ωnω1][u0L2(Ω)+1ω1n1k=1pkfL2(Ω)]+1ω1pnfL2(Ω). (3.8)

    Hence, we have unH1(Ω)u0L2(Ω)+1ω1nk=1pkfL2(Ω).

    Now, we give the convergence estimate between the exact solution and the solution of the semi-discretized problem.

    Theorem 3.2. Let u(x,t) be the exact solution of Eq (1.1), {un}Nn=0 be the numerical solution of semi-discretized Eq (3.5), and u(x,t0)=u0. Then there holds the following error estimates

    u(x,tn)unH1(Ω)c2ω1Tα1Δt2α1q1(1α1)Γ(2α1)+sj=2qjTα1αj(1αj)Γ(2αj),

    where c2 is a constant depending on u,aj and T.

    Proof. We will verify that the following estimate

    u(x,tj)ujH1(Ω)c2ω1ω1jΔt2α1,j=1,2,,N. (3.9)

    Here, we will prove it by the mathematical induction, let ˆen=u(x,tn)un. For j=1, we have, by combining Eq (1.1), Eq (3.5), Eq (3.1), the error equation

    (ˆe1,ν)+1ω1(ˆe1x,νx)=(ˆe0,ν)+(r1,ν)νH10(Ω). (3.10)

    Taking ν=ˆe1, we have

    u(x,t1)u1H1(Ω)c2ω1ω11Δt2α1.

    Hence, Eq (3.9) is proven for the case j = 1. Now, Suppose the inequality Eq (3.9) holds for all j=1,2,,n1. Then we need to prove that it holds also for j=n.

    Similar to j=1, we have νH10(Ω)

    (ˆen,ν)+1ω1(ˆenx,νx)=n1k=1ωkωk+1ω1(ˆenk,ν)+ωnω1(ˆe0,ν)+(rn,ν). (3.11)

    Taking ν=ˆen, and using Eq (3.9). We have

    ˆenH1(Ω)n1k=1ωkωk+1ω1ˆenkL2(Ω)+ωnω1ˆe0L2(Ω)+rnL2(Ω)[n1k=1ωkωk+1ω1+ωnω1]c2ω1ω1nΔt2α1.

    Then we obtain

    ˆenH1(Ω)c2ω1ω1nΔt2α1.

    The estimate Eq (3.9) is proved. Notice that limkkΔt=T is clearly true. Then, we have

    u(x,tn)unH1(Ω)c2ω1ω1nΔt2α1=c2ω1nα1ω1nΔtα1Tα1Δt2α1. (3.12)

    Then by Eq (3.12) and Lemma 2.8, we deduce that for sufficently large n, the assertion of the theorem is valid.

    First, we set

    VM={νH10(Ω)|νPM(Ω)},

    as the space of polynomials of degree M with respect to x. Let P be the orthogonal projection operator from H10(Ω) into VM. Then for all ψH10(Ω), we have PψVM. Further, we have νMVM such that

    (Pψ,νM)+1ω1(Pψx,νMx)=(ψ,νM)+1ω1(ψx,νMx). (3.13)

    This the orthogonal projection operator satisfy the inequality [42]

    ψPψH1(Ω)cM1mψHm(Ω), (3.14)

    where ψHm(Ω)H10(Ω),m1.

    Now, we consider the Galerkin scheme as follow: find unMVM, such that for all νMVM

    (unM,νM)+1ω1(unMx,νMx)=n1k=1ωkωk+1ω1(unkM,νM)+ωnω1(u0M,νM)+pnω1(fM,νM). (3.15)

    Theorem 3.3. Let u(x,t)H1((0,T),Hm(Ω)H10(Ω)),m>1 be the exact solution of Eq (1.1), {unM}Nn=0 be the numerical solution of full discretization Eq (3.15), and u0M=Pu0. Then there holds the following error estimates

    u(x,tn)unMH1(Ω)ω1Tα1(c2Δt2α1+cM1muL(Hm))q1(1α1)Γ(2α1)+sj=2qjTα1αj(1αj)Γ(2αj),

    where c2 is a constant depending on u,aj,T, and uL(Hm):=supt(0,T)u(x,t)Hm(Ω).

    Proof. From Eq (1.1), Eq (3.1) satisfy νH10(Ω),

    (u(x,tn),ν)+1ω1(u(x,tn)x,νx)=n1k=1ωkωk+1ω1(u(x,tnk),ν)+ωnω1(u(x,t0),ν)+(rn,ν)+pnω1(f,ν). (3.16)

    By projecting u(x,tn) in Pu(x,tn)VM, applying Eq (3.13), we obtain for all νMVM

    (Pu(x,tn),νM)+1ω1(Pu(x,tn)x,νMx)=n1k=1ωkωk+1ω1(u(x,tnk),νM)+ωnω1(u(x,t0),νM)+(rn,νM)+pnω1(f,νM). (3.17)

    Let ˆϵnM=Pu(x,tn)unM,ϵnM=u(x,tn)unM by subtracting Eq (3.15) from Eq (3.17), we have

    (ˆϵnM,νM)+1ω1(ˆϵnMx,νMx)=n1k=1ωkωk+1ω1(ϵnkM,νM)+ωnω1(ϵ0M,νM)+(rn,νM) (3.18)

    Taking νM=ˆϵnM in above equation and using the triangular inequality ϵnMH1(Ω)ˆϵnMH1(Ω)+u(x,tn)Pu(x,tn)H1(Ω), we obtain

    ϵnMH1(Ω)n1k=1ωkωk+1ω1ϵnkMH1(Ω)+rnH1(Ω)+u(x,tn)Pu(x,tn)H1(Ω)n1k=1ωkωk+1ω1ϵnkMH1(Ω)+c2Δt2α1+cM1mu(x,tn)Hm(Ω). (3.19)

    Here, we use the mathematical induction, by Eq (3.19), we can easy to verify when case of n = 1.

    ϵ1MH1(Ω)ω1ω11(c2Δt2α1+cM1mu(x,t1)Hm(Ω)) (3.20)

    Suppose we have proven

    ϵjMH1(Ω)ω1ω1j(c2Δt2α1+cM1mu(x,tj)Hm(Ω))j=1,2,,n1. (3.21)

    And then we need prove ϵnMH1(Ω)ω1ω1n(c2Δt2α1+cM1mu(x,tn)Hm(Ω)) By Eq (3.19), Eq (3.21) we have

    ϵnMH1(Ω)[n1k=1ωkωk+1ω1+ωnω1]ω1ω1n[c2Δt2α1+cM1mu(x,tn)Hm(Ω)]. (3.22)

    Finally, by limkkΔt=T and the last equality have

    u(x,tn)unMH1(Ω)=ω1ω1nnα1Δtα1Tα1(c2Δt2α1+cM1muL(Hm)). (3.23)

    Then by Eq (3.23) and Lemma 2.8, we deduce that for sufficently large n, the assertion of the theorem is valid.

    Now we give error estimates for semi-discrete and full discrete problems in the following result.

    Theorem 3.4. Let {un}Nn=0Hm(Ω)H10(Ω),m>1 is solution of semi-discrete problem Eq (3.5), {unM}Nn=0 is the solution of full discrete problem Eq (3.15), and u0M=Pu0. Then we have

    ununMH1(Ω)cω1Tα1M1mmax1jnujHm(Ω)q1(1α1)Γ(2α1)+sj=2qjTα1αj(1αj)Γ(2αj),

    Proof. According to the definition of P, for solution un of semi-discrete, we have

    (Pun,νM)+1ω1(Punx,νMx)=n1k=1ωkωk+1ω1(unk,νM)+ωnω1(u0,νM)+pnω1(f,νM). (3.24)

    Let ˆenM=PununM,enM=ununM, by subtracting Eq (3.15) from Eq (3.24), we have

    (ˆenM,νM)+1ω1(ˆenMx,νMx)=n1k=1ωkωk+1ω1(enkM,νM)+ωnω1(e0M,νM). (3.25)

    Taking νM=ˆenM in above equation and using the triangular inequality enMH1(Ω)ˆenMH1(Ω)+unPunH1(Ω), we obtain

    enMH1(Ω)n1k=1ωkωk+1ω1enkMH1(Ω)+unPunH1(Ω)n1k=1ωkωk+1ω1ϵnkMH1(Ω)+cM1munHm(Ω). (3.26)

    Here, we use the mathematical induction similar to Theorem 3.3, and by limkkΔt=T. We have follow result

    ununMH1(Ω)=ω1ω1nnα1Δtα1Tα1cM1mmax1jnujHm(Ω). (3.27)

    Then by Eq (3.27) and Lemma 2.8, we deduce that for sufficently large n, the assertion of the theorem is valid.

    In this section, we discuss the uniqueness of solution and the ill-posed analysis of the unknown source identification problem.

    Denote the eigenvalues of the operator L as λn and the corresponding eigenfunctions as φnH2(Ω)H10(Ω), then we have Lφn=λnφn. Without loss of generality, suppose a family of eigenvalues 0<λ1λ2λn,limnλn=+, then {φn}n=1 constituting an orthonormal basis in space L2(Ω).

    Definition 4.1. ([43]) For any γ>0, define

    D((L)γ)={ψL2(Ω);n=1λ2γn|(ψ,φn)|2<},

    where (,) is the inner product in L2(Ω), and define its norm

    ψD((L)γ)={n=1λ2γn|(ψ,φn)|2}12.

    Lemma 4.2. ([34]) For fixed αA,qQ, let f(x)Ld(0,T),p(t)D((L)γ) with some d[1,] and γ[0,1]. Then Eq (1.1) exists a unique solution u(x,t) given by

    u(x,t)=n=1fnQn(t)φn(x), (4.1)

    where Qn(t)=t0p(s)(ts)α11E(n)α,α1(ts)ds, and fn=(f,φn).

    From Lemma 4.2, we know there exists a unique weak solution uL2((0,T);H2(Ω)H10(Ω)) for the direct problem Eq (1.1), if we know a source function f(x)p(t)L(0,T;L2(Ω)). And we know the formal solution for Eq (1.1) can be expressed by

    u(x,t)=n=1fnt0p(s)(ts)α11E(n)α,α1(ts)dsφn(x). (4.2)

    Applying u(x,T)=g(x), we have

    n=1gnφn(x)=g(x)=n=1fnQn(T)φn(x), (4.3)

    where gn=(g,φn)andQn(T)=t0p(s)(Ts)α11E(n)α,α1(Ts)ds.

    Theorem 4.3. If p(t)C[0,T] satisfying p(t)p0>0 for all t[0,T], then the solution u(x,t),f(x) of problem Eq (1.1) is unique.

    Proof. By Lemma 2.3, we have E(n)α,α1(t)>0, for t>0. From Lemma 2.2 and Proposition 2.5, we know

    Qn(T)p0T0(Ts)α11E(n)α,α1(Ts)ds=p0Tα1E(n)α,1+α1(T)>0. (4.4)

    Hence, we know if g(x)=0, we have f(x)=0, and by Eq (4.2) we know u(x,t)=0. The proof is completed.

    In the following, we show the inverse source problem Eq (1.1) is ill-posed. Without loss of generality, we take a final data gk(x)=φk(x)λk in Eq (1.1) and then know gk=1λk0,ask. The corresponding source terms are fk(x)=φk(x)Qk(T)λk, and we have fk=1Qk(T)λk. By Lemma 2.7, we have

    Qk(T)pC[0,T]Tα1E(k)α,1+α1(T)pC[0,T]λk. (4.5)

    Hence fkλkqC[0,T], as k. Then we know the inverse source problem Eq (1.1) is ill-posed.

    For the sake of argument, we just talk about the one-dimensional case.

    First, let LM(x) denotes the Legendre polynomial of degree M. ξj,j=1,2,,M, are the Legendre-Gauss-Lobatto points, and these points satisfy (1x2)LM(x)=0. we consider the follow equation

    {(unM,νM)+1ω1(aunMx,νMx)1ω1(cunM,νM)=n1k=1ωkωk+1ω1(unkM,νM)+ωnω1(u0M,νM)+pnω1(fM,νM),(uNM,νM)=(gM,νM). (5.1)

    We set

    unM(x)=uM(x,nΔt)=Mi=0δnihi(x), (5.2)

    Here, hi is the Lagrangian polynomial defined in Ω i.e. hiVM,hi(ξj)=εij, with εij: the Krinecker-delta symbol. Since first equation in Eq (1.1) satisfies Dirichlet boundary condition, therefore, we have δn0=δnM=0. Take νM=hj(x),j=1,2,,M1, we can obtain

    {(M1i=1δnihi,hj)+1ω1(aM1i=1δnihix,hjx)1ω1(cM1i=1δnihi,hj)=n1k=1ωkωk+1ω1(M1i=1δnkihi,hj)+ωnω1(M1i=1δ0ihi,hj)+pnω1(fM,hj)(M1i=1δNihi,hj)=(gM,hj).

    Denote

    F=((fM,h1),,(fM,hM1))T,G=((gM,h1),,(gM,hM1))T,

    and set

    δn=(δn1,,δnM1)T,

    By above equation, we have

    {Kδn=1ω1pnF+n1k=1ωkωk+1ω1ˆKδnkˆKδN=G, (5.3)

    where i,j=1,2,,M1, and

    (K)ij=(hi,hj)+1ω1(ahix,hjx)1ω1(chi,hj),(ˆK)ij=(hi,hj). (5.4)

    Here, we use the Legendre Gauss-type quadratures to approximate Eq (5.4), and by the first equation in Eq (5.3) and Eq (5.2), we can solve the forward problem.

    For the inverse problem, by Eq (5.3), we have

    {δ1=K11ω1p1F=A1Fn=1δ2=K1(1ω1p2F+(1ω2ω1)ˆKδ1)=A2Fn=2δN1=AN1Fn=N1. (5.5)

    For n=N, we have

    {KδN=1ω1pNF+N1k=1ωkωk+1ω1ˆKδNkˆKδN=G, (5.6)

    By Eq (5.5), Eq (5.6), we have

    AF=G. (5.7)

    And then by solving linear equation Eq (5.7), we can obtain F. Further, we recover fM as an approximation of the source term by using a suitable integration rule.

    In order to prove the ill-posedness of equation Eq (5.7), an upper bound of A12 in Eq (5.7) is given below.

    Theorem 5.1. Let the orthogonal basis function {hj}M1j=1 satisfy hj(x)H1(Ω). Then the inverse matrix of the coefficient matrix A in Eq (5.7) have

    A12C(M1)2max1jM1hj2H1(Ω)min1jM1hj2H1(Ω) (5.8)

    Proof. From the first equation of Eq (5.6), we can get

    δN=pNω1K1F+N1k=1ωkωk+1ω1K1ˆKδNk. (5.9)

    Combining the second equation of Eq (5.6), we have

    pNω1ˆKK1F+N1k=1ωkωk+1ω1ˆKK1ˆKδNk=G. (5.10)

    Let k=Nk(alsodenotek), we have

    pNω1ˆKK1F+N1k=1ωNkωNk+1ω1ˆKK1ˆKδk=G. (5.11)

    On the other hand, by Eq (5.5) we have

    δk=(pkω1K1+k1m=1ωmωm+1ω1K1ˆKAkm)F=:AkF,k=1,2,,N1. (5.12)

    By substituting Eq (5.12) into Eq (5.11), the equation AF=G is obtained, where

    A=pNω1ˆKK1+N1k=1ωNkωNk+1ω1ˆKK1ˆKAk. (5.13)

    Since ωk is non-negative and decreasing with respect to k, therefore we have

    ωNkωNk+1ω1<1. (5.14)

    The definition of Ak tells us that the second part of A is the higher order term of the first part. Therefore, we only need to estimate the matrix

    A˜A:=pNω1ˆKK1. (5.15)

    Due to F˜A1G, we only have to estimate the upper bound of ˜A12. According to Eq (5.4), ˆK is a diagonal matrix, and

    (KˆK1)ij=(hi,hj)hj2+1ω1hj2(ahix,hjx)1ω1hj2(chi,hj). (5.16)

    Therefore, we have

    ˜A12=ω1pNKˆK12ω1pNKˆK1F=ω1pN(M1i,j=1(KˆK1)2ij)12=ω1pN(M1i,j=11hj4((hi,hj)hj2+1ω1(ahix,hjx)1ω1(chi,hj))2)123ω1min1jM1hj2H1(Ω)pN(M1i,j=1(hi,hj)2+1ω21(ahix,hjx)2+1ω21(chi,hj)2)12C(M1)2max1jM1hj2H1(Ω)min1jM1hj2H1(Ω), (5.17)

    where C is dependent on a,c,M,α,qj.

    Remark 1. It can be seen from Theorem 5.1 that A12 is completely determined by the projection dimension M1 of the spectral method and the maximum and minimum norm of the basis function hj(x). It can be seen from Eq (5.8) that as the projection dimension increases, the numerator increases with respect to j, while the denominator decreases with respect to j. Therefore, the projection dimension M1 can be used as regularization parameter. When the appropriate M is selected, the minimum value of A12 can be guaranteed, thus reducing the ill-posedness of the inverse problem.

    Remark 2. Here we emphasize that the Galerkin spectral method is a projection method. In the projection method, as long as the appropriate dimension of the projection space is selected, the regularization effect can be effectively generated for the linear ill-posed problem, and it is no longer necessary to adopt other regularization techniques for the problem. This phenomenon is sometimes called self-regularization or regularization by projection. A series of descriptions and proofs of this phenomenon can be found in the literature [44].

    In order to compare the finite difference method the and Galerkin spectral method in application of the multi-term time-fractional diffusion equation. We derive follow the scheme similar to reference [45].

    In the finite difference algorithm, we can denote uniu(xi,tn), where xi=iΔx,i=1,...,C. and tn=nΔt,n=0,1,...,N.

    The space has the following discrete form:

    Lu(xi,tn)1(Δx)2(ai+12uni+1)((ai+12+ai12)uni+ai12uni1)+c(xi)uni,

    for i=1,...,C1,n=1,...,N where ai+12=a(xi+12) with xi+12=(xi+xi+1)/2. In problem Eq (1.1), according to initial condition and boundary condition, we can get a numerical solution for forward problem Eq (1.1) from the finite difference scheme

    nk=1ωk(unk+1unk)=1(Δx)2(ai+12uni+1)((ai+12+ai12)uni+ai12uni1)+c(xi)uni+f(xi)p(tn). (5.18)

    Denote Un=(un1,un2,...,unC1)T,Y=(f(x1),f(x2),...,f(xC1))T, then the scheme Eq (5.18) leads to the following iterative scheme

    BU1=U0+1ω1p(t1)Y,BUn=c1Un1+c2Un2+...+cn1U1+ωnω1U0+1ω1p(tn)Y, (5.19)

    where cn=(ωnωn+1)/ω1 and B is a tridiagonal matrix given by Bii=ai+12+ai12(Δx)2ω1+1c(xi)ω1 for i=1,2,,C1 and Bi,i1=ai12(Δx)2ω1 for i=2,3,,C1 and Bi,i+1=ai+12(Δx)2ω1 for i=1,3,,C2.

    Inverse source problem based on finite difference. Similar to Section 5.1, we can derive a linear equation.

    ZY=D, (5.20)

    where D=(g(x1),g(x2),...,g(xC1))T, and Z is a matrix.

    In this section, first, we verify the stability and validity of the proposed numerical methods. Without lose of generality, let the maximum time is T=1. The noisy data are generated by adding random perturbations, i.e,

    gδ(x)=g(x)+εg(x)(2rand(size(g(x))1), (6.1)

    where ε is relative noise level and rand() generate random numbers uniformly distributed on [0,1]. The corresponding noise level is calculated by δ=ggδ. To show the accuracy of numerical solution, we compute the approximate L2-norm error denoted by

    e(f,ε)=ffδ, (6.2)

    and the approximate L2-norm relative error as

    er(f,ε)=ffδ/f, (6.3)

    where fδ is term reconstructed and f is the exact solution.

    In addition, for the convenience of writing, we will abbreviate the Galerkin spectral method as GSM and the finite difference method as FDM.

    Example 1. In this example, we firstly consider a three-term time fractional diffusion equation with a exact solution. Let ΩT=(1,1)×(0,T) and T=1,qj1. Take a source function p(t)=2Γ(3α1)t2α1+2Γ(3α2)t2α2+2Γ(3α3)t2α3+4π2t2, and f(x)=sin(2πx). Further more, we have exact analytical solution

    u(x,t)=t2sin(2πx). (6.4)

    In the first three graphs of Figure 1, namely (a), (b) and (c), we show the errors in L2-norm and L-norm for different time step in case of α1=(0.9,0.8,0.7),α2=(0.6,0.5,0.3), and α3=(0.3,0.2,0.1). Here fix the polynomial degree M=24. In the same case, the graphs (d), (e) and (f) from Figure 1 show the the errors in L2-norm and L-norm for different polynomial degree M. Here we also take time step Δt=1e4. From Figure 1, we can find that the theoretical convergence accuracy (see Theorem 3.3) is in good agreement with the numerical results.

    Figure 1.  The result of the forward problem for Example 1.

    In Table 1 and Figure 2, we make a comparison with FDM in terms of the forward problem by using the scheme Eq (5.1). Clearly, we observe that the scheme Eq (5.1) has higher accuracy than FDM on the forward problem. As a further step, we choose Δt=1e4 to obtain satisfactory results in time and accuracy.

    Table 1.  Error of two approximation methods of Example 1 for α=(0.3,0.2,0.1).
    timestep approximation method L2 error L error run time
    Δt=1e1 GSM 0.002121690 2.121690e-04 88.254s
    FDM 0.010039748 0.001419834 1.341s
    Δt=1e2 GSM 4.382322e-05 4.382322e-06 90.043s
    FDM 0.0085668512 0.001211535 10.076s
    Δt=1e3 GSM 8.535546e-07 8.535546e-08 90.622s
    FDM 0.008536392 0.001207228 96.931s
    Δt=1e4 GSM 1.913459e-08 1.913459e-09 198.829s
    FDM 0.008535798 0.001207144 1348.549s
    Δt=1e5 GSM 1.007804e-08 1.007804e-9 10834.238s
    FDM 0.008535787 0.001207142 33733.606s

     | Show Table
    DownLoad: CSV
    Figure 2.  Error of two approximation methods of Example 1 for α=(0.3,0.2,0.1).

    On the inverse source problem aspect. In order to avoid the impact of 'inverse crime'. We take M=20,Δt=1e4 to solve the direct problem by using the GSM Eq (5.1), and take M=20,Δt=1e3 to solve the inverse source problem by the scheme Eq (5.7). The related results are shown in (e), (f) of Figure 3 and Table 2. For the inverse source problem based on FDM. We take C=100,Δt=5e3. Numerical result for α2=(0.3,0.2,0.1) with various noise levels ε=5%,1%,0.5%,0.1% are shown in (c), (d) of Figure 3.

    Figure 3.  Result of Example 1 with α=(0.3,0.2,0.1).
    Table 2.  Relative error of the GSM for Example 1 for different α with ε=5%,1%,0.5%,0.1%.
    ϵ α er(f,ε)
    ε=5% (0.3,0.2,0.1) 0.082339
    (0.9,0.8,0.7) 0.074294
    (0.6,0.5,0.3) 0.059221
    ε=1% (0.3,0.2,0.1) 0.030306
    (0.9,0.6,0.3) 0.034862
    (0.6,0.5,0.3) 0.049156
    ε=0.5% (0.3,0.2,0.1) 0.014814
    (0.9,0.8,0.7) 0.043858
    (0.6,0.5,0.3) 0.043187
    ε=0.1% (0.3,0.2,0.1) 0.013495
    (0.9,0.8,0.7) 0.010624
    (0.6,0.5,0.3) 0.013953

     | Show Table
    DownLoad: CSV

    Example 2. In this example, let ΩT=(0,1)×(0,T),s=3,qj1,a(x)=1,c(x)=0. Take a exact source function p(t)=1 and f(x)=x(x0.4)(x0.6)(x0.8)(x1). We first solve the direct problem by using the GSM Eq (5.1) to obtain the additional data g(x). Then we use the scheme Eq (5.7) to solve inverse source problem. Numerical result for α1=(0.9,0.8,0.7),α2=(0.3,0.2,0.1), and α3=(0.6,0.5,0.3), with various noise levels ε=1%,0.5%,0.1%,0.01% are presented in (e), (f) of Figure 4 and Table 3. For the inverse source problem based on FDM. We take C=100,Δt=5e3. Numerical result for α2=(0.3,0.2,0.1) with various noise levels ε=1%,0.5%,0.1%,0.01% are shown in (c), (d) of Figure 4.

    Figure 4.  Result of Example 2 with α=(0.3,0.2,0.1).
    Table 3.  Relative error of the GSM for Example 2 for different α with ε=1%,0.5%,0.1%,0.01%.
    ε α er(f,ε)
    ε=1% (0.3,0.2,0.1) 0.085368
    (0.9,0.8,0.7) 0.092676
    (0.6,0.5,0.3) 0.088179
    ε=0.5% (0.3,0.2,0.1) 0.076024
    (0.9,0.8,0.7) 0.082049
    (0.6,0.5,0.3) 0.079423
    ε=0.1% (0.3,0.2,0.1) 0.073424
    (0.9,0.8,0.7) 0.0800324
    (0.6,0.5,0.3) 0.077136
    ε=0.01% (0.3,0.2,0.1) 0.0652520
    (0.9,0.8,0.7) 0.066370
    (0.6,0.5,0.3) 0.066352

     | Show Table
    DownLoad: CSV

    Example 3. In this example, let ΩT=(0,1)×(0,T),s=3,qj1,a(x)=1,c(x)=0. We consider a continuous piecewise smooth function.

    f(x)={2x,0x0.5,2x+2,0.5x<1. (6.5)

    We first solve the direct problem by using the GSM Eq (5.1) to obtain the additional data g(x). Then we use the scheme Eq (5.7) to solve inverse source problem. Numerical result for α1=(0.9,0.8,0.7),α2=(0.3,0.2,0.1), and α3=(0.6,0.5,0.3), with various noise levels ε=1%,0.5%,0.1%,0.01% are presented in (e), (f) of Figure 5 and Table 4. For the inverse source problem based on FDM. We take C=100,Δt=5e3. Numerical result for α2=(0.3,0.2,0.1) with various noise levels ε=1%,0.5%,0.1%,0.01% are shown in (c), (d) of Figure 5.

    Figure 5.  Result of Example 3 with α=(0.3,0.2,0.1).
    Table 4.  Relative error of the GSM for Example 3 for different α with ε=1%,0.5%,0.1%,0.01%.
    ε α er(f,ε)
    ε=1% (0.3,0.2,0.1) 0.045909
    (0.9,0.8,0.7) 0.058139
    (0.6,0.5,0.3) 0.058454
    ε=0.5% (0.3,0.2,0.1) 0.042321
    (0.9,0.8,0.7) 0.048615
    (0.6,0.5,0.3) 0.047175
    ε=0.1% (0.3,0.2,0.1) 0.040261
    (0.9,0.8,0.7) 0.048302
    (0.6,0.5,0.3) 0.046372
    ε=0.01% (0.3,0.2,0.1) 0.012787
    (0.9,0.8,0.7) 0.011516
    (0.6,0.5,0.3) 0.011305

     | Show Table
    DownLoad: CSV

    Example 4. Let ΩT=(0,1)×(0,T),s=3,qj1,a(x)=x2+1,c(x)=1. Take a exact source function p(t)=1 and f(x)=exsin(7πx). We first solve the direct problem by using the GSM Eq (5.1) to obtain the additional data g(x). Then we use the scheme Eq (5.7) to solve inverse source problem. Numerical result for α1=(0.9,0.8,0.7),α2=(0.3,0.2,0.1), and α3=(0.6,0.5,0.3), with various noise levels ε=10%,5%,1%,0.1% are presented in (e), (f) of Figure 6 and Table 5. For the inverse source problem based on FDM. We take C=100,Δt=5e3. Numerical result for α2=(0.3,0.2,0.1) with various noise levels ε=1%,0.5%,0.1%,0.01% are shown in (c), (d) of Figure. 6.

    Figure 6.  Result of Example 4 with α=(0.3,0.2,0.1).
    Table 5.  Relative error of the GSM for Example 4 for different α with ε=10%,5%,1%,0.1%.
    ε α er(f,ε)
    ε=10% (0.3,0.2,0.1) 0.100485
    (0.9,0.8,0.7) 0.111043
    (0.6,0.5,0.3) 0.107902
    ε=5% (0.3,0.2,0.1) 0.093797
    (0.9,0.8,0.7) 0.068835
    (0.6,0.5,0.3) 0.064510
    ε=1% (0.3,0.2,0.1) 0.072740
    (0.9,0.8,0.7) 0.073513
    (0.6,0.5,0.3) 0.072786
    ε=0.1% (0.3,0.2,0.1) 0.034704
    (0.9,0.8,0.7) 0.037078
    (0.6,0.5,0.3) 0.035328

     | Show Table
    DownLoad: CSV

    By observing Figures 36, we can see that (a) denotes the graph of the numerical solution of the forward problem, namely the input data. (b) is the numerical approximation of the inverse source problem without noisy data. (c) and (d) are results of the inverse source problem based on the FDM with and without regularization method, and (e) and (f) are results of the inverse source problem based on the GSM under the same case.

    From Tables 69, we show the error levels of two numerical reconstruction method with and without regularization method. We can find that the numerical results of the inverse problem based on FDM with regularization method are obviously better than those without regularization method. However, the numerical results of the inverse problem based on the GSM with and without regularization method do not show a big difference compared with the numerical results of the FDM. There may even be cases where regularization method is not required to get a better result. Hence, we can conclude that scheme Eq (5.7) itself has regularization function compared to the FDM in the inverse source problem. Next, we explain why scheme Eq (5.7) is not sensitive to noise level.

    Table 6.  Comparison of relative error of numerical results for Example 1 for α=(0.3,0.2,0.1) with ε=5%,1%,0.5%,0.1%.
    method ε=5% ε=1% ε=0.5% ε=0.1%
    GSM (no regularization) 0.2501 0.0832 0.0290 0.0079
    GSM (regularization) 0.0823 0.0303 0.0134 0.0148
    FDM (no regularization) 4.3792 0.7871 0.4036 0.0965
    FDM (regularization) 0.0301 0.0052 0.0025 6.273e-4

     | Show Table
    DownLoad: CSV
    Table 7.  Comparison of relative error of numerical results for Example 2 for α=(0.3,0.2,0.1) with ε=1%,0.5%,0.1%,0.01%.
    method ε=1% ε=0.5% ε=0.1% ε=0.01%
    GSM (no regularization) 0.0944 0.0635 0.0119 0.0010
    GSM (regularization) 0.0853 0.0760 0.0734 0.0652
    FDM (no regularization) 2.6407 1.220 0.2478 0.027
    FDM (regularization) 0.2043 0.0520 0.0208 0.0164

     | Show Table
    DownLoad: CSV
    Table 8.  Comparison of relative error of numerical results for Example 3 for α=(0.3,0.2,0.1) with ε=1%,0.5%,0.1%,0.01%.
    method ε=1% ε=0.5% ε=0.1% ε=0.01%
    GSM (no regularization) 0.4327 0.2868 0.0365 0.0126
    GSM (regularization) 0.0459 0.0423 0.0402 0.0127
    FDM (no regularization) 11.808 5.0791 1.4387 0.1025
    FDM (regularization) 0.0353 0.0239 0.0116 0.0102

     | Show Table
    DownLoad: CSV
    Table 9.  Comparison of relative error of numerical results for Example 4 for α=(0.3,0.2,0.1)with ε=10%,5%,1%,0.1%.
    method ε=10% ε=5% ε=1% ε=0.1%
    GSM (no regularization) 0.2924 0.0878 0.0425 0.0233
    GSM (regularization) 0.1004 0.0937 0.0727 0.0347
    FDM (no regularization) 2.79904 1.1476 0.2951 0.1435
    FDM (regularization) 0.1435 0.0846 0.0181 0.0014

     | Show Table
    DownLoad: CSV

    In Table 10, we show norm of matrices G, A1 and A1G for Example 1. Obviously, A1G2 is small and A12 has a reasonable size. Therefore numerical solutions are not sensitive with respect to the perturbation in the initial data. But, In Table 11, The numerical format of FDM is sensitive to the perturbation in the initial datal. Therefore, GSM has better anti-interference than FDM. On the other hand, in Tables 12 and 13, We can find that although GSM has better anti-interference than FMD, it can not completely get rid of the regularization method when the unknown source has less regularity and the meased data has a high noise level.

    Table 10.  GSM of ill-posedness for Example 1 for different ε.
    ε G2 A12 A1G2
    0.05 0.356708 86.11457 0.369712
    0.01 0.360401 86.11457 0.360933
    0.005 0.360111 86.11457 0.360085
    0.001 0.360341 86.11457 0.360345
    0.0001 0.360348 86.11457 0.360348

     | Show Table
    DownLoad: CSV
    Table 11.  FDM of ill-posedness for Example 1 for different ε.
    ε D2 Z12 Z1D2
    0.05 7.0881 2.322e+02 27.8227
    0.01 7.0743 2.322e+02 8.8667
    0.005 7.0809 2.322e+02 7.5734
    0.001 7.0804 2.322e+02 7.0929
    0.0001 7.0796 2.322e+02 7.0712

     | Show Table
    DownLoad: CSV
    Table 12.  GSM of ill-posedness for Example 3 for different ε.
    ε G2 A12 A1G2
    0.01 0.013187 1.5104e+04 0.161897
    0.005 0.013226 1.5104e+04 0.164711
    0.001 0.013227 1.5104e+04 0.159407
    0.0001 0.013228 1.5104e+04 0.159405

     | Show Table
    DownLoad: CSV
    Table 13.  FDM of ill-posedness for Example 3 for different ε.
    ε D2 Z12 Z1D2
    0.01 0.459280 3.9992e+04 76.739424
    0.005 0.459857 3.9992e+04 27.761398
    0.001 0.459380 3.9992e+04 7.957420
    0.0001 0.459332 3.9992e+04 5.798087

     | Show Table
    DownLoad: CSV

    In this paper, we first obtain a high accuracy numerical solution by using the GSM, and give the error estimates between exact solution and semi-discrete solution as well as full-discrete one, and compare with the FDM. It is indicate that our method has a better accuracy. Secondly, the GSM is extended to solve the inverse source problem. Moreover, we find that this method can effectively reduce the ill-posednes of inverse source problem compared with the traditional FDM. Thus the spectral method itself can play a regularization role. It should be mentioned that the estimates given in the paper are also valid in two and three dimensional cases. In the following work, we will continue to try to optimize the GSM numerical scheme so that it can better play the role of regularization method. For example, replacing the original basis functions with smooth periodic functions or equidistant trigonometric functions to avoid Runge phenomenon.

    This work is supported by the NSF of China (grant no. 12201502), the Youth Science and Technology Fund of Gansu Province (grant no. 20JR10RA099), the Innovation Capacity Improvement Project for Colleges and Universities of Gansu Province (grant no. 2020B-088) and the Innovation star of Gansu Province (grant no. 2022CXZX-324).

    The authors declare there is no conflict of interest.



    [1] R. Metzler, J. Klafter, The random walk's guide to anomalous diffusion: a fractional dynamics approach, Phys. Rev. E, 61 (2000), 6308–6311. https://doi.org/10.1016/s0370-1573(00)00070-3 doi: 10.1016/s0370-1573(00)00070-3
    [2] B. Berkowitz, H. Scher, S. Silliman Anomalous transport in laboratory-scale, heterogeneous porous media, Water Resour. Res., 36 (2000), 149–158. https://doi.org/10.1029/2000wr900026 doi: 10.1029/2000wr900026
    [3] F. Mainardi, Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models, Singapore: World Scientific, 2010. https://doi.org/10.1142/p614
    [4] Y. Hatano, N. Hatano, Dispersive transport of ions in colum experiments: an explanation of long-tailed profiles, Water Resour. Res., 34 (1998), 1027–1033. https://doi.org/10.1029/98WR00214 doi: 10.1029/98WR00214
    [5] S. Benson, M. B. Meerschaert, Fracral mobile/immobile solute transport, Water Resour. Res., 39 (2003), 1–12. https://doi.org/2003WR002141
    [6] S. Rina, D. A. Benson, M. M. Mark, B. Boris, Fractal mobile/immobile solute transport, Water Resour. Res., 39 (2003), 1296. https://doi.org/10.1029/2003WR002141
    [7] R. R. Nigmatullin, The realization of the generalized transfer equation in a medium with fractal geometry, Phys. Stat. Sol. B, 133 (1986), 425–430. https://doi.org/10.1002/pssb.2221330150 doi: 10.1002/pssb.2221330150
    [8] A. H. Bhrawy, M. A. Zaky, Highly accurate numerical schemes for multi-dimensional space variable-order fractional Schrödinger equations, Comput. Math. Appl., 73 (2017), 1100–1117. https://doi.org/10.1016/j.camwa.2016.11.019 doi: 10.1016/j.camwa.2016.11.019
    [9] Z. Y. Li, Y. K. Liu, M. Yamamoto, Initial-boundary value problems for multi-term time-fractional diffusion equations with positive constant coefficients, Appl. Math. Comput., 257 (2015), 381–397. https://doi.org/10.1016/j.amc.2014.11.073 doi: 10.1016/j.amc.2014.11.073
    [10] Y. Luchko, Initial-boundary-value problems for the generalized multi-term time-fractional diffusion equation, J. Math. Anal. Appl., 374 (2011), 538–548. https://doi.org/10.1016/j.jmaa.2010.08.048 doi: 10.1016/j.jmaa.2010.08.048
    [11] X. L. Ding, J. J. Nieto, Analytical solutions for multi-term time-space fractional partial differential equations with nonlocal damping terms, Fract. Calc. Appl. Anal., 21 (2018), 312–335. https://doi.org/10.1016/j.cnsns.2018.05.022 doi: 10.1016/j.cnsns.2018.05.022
    [12] C. S. Sin, G. I. Ri, M. C. Kim, Analytical solutions to multi-term time-space Caputo-Riesz fractional diffusion equations on an infinite domain, Adv. Difference Equ., 1 (2017), 306. https://doi.org/10.1186/s13662-017-1369-x doi: 10.1186/s13662-017-1369-x
    [13] G. S. Li, C. L. Sun, X. Z. Jia, D. H. Du, Numerical solution to the multi-term time fractional diffusion equation in a finite domain, Numer. Math. Theory Methods Appl., 9 (2016), 337–357. https://doi.org/10.4208/nmtma.2016.y13024 doi: 10.4208/nmtma.2016.y13024
    [14] M. R. Cui, Finite difference schemes for the two-dimensional multi-term time-fractional diffusion equations with variable coefficients, Comput. Appl. Math., 40 (2021), 167. https://doi.org/10.1007/s40314-021-01551-11 doi: 10.1007/s40314-021-01551-11
    [15] Y. L. Zhao, P. Zhu, X. M. Gu, A second-order accurate implicit difference scheme for time fractional reaction-diffusion equation with variable coefficients and time drift term, East Asian J. Appl. Math., 9 (2019), 723–754. https://doi.org/10.4208/eajam.200618.250319 doi: 10.4208/eajam.200618.250319
    [16] Z. Y. Li, O. Y. Imanuvilov, M. Yamamoto, Uniqueness in inverse boundary value problems for fractional diffusion equations, Inverse Problems, 32 (2016), 015004. https://doi.org/10.1088/0266-5611/32/1/015004 doi: 10.1088/0266-5611/32/1/015004
    [17] W. P. Bu, S. Shu, X. Q. Yue, A. G. Xiao, W. Zeng, Space-time finite element method for the multi-term time-space fractional diffusion equation on a two-dimensional domain, Comput. Math. Appl., 75 (2019), 1367–1379. https://doi.org/10.1016/j.camwa.2018.11.033 doi: 10.1016/j.camwa.2018.11.033
    [18] J. Zhou, D. Xu, H. B. Chen, A weak Galerkin finite element method for multi-term time-fractional diffusion equations, East Asian J. Appl. Math., 8 (2018), 181–193. https://doi.org/10.4208/eajam.260617.151117a doi: 10.4208/eajam.260617.151117a
    [19] L. L. Wei, Stability and convergence of a fully discrete local discontinuous Galerkin method for multi-term time fractional diffusion equations, Numer. Algorithms, 76 (2017), 695–707. https://doi.org/10.1007/s11075-017-0277-1 doi: 10.1007/s11075-017-0277-1
    [20] S. M. Guo, L. Q. Mei, Z. Q. Zhang, Y. T. Jiang, Finite difference/spectral-Galerkin method for a two-dimensional distributed-order time-space fractional reaction-diffusion equation, Appl. Math. Lett., 85 (2018), 157–163. https://doi.org/10.1016/j.aml.2018.06.005 doi: 10.1016/j.aml.2018.06.005
    [21] R. M. Zheng, F. W. Liu, X. Y. Jiang, A Legendre spectral method on graded meshes for the two-dimensional multi-term time-fractional diffusion equation with non-smooth solutions, Appl. Math. Lett., 104 (2020), 106247. https://doi.org/10.1016/j.aml.2020.106247 doi: 10.1016/j.aml.2020.106247
    [22] Y. Q. Liu, X. L. Yin, F. W. Liu, X. Y. Xin, Y. F. Shen, L. B. Feng, An alternating direction implicit Legendre spectral method for simulating a 2D multi-term time-fractional Oldroyd-B fluid type diffusion equation, Comput. Math. Appl., 113 (2022), 160–173. https://doi.org/10.1016/j.camwa.2022.03.020 doi: 10.1016/j.camwa.2022.03.020
    [23] M. Zheng, F. Liu, V. Anh, I. Turner, A high-order spectral method for the multi-term time-fractional diffusion equations, Appl. Math. Model., 40 (2016), 4970–4985. https://doi.org/10.1016/j.apm.2015.12.011 doi: 10.1016/j.apm.2015.12.011
    [24] M. A. Zaky, A Legendre spectral quadrature tau method for the multi-term time-fractional diffusion equations, Comput. Appl. Math., 37 (2018), 3525–3538. https://doi.org/10.1007/S40314-017-0530-1 doi: 10.1007/S40314-017-0530-1
    [25] Y. Zhang, X. Xu, Inverse source problem for a fractional diffusion equation, Inverse Probl, 27 (2011), 538–548. https://doi.org/10.1088/0266-5611/27/3/035010 doi: 10.1088/0266-5611/27/3/035010
    [26] T. Wei, L. L. Sun, Y. S. Li, Uniqueness for an inverse space-dependent source term in a multi-dimensional time-fractional diffusion equation, Appl. Math. Lett., 61 (2016), 108–113. https://doi.org/10.1016/j.aml.2016.05.004 doi: 10.1016/j.aml.2016.05.004
    [27] X. B. Yan, T. Wei, Determine a space-dependent source term in a time fractional diffusion-wave equation, Acta Appl. Math., 165 (2020), 163–181. https://doi.org/10.1007/s10440-019-00248-2 doi: 10.1007/s10440-019-00248-2
    [28] L. L. Sun, X. B. Yan, K. Liao, Simultaneous inversion of a fractional order and a space source term in an anomalous diffusion model, J Inverse Ill Posed Probl, 30 (2022), 791–805. https://doi.org/10.1515/jiip-2021-0027 doi: 10.1515/jiip-2021-0027
    [29] S. Yeganeh, R. Mokhtari, J. S. Hesthaven, Space-dependent source determination in a time-fractional diffusion equation using a local discontinuous Galerkin method, BIT Numer. Math., 57 (2017), 685–707. https://doi.org/10.1007/s10543-017-0648-y doi: 10.1007/s10543-017-0648-y
    [30] D. J. Jiang, Z. Y. Li, Y. K. Liu, M. Yamamoto, Weak unique continuation property and a related inverse source problem for time-fractional diffusion-advection equations, Inverse Problems, 33 (2017), 055013. https://doi.org/10.1088/1361-6420/aa58d1 doi: 10.1088/1361-6420/aa58d1
    [31] Y. S. Li, L. L. Sun, Z. Q. Zhang, T. Wei, Identification of the time-dependent source term in a multi-term time-fractional diffusion equation, Numer. Algorithms, 82 (2019), 1279–1301. https://doi.org/10.1007/s11075-019-00654-5 doi: 10.1007/s11075-019-00654-5
    [32] L. L. Sun, X. B. Yan, An inverse source problem for distributed order time-fractional diffusion equation, Adv. Math. Phys., 2020 (2020), 1825235. https://doi.org/10.1155/2020/1825235 doi: 10.1155/2020/1825235
    [33] S. A. Malik, A. Ilyas, A. Samreen, Simultaneous determination of a source term and diffusion concentration for a multi-term space-time fractional diffusion equation, Math. Model. Anal., 26 (2021), 411–431. https://doi.org/10.3846/mma.2021.11911 doi: 10.3846/mma.2021.11911
    [34] L. L. Sun, Y. S. Li, Y. Zhang, Simultaneous inversion of the potential term and the fractional orders in a multi-term time-fractional diffusion equation, Inverse Probl, 37 (2021), 055007. https://doi.org/10.1088/1361-6420/abf162 doi: 10.1088/1361-6420/abf162
    [35] Y. H. Lin, H. Y. Liu, X. Liu, S. Zhang, Simultaneous recoveries for semilinear parabolic systems, Inverse Probl, 38 (2022), 115006. https://doi.org/10.1088/1361-6420/ac91ee doi: 10.1088/1361-6420/ac91ee
    [36] H. Y. Liu, G. Uhlmann, Determining both sound speed and internal source in thermo-and photo-acoustic tomography, Inverse Probl, 31 (2015), 105005. https://doi.org/10.1088/0266-5611/31/10/105005 doi: 10.1088/0266-5611/31/10/105005
    [37] X. L. Cao, H. Y. Liu, Determining a fractional Helmholtz equation with unknown source and scattering potential, Commun. Math. Sci., 17 (2019), 1861–1876. https://doi.org/110.4310/CMS.2019.v17.n7.a5
    [38] L. L. Sun, Y. Zhang, T. Wei, Recovering the time-dependent potential function in a multi-term time-fractional diffusion equation, Appl. Numer. Math., 135 (2019), 228–245. https://doi.org/10.1016/j.apnum.2018.09.001 doi: 10.1016/j.apnum.2018.09.001
    [39] E. Bazhlekova, Properties of the fundamental and the impulse-response solutions of multi-term fractional differential equations, Complex Analysis and Applications, 2 (2013), 55–64.
    [40] C. L. Sun, G. S. Li, X. Z. Jia, Numerical inversion for the initial distribution in the multi-term time-fractional diffusion equation using final observations, Adv. Appl. Math. Mech., 9 (2017), 1525–1546. https://doi.org/10.4208/aamm.OA-2016-0170 doi: 10.4208/aamm.OA-2016-0170
    [41] Y. M. Lin, C. J. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, J. Comput. Phys., 225 (2007), 1533–1552. https://doi.org/10.1016/j.jcp.2007.02.001 doi: 10.1016/j.jcp.2007.02.001
    [42] C. Bernardi, Y. Maday, Approximations spectrales de problèmes Aux Limites Elliptiques, Paris: Springer-Verlag, 1992.
    [43] K. Sakamoto, M. Yamamoto, Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems, J. Math. Anal. Appl., 382 (2011), 538–548. https://doi.org/10.1016/j.jmaa.2011.04.058 doi: 10.1016/j.jmaa.2011.04.058
    [44] N. Frank, Regularisierung schlecht gestellter Probleme durch Projektionsverfahren, Numer. Math., 28 (1977), 329–341. https://doi.org/10.1007/BF01389972 doi: 10.1007/BF01389972
    [45] T. Wei, J. G. Wang, A modified quasi-boundary value method for the backward time-fractional diffusion problem, ESAIM Math. Model. Numer. Anal., 48 (2014), 603–621. https://doi.org/10.1051/m2an/2013107 doi: 10.1051/m2an/2013107
  • This article has been cited by:

    1. Liangliang Sun, Zhaoqi Zhang, Exponential Tikhonov Regularization Method for an Inverse Source Problem in a Sub-diffusion Equation, 2024, -1, 1027-5487, 10.11650/tjm/240901
    2. Liangliang Sun, Yuxin Wang, Maoli Chang, A Fractional-order Quasi-reversibility Method to a Backward Problem for the Multi-term Time-fractional Diffusion Equation, 2023, 27, 1027-5487, 10.11650/tjm/230801
    3. A. Oulmelk, L. Afraites, A. Hadri, M.A. Zaky, A.S. Hendy, Alternating direction multiplier method to estimate an unknown source term in the time-fractional diffusion equation, 2024, 156, 08981221, 195, 10.1016/j.camwa.2023.12.027
    4. Ivan Polansky, Kirill Loginov, Optimal Nonlinear Filtering of Information Impact Estimates in a Stochastic Model of Information Warfare, 2023, 22, 2713-3206, 745, 10.15622/ia.22.4.2
    5. Liangliang Sun, Zhaoqi Zhang, Yunxin Wang, The quasi-reversibility method for recovering a source in a fractional evolution equation, 2025, 1311-0454, 10.1007/s13540-025-00370-z
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1775) PDF downloads(105) Cited by(5)

Figures and Tables

Figures(6)  /  Tables(13)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog