# Inverse Source Problem for a Multiterm Time-Fractional Diffusion Equation with Nonhomogeneous Boundary Condition

Authors: L. L. Sun; X. B. Yan
Journal: Advances in Mathematical Physics
 (2020)

---

## Abstract

This paper is devoted to identify a space-dependent source function in a multiterm time-fractional diffusion equation with nonhomogeneous boundary condition from a part of noisy boundary data. The well-posedness of a weak solution for the corresponding direct problem is proved by the variational method. We firstly investigate the uniqueness of an inverse initial problem by the analytic continuation technique and the Laplace transformation. Then, the uniqueness of the inverse source problem is derived by employing the fractional Duhamel principle. The inverse problem is solved by the Levenberg-Marquardt regularization method, and an approximate source function is found. Numerical examples are provided to show the effectiveness of the proposed method in one- and two-dimensional cases.


---

## Body

## 1. Introduction

It is well known that the standard diffusion equation has been used to describe the Gaussian process of particle motion. Time-fractional diffusion equations (TFDEs in short) are deduced by replacing the standard time derivative with a time-fractional derivative and can be used to describe anomalous diffusion phenomena. Anomalous diffusion deviates from the standard Fickian description of Brownian motion, the main character of which is that its mean squared displacement is a nonlinear growth with respect to time, such asx2t∼tα. However, as the research problems become more and more complex, the differential order of a time-fractional diffusion equation is no longer limited to a fixed number but distributed over the unit interval. This creates TFDEs of distributed order of which a particular case is the multiterm time-fractional diffusion equations (MTFDEs in short).Direct problems, i.e., initial value problems and initial boundary value problems for MTFDEs, have attracted much more attentions in recent years. For example, the maximum principle in [1, 2] and the well-posedness and long-time asymptotic behavior in [3] for general multiterm time-fractional diffusion equations are investigated. The numerical solutions are shown in [4, 5] by the high-order space-time spectral method and the Galerkin finite element method.However, in most cases, the parameters that characterize the diffusion process cannot be measured directly or easily. The inverse source problem of diffusion process is aimed at detecting the source function of a physical field from some indirect measurements (such as final time information or boundary measurement) and is of great importance in engineering. There are many ripe regularization theories on the inverse source problems at present, e.g., the reproducing kernel Hilbert space method [6], the Fourier truncation method [7], the Tikhonov regularization method [8], and the modified quasiboundary value method and quasireversibility method [9, 10].In this paper, we investigate an inverse space-dependent source problem in a multiterm time-fractional diffusion equation with nonhomogeneous boundary condition. LetΩ be a bounded domain in ℝd with sufficient smooth boundary ∂Ω. For a fixed positive integer m, let qj and α=αj:j=1,2,⋯,m be the positive constants such that 0<αm<⋯<α1<1. Consider the following initial boundary value problem (IBVP) for a time-fractional diffusion equation
(1)∑j=1mqj∂0+αjux,t−Δux,t=fxrt,x,t∈ΩT≔Ω×I,ux,0=ϕx,x∈Ω¯,∂u∂nx,t=bx,t,x,t∈∂Ω×I,where I=0,T and n is the unitary outer normal vector of ∂Ω. We can assume q1=1 without loss of generality.Here,∂0+αj denotes the Caputo derivative defined by
(2)∂0+αjux,t=1Γ1−αj∫0t∂ux,s∂sdst−sαj,t>0,where Γ· is the Gamma function.If all functionsfx, rt, ϕx, and bx,t, and the parameters m, qj and αj are given appropriately, problem (1) is a direct problem. The inverse problem here is to determine the source term fx based on problem (1) and additional data
(3)ux,tΓ=hx,t,t∈I,where Γ is a nonempty open part of ∂Ω.To our best knowledge, it has made significant progress in determining the source problem for (1) with part of boundary data. For a single term (i.e. m=1), Wei et al. identified a time-dependent source term in a multidimensional TFDE from boundary Cauchy data. Zhang et al. [11] proved a uniqueness result to inverse the space-dependent source term in one-dimensional case by using one-point Cauchy data and provided an efficient numerical method. For a general area of high dimension, Wei et al. [12] proved the uniqueness of inverse space-dependent source but no numerical computation. Yan et al. [13] studied an inverse spatial-dependent source problem by noisy boundary data in a time-fractional diffusion-wave equation and carried out the numerical inversion experiments by a nonstationary iterative Tikhonov regularization method. For a multiterm case, Jiang et al. [14] established a weak unique continuation property for time-fractional diffusion-advection equations and studied an inverse problem on determining the spatial component in the source term by interior measurements. Very recently, Li et al. [15] investigated an inverse time-dependent source term problem in a MTFDE from the boundary Cauchy data and employed the conjugate gradient method to find the approximate source term.Nevertheless, as far as the authors know, the most of the existing literature only treated the inverse source problems with homogeneous boundary condition, in which the series expression provides convenience for the argument. The case of nonhomogeneous boundary is also significant from the practical point of view but difficult whether the forward problem or the inverse problem. Although the weak unique continuation property in [14] is also valid for the boundary measurements only taking a zero extension, however, it does not hold true for the nonhomogeneous boundary case. So far, there is no publication on inversion space-dependent source term by the Cauchy data in a MTFDE with nonhomogeneous boundary condition.Generally speaking, the problems of recovering spatial information from data along a time trace are notoriously severely ill posed. In this paper, we focus on an inverse space-dependent source problem of (1) from the Cauchy data in a general domain. This is an extension and improvement of [12], which dealt with an inverse space-dependent source problem by the boundary measured data in an infinite time interval for a single-term equation with homogeneous boundary. In this study, however, only part of the boundary Cauchy data in a finite time interval t∈0,TT<∞ is enough to obtain the uniqueness of inverse problem. We firstly obtain the well-posedness for the direct problem by the variational method for studying the inverse problem. For inverse problem, roughly speaking, we first transfer the original forward problems into three IBVPs by the superposition principle and turn the inverse source problem into an inverse initial problem by employing the fractional Duhamel principle. Then, we prove the uniqueness of the inverse initial problem by the analytic continuation technique and the Laplace transformation. Finally, we employ the Levenberg-Marquardt regularization method to solve numerically the inverse source problem. Here, we point out that the Levenberg-Marquardt method is modified to a more efficient algorithm without bringing in a differential step in the present paper because of the linearity of the operator. The numerical results for three examples in one- and two-dimensional cases are provided, and the numerical implementation shows the effectiveness and robustness of the proposed methods.The main result in this paper is the following uniqueness result for the inverse space-dependent source problem.Theorem 1.
Assume thatf∈D−∆+1γ for γ>max0,d/2−1 and r∈C10,T with r0≠0, the initial function ϕ∈H2Ω, and the boundary function b∈L2∂Ω×0,T. Suppose uix,t,i=1,2 is the solutions of (1) with f=fi, then u1x,t=u2x,t for x,t∈Γ×I implies f1x=f2x.The remainder of this paper is organized as follows. Some preliminaries are presented in Section2. The well-posedness for the direct problem is proved in Section 3. In Section 4, we present the uniqueness result of inverse space-dependent source problem. In Section 5, we use the Levenberg-Marquardt regularization method to find the approximate space-dependent source function. Numerical results for three examples are provided to illustrate the efficiency of our method in Section 6. Finally, we give a brief conclusion in Section 7.


## 2. Preliminaries

We firstly introduce some preliminaries as follows in this section.Definition 2 (see [16, 17]).
Ifft∈L10,T, then for α>0, the Riemann-Liouville fractional left-sided integral I0+αf and right-sided integral IT−αf are defined by
(4)I0+αft=1Γα∫0tfsdst−s1−α,0<t≤T,IT−αft=1Γα∫tTfsdss−t1−α,0≤t<T.Definition 3. (see [16, 17]).
Letzt∈AC0,T, then for 0<α<1, the Riemann-Liouville fractional left-sided derivative D0+αyt and right-sided derivative DT−αyt of order α are defined by
(5)D0+αyt=1Γ1−αddt∫0tysdst−sα≕ddtI0+1−αyt,0<t≤T,DT−αyt=−1Γ1−αddt∫tTysdss−tα≕−ddtIT−1−αyt,0≤t<T.Definition 4 (see [16, 18]).
The multinomial Mittag-Leffler function is defined by(6)Eθ1,⋯,θm,θ0z1,⋯,zm≔∑k=0∞∑k1+⋯+km=kk;k1,⋯,kmΠj=1mzjkjΓθ0+∑j=1mθjkj,where θ0,θj∈ℝ, and zj∈ℂj=1,⋯,m, and k;k1,⋯,km denotes the multinomial coefficient
(7)k;k1,⋯,km≔k!k1!⋯km!withk=∑j=1mkj.
For later use, we adopt the abbreviation(8)Eα′,βλnt≔Eα1,α1−α2,⋯,α1−αm,β−λntα1,−q2tα1−α2,⋯,−qmtα1−αm,t>0.
We have the following three results on the multinomial Mittag-Leffler function [3].Lemma 5.
Letθ0>0, 0<θj<1j=1,⋯,m and zj∈ℂj=1,⋯,m be fixed.
Then,(9)1Γθ0+∑j=1mzjEθ1,⋯,θm,θ0+θjz1,⋯,zm=Eθ1,⋯,θm,θ0z1,⋯,zm.Lemma 6.
Letθ0>0,1>α1>⋯>αm>0 be given. Assume that α1π/2<μ<α1π, μ≤∣argz1∣≤π and π−ϵ≤∣argzj∣≤π, j=2,⋯,m with ϵ small enough. Then, there exists a constant C>0 depending only μ, K, αj(j =1, · · ·,m) and β such that
(10)Eα1,α1−α2,⋯,α1−αm,βz1,⋯,zm≤C1+z1.Lemma 7.
Let1>α1>⋯>αm>0. Then,
(11)ddttα1Eα′,1+α1λnt=tα1−1Eα′,α1λnt,t>0.Lemma 8.
Assumer∈L∞0,T denote
(12)rnt=∫0trst−sα1−1Eα′,α1λnt−sds,t∈0,T,n=1,2,⋯,then rnt∈C0,T and rn0=0.Proof.
Letent=tα1−1Eα′,α1λnt. Then, rnt=r∗ent. By Lemma 6, we have
(13)enL10,T=∫0Ttα1−1Eα′,α1λntdt≤C∫0Ttα1−1dt<∞.
Therefore, we arrive aten∈L10,T.
Using the Young theorem, we obtainrnt∈C0,T. Combining Lemma 7, we have
(14)rn0=limt→0rnt≤limt→0rL∞0,ttα1Eα′,1+α1λnt,∀t∈0,T.
Thus, we havern0=0.Lemma 9.
Denoteant=1−λntα1Eα′,1+α1λnt, t∈0,+∞, n=1,2,⋯, then
(15)Lans=∑j=1mqjsαj−1∑j=1mqjsαj+λn,s∈ℂwithRes>0,where Lan denotes the Laplace transform of function ant.Proof.
First, we can take the Laplace transforms termwise for∫0∞e−ztλntα1Eα′,1+α1λntdt with Rez>λn1/α1 by its definition (8). Since λntα1Eα′,1+α1λnt≤C by Lemma 6, we see that ∫0∞e−ztλntα1Eα′,1+α1λnt is analytic with respect to z in Rez>0. Therefore, the analytic continuation can be done for Rez>0. By Theorem 4.1 in [18], an(t) satisfies an initial value problem for a fractional ODE
(16)∑j=1mqj∂0+αjant+λnant=0,t>0,an0=1.
Taking the Laplace transform to (16) and using the formula
(17)L∂0+αfs=sαLfs=sα−1f0+,we obtain
(18)Lans=∑j=1mqjsαj−1∑j=1mqjsαj+λn.


## 3. Existence and Uniqueness of a Weak Solution for the Direct Problem

In this section, we will obtain the existence and uniqueness of a weak solution for direct problem (1). First, we define a Hilbert scale space D−∆+1γ in L2Ω. Denote the eigenvalues of −∆ with homogeneous Neumann boundary condition as λn and the corresponding eigenfunctions as φn∈ψ∈H2Ω;∂ψ/∂n∂Ω=0 such that ∆φn=−λnφn.We know that0=λ1<λ2≤λ3≤⋯≤λn≤⋯ counting according to the multiplicities, and φn is an orthonormal basis in L2Ω. Define the Hilbert scale space
(19)D−Δ+1γ=ψ∈L2Ω;∑n=1∞λn+12γψ,φn2<∞,where ·,· is the inner product in L2Ω and define its norm
(20)ψD−Δ+1γ=∑n=1∞λn+12γψ,φn21/2.Here, one points out that the number 1 in the operator−∆+1 is not necessary. It may be any nonzero positive number.In the following, we assume that the initial dataϕx∈H2Ω, the boundary data bx,tϵL2∂Ω×I, and the source functions,fxϵL2Ωrt∈L20,T. We give a weak formulation for direct problem (1) and prove that its solution exists uniquely.Letux,t be the smooth solution of (1). Denote vx,t=ux,t−ϕx, then v solves the following IBVP
(21)∑j=1mqj∂0+αjvx,t−Δvx,t=f¯x,t,x,t∈ΩT,vx,0=0,x∈Ω¯,∂v∂nx,t=b−∂ϕ∂n,x,t∈∂Ω×I,where f¯x,t=fxrt+∆ϕx.LetHα10,T be the fractional Sobolev space on the interval (0,T) (see, e.g., Adams [19]).We set(22)0Ha10,T=u∈Ha10,T:u0=0,if 1/2<α1≤1 and we identify 0Ha10,T with Hα10,T for 0<α1<1/2. For 0<α1<1, we define a space
(23)Bα1ΩT≔H0α10,T;L2Ω∩L20,T;H1Ω,equipped with the norm
(24)uBα1ΩT=uHα10,T;L2Ω2+uL20,T;H1Ω21/2.Referring to [20], we know that Bα1ΩT is a Hilbert space, and for α1∈0,1/2, we have uHα0,T;L2Ω2≅uL20,T;L2Ω2+D0+αuL20,T;L2Ω2≅uL20,T;L2Ω2+DT−αuL20,T;L2Ω2, where D0+αu and DT−αu denote the Riemann-Liouville fractional left derivative and right derivative of order α (see Definition 3). Based on [20–22], we can deduce a weak formulation for problem (21) as follows.Findν∈Bα1/2ΩT such that
(25)Av,w=Fw,∀w∈Bα1/2ΩT,where
(26)Av,w≔∑j=1mqjD0+αj/2v,DT−αj/2wΩT+∇v,∇wΩT,Fw≔f¯,wΩT+b−∂ϕ∂n,w∂Ω×I.Now, we could obtain the following theorem.Theorem 10.
Letf∈L2Ω, rt∈L20,T, and ϕ∈H2Ω. Then, there exists a unique solution vx,tto (25) and the solution satisfies
(27)vBα1/2ΩT≤Cα1,T,qjfL2ΩrL20,T+ϕH2Ω+bL2∂Ω×I,where Cα1,T,qj>0 is a constant independent of v.Proof.
The well-posed of problem (25) is guaranteed by the well-known Lax-Milgram theorem. The continuity of the functional F is obvious. Here, we give the continuity and coercivity of the bilinear form A. By the Cauchy-Schwarz inequality, we have
(28)Av,w≤∑j=1mqjD0+αj/2vL2ΩTDT−αj/2wL2ΩT+∇vL2ΩT∇wL2ΩT≤Cmmax1≤j≤mqjvBα1/2ΩTwBα1/2ΩT+vBα1/2ΩTwBα1/2ΩT≤1+Cmmax1≤j≤mqjvBα1/2ΩTwBα1/2ΩT.
On the other hand, we arrive at(29)Aw,w≥∑j=1mqjcosπαj/2D0+α1/2wL2ΩT2+∇wL2ΩT2≥cosπα1/2D0+α1/2wL2ΩT22+cosπα1/2wL2ΩT22+∇wL2ΩT2≥ccosπα1/22wBα1/2ΩT2,where we use Lemma 2.2 of [20] in the first inequality and the embedding theorem in the second inequality. In order to obtain the stability of (27), we can take w=v in (25) and use (29) to obtain
(30)vBα1/2ΩT≤Ccosπα1/2fL2ΩrL20,T+ϕH2Ω+bL2∂Ω×I.


## 4. Uniqueness for the Inverse Problem

In order to prove our main result, we firstly introduce the following famous fractional Duhamel principle.Assumeu solves equation (1) with ϕ=0 and b=0, and let v be the solution of the following problem
(31)∑j=1mqj∂0+αjvx,t−Δvx,t=0,x,t∈ΩT,vx,0=fx,x∈Ω¯,∂v∂nx,t=0,x,t∈∂Ω×I.Lemma 11. (see Lemma 4.2 of [2]).
Letu be the solution to (1) with ϕ=0 and b=0, where r∈C10,T and f∈D−∆+1ϵ with some ϵ>0. Then, u allows the representation
(32)u⋅,t=∫0tμt−svf⋅,sds,0<t≤T,where vf solves homogeneous problem (31) with f as the initial data, and μ satisfies
(33)∑j=1mqjI0+1−αjμt=rt,0<t≤T.
Moreover, there exists a uniqueμ∈L10,T satisfying (33), and there is a constant C>0 independent of t such that
(34)∣μt∣≤Ctα1−1a.e.t∈0,T.
Based on the method of separation of variables, we could obtain a formal solution to the homogeneous boundary problem of (1) given by
(35)ux,t=∑n=1∞antϕnφnx+∑n=1∞rntfnφnx,where ant=1−λntα1Eα′,1+α1λnt,rnt=∫0trst−sα1−1Eα′,1+α1λnt−sds,ϕn=ϕ,φn and fn=f,φn.
Now, we proceed to the proof of the uniqueness of the inverse initial problem for the homogeneous initial-boundary value problem (31).Theorem 12.
Assume thatf and f~∈D−∆+1γ for γ>max0,d/2−1. Suppose vx,t and v~x,t are the weak solutions of (31) with respect to initial conditions f and f~, respectively, then vx,t=v~x,t for x,t∈Γ×I implies fx=f~xinL2Ω.Proof.
By (35), the expressions vx,t and v~x,t are given by
(36)vx,t=∑n=1∞antfnφnx,v~x,t=∑n=1∞antf~nφnx,where f~n=f~,φn. Next, we prove the uniform convergence of above series. By the Sobolev embedding theorem, we have
(37)φnCΩ¯≤CφnH2mΩ,for m>d/4. Moreover, we have
(38)φnH2mΩ≤C−∆+1mφnL2Ω=Cλn+1m.
On the one hand, whenz∈S1≔z∈ℂ;argz≤minπ−μ/α1,ϵ/α1−αm,π, by Lemma 5 and Lemma 6, combining (37) and (38) and the Cauchy-Schwarz inequality, we arrive at
(39)∑n=2∞maxx∈Ω¯1−λnzα1Eα′,1+α1λnzfnφnx=∑n=2∞maxx∈Ω¯Eα′,1λnz+∑j=2mqjzα1−αjEα′,1+α1−α2λnzfnφnx≤C∑n=2∞∑j=1mzα1−αj1+λnzα1fnφnxL∞Ω¯≤C∑n=2∞fnλn+1m−1≤C∑n=2∞1λn+12γ+1−m1/2∑n=2∞λn+12γfn1/2.
By the book [23], we know that λn≥Cn2/d,n∈ℕ. If we choose γ>m+d/4−1, i.e.,γ>d/2−1, then the first series of the last term is convergent. By f∈D−∆+1γ, we know that the above series is convergent uniformly over z∈S1 for x∈Ω¯ . On the other hand, we know that the multinomial Mittag-Leffler function Eα′,1+α1λnz is analytic over the domain S2≔z∈ℂ;z>0,argz<π based on its definition. Thus, the expressions vx,tand v~x,t are analytic over the domain z∈S1∩S2 for x∈Ω¯ based on the Weierstrass theorem. Especially, ux,t and vx,t can be extended to Ω¯×0,∞.
Therefore, by the analytic continuation, the conditionvx,t=v~x,t for x,t∈Γ×I gives
(40)∑n=1∞antfnφnx=∑n=1∞antf~nφnx,x∈Γ,t>0.
We setσ−∆=μnn∈ℕ with 0=μ1<μ2<μ3<⋯ and denote by φnk1≤k≤mn an orthonormal basis of Kerμn+∆. Now, we consider σ−∆ as a set, not as a sequence with multiplicities. Then, we can rewrite (40) by
(41)∑n=1∞∑k=1mnfnkφnkxa~nt=∑n=1∞∑k=1mnf~nkφnkxa~nt,x∈Γ,t>0,where a~nt=1−μntα1Eα′,1+α1μnt and fnk=f,φnk. From (39), one obtains
(42)e−tRez∑n=2∞fnantφnx≤Ce−tRezfD−Δ+1γ,and e−tRez is integrable in t∈0,+∞ for fixed z satisfying Rez>0. By the Lebesgue convergence theorem, we can take the Laplace transform for (41). By Lemma 9, we have
(43)Lans=∑j=1mqjsαj−1∑j=1mqjsαj+λn,Res>0.
Therefore, we arrive at(44)∑n=1∞∑k=1mnfnkφnkx∑j=1mqjsαj−1∑j=1mqjsαj+μn=∑n=1∞∑k=1mnf~nkφnkx∑j=1mqjsαj+1∑j=1mqjsαj+μn,x∈Γ,Res>0,which implies
(45)∑n=1∞∑k=1mnfnkφnkxηη+μη=∑n=1∞∑k=1mnf~nkφnkxηη+μn,x∈Γ,η∈S3⊃ℝ+,where η=∑j=1mqjsαj. As
(46)∑n=1∞∑k=1mnf,φnkηη+μnφnkx≤∑n=1∞f,φnφnL∞Ω¯ηη+λn,one can see that the above series is internally closed uniform convergence in η∈ℂ\−μnn=1∞ from (39). Using the Weierstrass theorem, both sides of (45) are analytic in η∈ℂ\−μnn=1∞. Therefore, one can analytically continue η such that (45) holds for η∈ℂ\−μnn=1∞.
Now, we can take a suitable disk which only includesμl. Using the Cauchy integral theorem, integrating (45) along this disk, we obtain that μlx=∑k=1mlf−f~φlkφlkx=0,x∈Γ. Since Δ+μlul=0 in Ω, and ul=∂ul/∂n=0, on Γ the uniqueness of the Cauchy problem for elliptic equations (e.g., see Theorem 3.3.1 of Isakov [24], p. 58) implies ul=0 in Ω for each l∈ℕ. Combining the linearly independence of φlk1≤k≤ml in Ω, we obtain that f−f~,φlk=0 for 1≤k≤ml,l∈ℕ. Therefore f=f~ in L2Ω.Remark 13.
Assume thatϕ=0 and b=0, f∈D−Δ+1γ for γ>max0,d/2−1 and r∈L∞0,∞, then the solution of (1) u∈CΩ¯×0,+∞.Proof.
From (35), the solution ux,t is given by
(47)ux,t=∑n=1∞rntfnφnx.
Lett0>0 be arbitrarily fixed. For t≥t0, by Lemma 7 and Lemma 6, we have
(48)∑n=2∞maxx∈Ω¯∫0trst−sa1−1Eα′,a1λnt−sdsfnφnx≤r∞∑n=2∞fnφnL∞Ω¯Ctα11+λntα1≤Cr∞∑n=2∞fnλn+1m−1≤Cr∞∑n=2∞1λn+12γ+1−m1/2∑n=2∞λn+12γfn21/2≤Cr∞fD−Δ+1γ,where γ>max0,d/2−1. Thus, the series (48) is also convergent onΩ¯×t0,∞ uniformly. By Lemma 8, we know u in (48) is continuous on Ω¯×t0,∞.Finally, we can proceed to prove the main result Theorem1 raised at the end of the introduction.Proof of Theorem 1.
As we know, by the principle of linear superposition, we haveu=uf+uϕ+ub and they solve the following three IBVPs:
(49)∑j=1mqj∂0+αjufx,t−Δfx,t=fxrt,x,t∈ΩT,ufx,0=0,x∈Ω¯,∂uf∂nnx,t=0,x,t∈∂Ω×I,(50)∑j=1mqj∂0+αjuϕx,t−Δuϕx,t=0x,t∈ΩT,uϕx,0=ϕx,x∈Ω¯,∂uϕ∂nx,t=0,x,t∈∂Ω×I,(51)∑j=1mqj∂0+αjubx,t−Δubx,t=0,x,t∈ΩT,ubx,0=0,x∈Ω¯,∂ub∂nx,t=bx,t,x,t∈∂Ω×I.
It is not hard to find that the solutions of IBVPs (50) and (51) are independent of the source term f, and then, the inverse source problem becomes more precise, i.e., can we determine fx from ufx,t for x,t∈Γ×I uniquely. According to Remark 13, we know that uf∈CΩ¯×0,∞. By Lemma 11, we know that ufx,t=0, x∈Γ, 0<t≤T deduce
(52)∫0tμt−svfx,sds=0,x∈Γ,0<t≤T.
By the Fibini theorem, we obtain from (33) that
(53)∑j=1mqjI0+1−αj∫0tμt−svfx,sds=∑j=1mqjΓ1−αj∫0tdτt−ταj∫0τμτ−svfx,sds,=∑j=1mqjΓ1−αj∫0tvfx,sds∫0tμτ−st−ταjdτ,=∫0trt−svfx,sds.
So, we have∫0trt−svfx,sds=0, x∈Γ, 0<t≤T. Differentiating the above equality with respect to t, we arrive at r0vfx,t+∫0tr′t−svfx,sds=0. As the assumption r0≠0, we have by the generalized Minkowski inequality that
(54)vf⋅,tL2Γ≤1r0∫0tr′t−svf⋅,sL2Γds≤rC10,Tr0∫0tvf⋅,sL2Γds.
Using the Gronwall inequality, we havevfx,t=0, x∈Γ, 0<t≤T. Therefore, we know fx=0 by Theorem 12. Thus, we complete the proof.


## 5. Levenberg-Marquardt Regularization Method

In this section, we solve numerically the inverse problem of identifying the space-dependent sourcefx. As we know, most of inversion algorithms are based on regularization strategies so as to overcome ill-posedness of inverse problems, and different kinds of inverse problems may need different approximate methods on the basis of conditional well-posedness analysis. In this paper, we employ the Levenberg-Marquardt regularization technique to obtain an approximation to the source term.Based on Theorem10, we can define a forward linear operator
(55)F:fx∈L2Ω→ux,t;f∈L2Γ×0,T.Thus, the inverse problem is translated into solving the following abstract operator equation(56)Ff=ux,tΓ≜hx,t,t∈I.From Theorem10, we know that ux,t;f∈Bα1/2ΩT. By the trace theorem in [25], we know ux,t;fΓ×I∈L20,T;H1/2Γ∩0Ha1/40,T;L2Γ↬L2Γ×0,T compactly, then the operator F:L2Ω→L2Γ×0,T is compact and thus the inverse source problem is ill-posed. In order to ensure a stable numerical reconstruction of fx, we give the following minimization problem with the Tikhonov-type regularization term
(57)minf∈D−Δ+1γJμf=Ff−hδL2Γ×0,T2+μf−f∗D−Δ+1γ2,where γ>d/2−1, μ>0 is a regularization parameter, f∗∈D−Δ+1γ is a suitable guess of f, and hδ is the noisy function of h.Lemma 14. (see [26]).
Assume thatE is a uniformly convex Banach space. Let (xn) be a sequence in E such that xn⇀x weakly σE,E′ andlimsup∥xn∥≤∥x∥. Then, xn⇀x strongly.
By the canonical book [27], we know that with above minimization problem exists a unique stable solution of regularization. Moreover, the convergence result can be obtained in the following theorem.Theorem 15.
Suppose thatFf=h and the noisy data hδ∈L2Γ×0,T satisfying hδ−hL2Γ×0,T≤δ and let μδ satisfy μδ→0 and δ2/μδ→0 as δ→0; then, the minimizer fμδ of variational problem (57) is convergent, i.e., fμδ→f in D−Δ+1γ as δ→0.Proof.
Letδk be any sequence such that δk→0 and denote μk=μδk. Since
(58)Ffμkδk−hδkL2Γ×0,T2+μkfμkδk−f∗D−Δ+1γ2≤Ff−hδkL2Γ×0,T2+μkf−f∗D−Δ+1γ2≤δk2−μkf−f∗D−Δ+1γ2,we have fμkδk−f∗D−Δ+1γ2≤δk2/μk+f−f∗D−Δ+1γ2 and limsupδk→0fμkδk−f∗D−Δ+1γ≤f−f∗D−Δ+1γ by the condition for μ. So we know that fμkδkDΔ+1γ is bounded and has a weak convergence subsequence by the reflexivity of D−Δ+1γ and also denotes fμkδk such that fμkδk⇀zinD−Δ+1γ.
Since the linear operatorF is bounded, we know that Ffμkδk⇀FzinL2Γ×0,T.
From (58), we have Ffμkδk−hδkL2Γ×0,T≤δk2+μkf−f∗D−Δ+1γ2, and thus limδk→0fμkδk=h and Fz=h. By the uniqueness of Ff=h, we know that z=f.
By the weak lower semicontinuity of the norm in the Hilbert space, we havef−f∗D−Δ+1γliminfδk→0fμkδk−f∗D−Δ+1γ≤limsupδk→0fμkδk−f∗D−Δ+1γ≤f−f∗D−Δ+1γ, and hence, we know that limδk→0fμkδk−f∗D−Δ+1γ=f−f∗D−Δ+1γ. Combing with the weak convergence of fμkδk, we have by Lemma 14 that fμkδk→f. From the uniqueness of Ff=h, we have that fμδδ→f in D−Δ+1γ.In the following, we use the Levenberg-Marquardt method to minimize problem (57). The Levenberg-Marquardt method was first introduced by [28, 29], and it is a kind of the Newton-type method. Moreover, it has been well applied to fractional equations, for example, see [30–32]. From physical considerations, we know that hδ is a reasonably close approximation of some ideal h=Ff in the range of F. Let f∗ be an approximation of f. By the linearity expansion for linear operator F at f∗, we can obtain
(59)Ff=Ff∗+Ff′f∗f−f∗.Then, the inverse problemFf=hδ can be transformed into
(60)Ff′f∗f−f∗=hδ−Ff∗.Therefore, it is easily seen that the problem of minimizing (57) is equivalent to minimizing
(61)minδf∈D−Δ+1γJμf=Ff′f∗δf−hδ−Ff∗L2Γ×0,T2+μδfD−Δ+1γ2,where δf=f−f∗.Now, we consider the discretization of the minimization problem. We define an admissible set of unknown sourcefx. Suppose that φsx,s=1,2,⋯,∞ is a basis in D−Δ+1γ, let
(62)fx≈f^x=∑s=1sasφsx,f∗x≈f∗^x=∑s=1sas∗φsx,where f^ is the S-dimensional approximation solution to f and S∈ℕ is a truncated level of f, and as,s=1,2,⋯, S are the expansion coefficients. It is convenient to set a finite dimensional space as
(63)ΦS=spanφ1,φ2,⋯,φS,and S-dimensional vector a=a1,a2,⋯,aS∈ℝS. We identify an approximation f^x∈ΦS with a vector a∈ℝS.Based on the above discussions, by settingux,t;a=ux,t;f^ as a unique solution of the forward problem, a feasible way for numerical solution to solve the following minimization problem
(64)minδa∈ℝSδa∇aTux,t;a∗−hδ−ux,t;a∗L2Γ×0,T2+μδaAδaT,where A=λiγφi,λjγφjL2ΩS×S, a∗=a1∗,a2∗,⋯,aS∗, δa=a−a∗, and aT denotes the transpose of a.Next, we give an inversion algorithm for determining coefficienta. Given the initial a0 and let a∗=a0. Suppose ak∈ℝS is the kth step iteration value, then we obtain (k+1)th step approximation by solving
(65)ak+1=ak+δak,k=0,1,⋯,where δak denotes a perturbation of ak for each k and k is the number of iterations. For convenience of writing, ak and δak are also abbreviated as a and δa, respectively.Now, we computeδa by the variational theory. As we know,
(66)∇aux,t;ak=ux,t;ak+τe1−ux,t;akτ,⋯,ux,t;ak+τeS−ux,t;akτ,whereei=0,⋯,0,1,0,⋯,0, i=1, ⋯, S. Since operator F is linear respect to the source f, we have
(67)∇aux,t;ak=ux,t;e1,⋯,ux,t;eS≕b1,⋯,bS.By the variational theory, problem (64) is equivalent to the following normal equations
(68)μA+BδaT=V,where B=b1,bjL2Γ×0,TS×S and V=hδx,t−ux,t;ak,biL2Γ×0,TS×1.


## 6. Numerical Experiments

In this section, we present the numerical results for three examples in one-dimensional and two-dimensional cases to show the effectiveness of the Levenberg-Marquardt method.The noisy data is generated by adding a random perturbation, i.e.,(69)hδ=h+∈h⋅2⋅randsizeh−1.The corresponding noise level is calculated byδ=hδ−hL20,T;L2Γ×0,zt.To show the accuracy of numerical solution, we compute the approximateL2 error denoted by
(70)rek=fkx−fxL2ΩfkxL2Ω,where fkx is the coefficient term reconstructed at the kth iteration, and fx is the exact solution.The residualEk at the kth iteration is given by
(71)Ek=ux,t;fk−hδx,tL2Γ×0,T.In an iteration algorithm, the most important work is to find a suitable stopping rule. In this study, we use the well-known discrepancy principle [33], i.e., we choose k satisfying the following inequality
(72)Ek≤ζδ<Ek−1,where ζ>1 is a constant and can be taken heuristically to be 1.01, as suggested by Hanke and Hansen [34].


### 6.1. One-Dimensional Case

Without loss of generality, we setm=3, qj=1j=1,2,3, and ΩT=0,1×0,T. The grid points on [0,1] and 0,T are both 201 when solving the direct problem by finite difference method in [35]. In this case, we take Γ as a point x=0 and ΦS is chosen as a subspace of eigenfunctions
(73)ΦS=span1,2cosπx,⋯,2cosS−1πx.Example 16.
We firstly test a smooth solution. Take an unknown sourcefx=5cos3πx+10x21−x2, the temporal source rt=expt, the initial ϕx=cosπx, and the boundary b0,t=1+t, b1,t=−1+t. The boundary data u0,t is obtained by solving the direct problem (1) by using the finite difference method. The numerical results for Example 16 by using the discrepancy principle for various noise levels in the cases of α=0.9,0.3,0,1 and α=0.7,0.6,0.2 are shown in Figures 1(a) and 1(b), respectively. We choose the initial guess as f0=f∗=0, the truncated level S=4.Example 17.
In the second example, we test a nonsmooth solution with a cusp. Take an unknown source(74)fx=1,0≤x<0.1,1+2.5x−0.1,0.5≤x<0.5,2+2.50.5−x,0.5≤x<0.9,1,0.9≤x≤1,and the rest of definite conditions are the same as Example 16. The boundary data u0,t is obtained by solving the direct problem (1) by using the finite difference method. The numerical results for Example 17 by using the discrepancy principle for various noise levels in the cases of α=0.9,0.6,0,1 and α=0.8,0.7,0.6 are shown in Figures 2(a) and 2(b), respectively. We choose the initial guess as f0=0 and f∗=0, the truncated level S=5.Figure 1
The numerical results for Example16 for various noise levels with μ=10−6δ3/4 and γ=0.
(a)
α=0.9,0.3,0.1(b)
α=0.7,0.6,0.2Figure 2
The numerical results for Example17 for various noise levels with μ=10−5δ3/4 and γ=1/10.
(a)
α=0.9,0.6,0.1(b)
α=0.8,0.7,0.6From Tables1 and 2, we can come to the following conclusions. First, it is observed that the dominate fractional order α1 has a greater influence on the accuracy of the numerical inversion results. For each fixed ϵ, we can see that the numerical results become worse as the dominate fractional order α1 becomes smaller, but it is not sensitive to the rest of the orders α2 or α3. On the other hand, the numerical results become worse as the index of penalty term γ becomes larger for each fixed ϵ. Since here we test a smooth function, the penalty term with high index may not have positive effect to the numerical inversion. So we choose γ=0 in the first example.Table 1
The relative errorrek with different noise levels ϵ and fractional orders α for Example 16.
α00.0010.0050.01rek∈0.9,0.3,0.20.00590.01350.01260.04840.8,0.7,0.60.00600.00820.03570.09410.7,0.5,0.20.00590.09000.22580.31040.6,0.3,0.20.00590.12170.60460.75690.5,0.3,0.20.00580.29470.76560.8595Table 2
The relative errorrek with different noise levels ϵ and index of penalty term γ for α=0.8,0.6,0.1 in Example 16.
γ00.0010.0050.01rek∈00.00590.02570.03480.02971/100.00590.01010.18510.23441/30.00590.03610.24620.48121/20.00590.09390.43630.5727From Table3, we can see that the regularity of the penalty term can do modify the computation result well for a nonsmooth numerical test. But this regularity should not be too strong, or it will have the opposite effect. We find that γ=1/10 is a better index by the numerical test, so we choose γ=1/10 in the second example.Table 3
The relative errorrek with different noise levels ϵ and index of penalty term γ for α=0.9,0.5,0.2 in Example 17.
α00.0010.0050.01rek∈−1/100.01940.06950.13720.188900.01940.02740.05320.0801−1/100.01940.02310.02610.02721/40.01940.02380.06870.07041/20.01940.03330.07600.0772


### 6.2. Two-Dimensional Case

In a two-dimensional case, we present one numerical example to show the effectiveness of our method. We also assumeT=1, m=3, qj=1j=1,2,3 without loss of generality. The space domain Ω is taken as 0,1×0,1. The grid points on [0,1] and 0,T are both 51 when solving the direct problem by finite difference method in [35]. We take Γ=x,yx=0,0≤y≤1. In this case, we take ΦS1×S2=spanφixφjy, where
(75)φix=1,i=1,2cosi−1πx,i>1,i=1,2,⋯,S1,φjy=1,j=1,2cosi−1πy,j>1,j=1,2,⋯,S2.Example 18.
Assumefx,y=cos2πxcos2πy+1/10exp−1/3x3−3x22y3−3y2 and rt=expt, ϕx,y=0, b=0. The boundary data u0,y,t is obtained by solving direct problem (1). The exact source solutions, reconstructed source terms, and the absolute errors between the exact source functions and the numerical solutions for α=0.9,0.5,0,1 and α=0.6,0.5,0,4 are shown in Figures 3 and 4 by taking noise level ϵ=0.01 and the regularization parameter μ=0.01δ3/4, respectively. Here, we choose the initial guess as f0=0 and f∗=0. The truncated levels are S1=S2=3.Figure 3
The numerical results for Example18 for α=0.9,0.5,0.1 with μ=0.01δ3/4.
(a)
Exact solution(b)
Regularized solution(c)
Absolute errorFigure 4
The numerical results for Example18 for α=0.6,0.5,0.4 with μ=0.01δ3/4.
(a)
Exact solution(b)
Regularized solution(c)
Absolute errorFrom Figures3 and 4, we can see the proposed method still yields relatively accurate numerical solutions for a two-dimensional example.Remark 19.
As we know, we need to compute the gradient of an abstract operator by the Levenberg-Marquardt method, and we usually use numerical differential method to approximately substitute the gradient when it is hard to obtain the exact expression for the gradient. So, we need to introduce a numerical differentiation stepτ. However, a lot of numerical examples are sensitive to τ. In this computation, we have successfully avoided the problem of selection of the differential step. This is because that the operator F is linear to f and the gradient in (66) can be replaced by (67) where the gradient is not dependent to τ. But when we calculate (66) by a finite difference, the initial and boundary conditions need to be homogeneous.


## 6.1. One-Dimensional Case

Without loss of generality, we setm=3, qj=1j=1,2,3, and ΩT=0,1×0,T. The grid points on [0,1] and 0,T are both 201 when solving the direct problem by finite difference method in [35]. In this case, we take Γ as a point x=0 and ΦS is chosen as a subspace of eigenfunctions
(73)ΦS=span1,2cosπx,⋯,2cosS−1πx.Example 16.
We firstly test a smooth solution. Take an unknown sourcefx=5cos3πx+10x21−x2, the temporal source rt=expt, the initial ϕx=cosπx, and the boundary b0,t=1+t, b1,t=−1+t. The boundary data u0,t is obtained by solving the direct problem (1) by using the finite difference method. The numerical results for Example 16 by using the discrepancy principle for various noise levels in the cases of α=0.9,0.3,0,1 and α=0.7,0.6,0.2 are shown in Figures 1(a) and 1(b), respectively. We choose the initial guess as f0=f∗=0, the truncated level S=4.Example 17.
In the second example, we test a nonsmooth solution with a cusp. Take an unknown source(74)fx=1,0≤x<0.1,1+2.5x−0.1,0.5≤x<0.5,2+2.50.5−x,0.5≤x<0.9,1,0.9≤x≤1,and the rest of definite conditions are the same as Example 16. The boundary data u0,t is obtained by solving the direct problem (1) by using the finite difference method. The numerical results for Example 17 by using the discrepancy principle for various noise levels in the cases of α=0.9,0.6,0,1 and α=0.8,0.7,0.6 are shown in Figures 2(a) and 2(b), respectively. We choose the initial guess as f0=0 and f∗=0, the truncated level S=5.Figure 1
The numerical results for Example16 for various noise levels with μ=10−6δ3/4 and γ=0.
(a)
α=0.9,0.3,0.1(b)
α=0.7,0.6,0.2Figure 2
The numerical results for Example17 for various noise levels with μ=10−5δ3/4 and γ=1/10.
(a)
α=0.9,0.6,0.1(b)
α=0.8,0.7,0.6From Tables1 and 2, we can come to the following conclusions. First, it is observed that the dominate fractional order α1 has a greater influence on the accuracy of the numerical inversion results. For each fixed ϵ, we can see that the numerical results become worse as the dominate fractional order α1 becomes smaller, but it is not sensitive to the rest of the orders α2 or α3. On the other hand, the numerical results become worse as the index of penalty term γ becomes larger for each fixed ϵ. Since here we test a smooth function, the penalty term with high index may not have positive effect to the numerical inversion. So we choose γ=0 in the first example.Table 1
The relative errorrek with different noise levels ϵ and fractional orders α for Example 16.
α00.0010.0050.01rek∈0.9,0.3,0.20.00590.01350.01260.04840.8,0.7,0.60.00600.00820.03570.09410.7,0.5,0.20.00590.09000.22580.31040.6,0.3,0.20.00590.12170.60460.75690.5,0.3,0.20.00580.29470.76560.8595Table 2
The relative errorrek with different noise levels ϵ and index of penalty term γ for α=0.8,0.6,0.1 in Example 16.
γ00.0010.0050.01rek∈00.00590.02570.03480.02971/100.00590.01010.18510.23441/30.00590.03610.24620.48121/20.00590.09390.43630.5727From Table3, we can see that the regularity of the penalty term can do modify the computation result well for a nonsmooth numerical test. But this regularity should not be too strong, or it will have the opposite effect. We find that γ=1/10 is a better index by the numerical test, so we choose γ=1/10 in the second example.Table 3
The relative errorrek with different noise levels ϵ and index of penalty term γ for α=0.9,0.5,0.2 in Example 17.
α00.0010.0050.01rek∈−1/100.01940.06950.13720.188900.01940.02740.05320.0801−1/100.01940.02310.02610.02721/40.01940.02380.06870.07041/20.01940.03330.07600.0772


## 6.2. Two-Dimensional Case

In a two-dimensional case, we present one numerical example to show the effectiveness of our method. We also assumeT=1, m=3, qj=1j=1,2,3 without loss of generality. The space domain Ω is taken as 0,1×0,1. The grid points on [0,1] and 0,T are both 51 when solving the direct problem by finite difference method in [35]. We take Γ=x,yx=0,0≤y≤1. In this case, we take ΦS1×S2=spanφixφjy, where
(75)φix=1,i=1,2cosi−1πx,i>1,i=1,2,⋯,S1,φjy=1,j=1,2cosi−1πy,j>1,j=1,2,⋯,S2.Example 18.
Assumefx,y=cos2πxcos2πy+1/10exp−1/3x3−3x22y3−3y2 and rt=expt, ϕx,y=0, b=0. The boundary data u0,y,t is obtained by solving direct problem (1). The exact source solutions, reconstructed source terms, and the absolute errors between the exact source functions and the numerical solutions for α=0.9,0.5,0,1 and α=0.6,0.5,0,4 are shown in Figures 3 and 4 by taking noise level ϵ=0.01 and the regularization parameter μ=0.01δ3/4, respectively. Here, we choose the initial guess as f0=0 and f∗=0. The truncated levels are S1=S2=3.Figure 3
The numerical results for Example18 for α=0.9,0.5,0.1 with μ=0.01δ3/4.
(a)
Exact solution(b)
Regularized solution(c)
Absolute errorFigure 4
The numerical results for Example18 for α=0.6,0.5,0.4 with μ=0.01δ3/4.
(a)
Exact solution(b)
Regularized solution(c)
Absolute errorFrom Figures3 and 4, we can see the proposed method still yields relatively accurate numerical solutions for a two-dimensional example.Remark 19.
As we know, we need to compute the gradient of an abstract operator by the Levenberg-Marquardt method, and we usually use numerical differential method to approximately substitute the gradient when it is hard to obtain the exact expression for the gradient. So, we need to introduce a numerical differentiation stepτ. However, a lot of numerical examples are sensitive to τ. In this computation, we have successfully avoided the problem of selection of the differential step. This is because that the operator F is linear to f and the gradient in (66) can be replaced by (67) where the gradient is not dependent to τ. But when we calculate (66) by a finite difference, the initial and boundary conditions need to be homogeneous.


## 7. Conclusions

In this paper, we investigate an inverse space-dependent source problem of a multiterm time-fractional diffusion equation with nonhomogeneous boundary condition in a general domain. The existence and uniqueness of a weak solution for the direct problem are obtained by the Lax-Milgram theorem. Then, the uniqueness for the inverse problem is provided by using the analytic continuation and the Laplace transformation. Finally, we employ the Levenberg-Marquardt regularization method to find an approximation of the source function.In fact, we have obtained the well-posedness of the direct problem and also proved the uniqueness of the inverse problem for the case of nonhomogeneous boundary in the present paper. But there are many other important problems still remain open. For theoretical aspect, the conditional stability of the inverse source problem is still not obtained. For numerical computation, on the one hand, we have indeed improved the original M-L algorithm. However, some nonsmooth source functions, e.g., discontinuous function, are not reconstructed well although the uniqueness holds true theoretically.



---

*Source: 1825235-2020-09-19.xml*