# A Point Source Identification Problem for a Time Fractional Diffusion Equation

Authors: Xiao-Mei Yang; Zhi-Liang Deng
Journal: Advances in Mathematical Physics
 (2013)

---

## Abstract

An inverse source identification problem for a time fractional diffusion equation is discussed. The unknown heat source is supposed to be space dependent only. Based on the use of Green’s function, an effective numerical algorithm is developed to recover both
the intensities and locations of unknown point sources from final measurements. Numerical results indicate that the proposed method is efficient and accurate.


---

## Body

## 1. Introduction

LetΩ be a bounded domain in ℝ2 and let ∂Ω be the boundary of Ω. Consider the following time fractional diffusion process:
(1)0Dtγu(x,t)-ℒu(x,t)=f(x),x∈Ω,t∈(0,T),u(x,0)=0,x∈Ω,ℬu(x,t)∶=μu(x,t)+β∂u∂v=0,x∈∂Ω,t∈(0,T),
where ℒ is the uniformly elliptic operator, ν is the outward normal at the boundary ∂Ω, and μ, β are known constants which are not simultaneously zero. Here, 0Dtγ stands for the Caputo fractional derivative operator of order 0<γ≤1 defined by
(2)0Dtγψ(t)∶={1Γ(1-γ)∫0tψ′(s)(t-s)γds,0<γ<1ψ′(t),γ=1,
where Γ(·) is the standard Γ-function and the prime denotes the general derivative.From the last few decades, fractional calculus grabbed great attention of not only mathematicians and engineers but also many scientists from all fields (e.g., see [1–4]). Fractional diffusion equations describe anomalous diffusions on fractals (physical objects of fractional dimension, like some amorphous semiconductors or strongly porous materials; see [5, 6] and references therein). Indeed, fractional derivatives provide an excellent tool for the description of memory and hereditary properties of various materials and processes. This is the main superiority of fractional derivatives in comparison with classical integer-order models, in which such effects are in fact neglected. For the detailed theory and application of fractional calculus, one can refer to [1–4] and references therein. Not only have differential equations of fractional order attracted people’s attention, but also theories and applications related to physics and geometry of fractal dimension have been well studied (e.g., [7–11]).If the initial condition is nonhomogeneous, that is,u(x,0)=φ(x), we are always able to simplify the system (1) into two components; that is, u=v+w, where v solves the homogeneous equation with nonhomogeneous initial condition and w satisfies the nonhomogeneous equation with homogeneous initial condition. As we know, the initial value/boundary value problem associated with v is well-posed and there exist many works on such forward problem, for example, [12, 13]. In the following, instead of nonhomogeneous initial condition, we only focus on the system (1) with homogeneous initial condition. Ordinarily, when f is a known function, we are asked to determine the solution function u(x,t) so as to satisfy (1). So posed, this is a direct problem. However, the source term f is not always known and has to be computed from some additional data. The additional information is mainly the following: the interior/boundary transient measurement values and the final measurement values. Here, we suppose that the measured data are given in final time t=T as u(x,T)+ϵ*rand (ϵ is the noise level).For most classical partial differential equations, the identification and reconstruction of source functions from the final data or the partial boundary data are an inverse problem with many applications (e.g., [14]). A number of articles address the solvability problem of source term identification. For parabolic-type differential equation, please see [15–26]. For elliptic-type differential equation, one can refer to [27–29], though the source identification problem has been well discussed in the classic framework, yet, to the best of the authors’ knowledge, there are rare researches in the aspect of the source identification problem associated with fractional differential equation in spite of the physical and practical importance. As indicated in [30–33], the source identification problem associated with the time fractional diffusion equation is also ill-posed. That means the solution does not depend continuously on the given data and any small perturbation in the given data may cause large change to the solution. In [33], when additional data is given on the partial boundary, the uniqueness in identifying a source term independent of time t is established for one-dimensional time fractional diffusion equation. In [30], if the final time temperature distribution is known, the existence and uniqueness results are proved. Murio and Mejía [31] propose a mollification regularization technique to reconstruct the unknown forcing term f(x,t). In this paper, we aim to deal with the special case that the sources and measurements are both point like. The main focus will be placed on the recovery of both intensities and locations of the unknown point source term. For this, we propose a method based on the use of Green’s function to solve the inverse source identification problems.The outline of the paper is as follows. In Section2, we provide a brief sketch on the considered identification problem. The reconstruction method by Green’s function is then given in Section 3. Numerical implementation of the proposed method is provided in Section 4. In Section 5, we summarize the results.


## 2. Statement of the Problem

In this paper, we deal with the special case that the source functionf(x) is of the form
(3)f(x)=∑1M‍λjδ(x-θj),j=1,2,…,M,
where θj denotes the location of the point source and λj is the intensity associated with each point source at θj. Thereby, the temperature distribution u=u(x,t) inside the domain Ω is generated by f(x) satisfying
(4)0Dtγu(x,t)-ℒu(x,t)=∑j=1Mλjδ(x-θj),x∈Ω,t∈(0,T),u(x,0)=0,x∈Ω,ℬu(x,t)∶=μu(x,t)+β∂u∂ν=0,x∈∂Ω,t∈(0,T),
where δ(·) is the Dirac delta function. Meanwhile, let N be a natural number and let {xi}i=1N be a group of points in Ω. Here, the points xi, i=1,…,N, scattered in Ω are the collocation points. Our goal is to determine the strength sources λj and the locations θj from user-input estimated position and the set of final measurement data
(5)uϵ(xi,T)∶=ui,T+υi,xi∈Ω,i=1,2,…,N,
where (υi)i=1N denotes the Gaussian variable with mean zero and variance ϵ. This magnitude ϵ also represents the level of noises.Let us first suppose that the locations of the point sources{θj} are given. Under this assumption, we come to the problem of the recovery of the intensity λj associated with the point sources θj from the N distinct final collocation data uϵ(xi,T). This recovery problem is ill-posed, which prompts us to use some regularization methods.Consequently, we assume that the locations of the point sourcesθj are not known but an initial guess location θ-j is given for each unknown point source. Moreover, we make the assumption that each point source belongs to a distinct ball inside the domain; that is,
(6)θj∈B(θ-j,ρj)∩Ω,j=1,2,…,M,B(θ-j,ρj)∩B(θ-k,ρk)=∅,1≤j<k≤M,
where B(θ,ρ) denotes the ball centered at θ with radius ρ. It should be pointed that if two or more point sources are concentrated in a sufficient small domain, the proposed method in the following section will treat them as one point source.


## 3. Methodology Based on Green’s Function

In this section, we discuss the identification method based on Green’s function. Green’s functionG(x,t;θ) can be defined as the solution of the following equations:
(7)0DtγG(x,t;θ)-ℒG(x,t;θ)=δ(x-θ),x,θ∈Ω,t∈(0,T),G(x,0;θ)=0,x,θ∈Ω,ℬG(x,t;θ)∶=μG(x,t;θ)+β∂G∂ν=0,x,θ∈∂Ω,t∈(0,T).
By applying Laplace transform technique, we have that
(8)G(x,t;θ)=∑n=1∞tγEγ,γ+1(-λn2tγ)φn(θ)φn(x),
where φn is the nth orthonormal eigenfunction and λn is the corresponding eigenvalue to the Sturm-Liouville problem
(9)ℒφ+λn2φ=0,μφ+β∂φ∂ν=0,
and Eγ,ζ(z) is the Mittag-Leffler function defined by
(10)Eγ,ζ(z)=∑k=0∞zkΓ(γk+ζ).
For the details of Mittag-Leffler function, one can refer to [2].Utilizing Green’s function, we then can write the solution of (4) as
(11)u(x,t)=∑j=1MλjG(x,t;θj).
Therefore, when the locations {θj} of point sources are known, once we obtain the final time measurement data specified in (5), we can solve the following linear algebraic equation:
(12)∑j=1MλjG(xi,T;θj)=ui,Tϵ,i=1,2,…,N
to get the unknown values of the intensities λj. Moreover, denoting λ=(λj) and ϖ=(ui,Tε), (12) can be rewritten as the following matrix form:
(13)Aλ=ϖ,
where A is an N×M matrix:
(14)Aij=G(xi,T;θj).
Taking N=M, then the system of (13) contains N linear equations with N unknowns. Subsequently, if the matrix A is invertible, one simply has
(15)λ=A-1ϖ.
However, due to the ill-posedness of the source identification problem, the system of (13) is ill-conditioned, and hence a direct solution, as given by (15), will be either impossible or will produce very inaccurate results. To obtain stable solutions to these kinds of ill-conditioning systems, various regularization techniques have been studied and applied extensively [34]. Here, a standard Tikhonov regularization technique is adopted to find the approximation solution of the matrix equation (13). By λα, we denote the Tikhonov regularized solution defined to be the minimal element of the following least square problem:
(16)minλJα(λ)∶=minλ∥Aλ-ϖ∥2+α2∥λ∥2,
where α>0 is the regularization parameter and ∥·∥ denotes the usual Euclidean norm. It is well known [34] that the minimal element of Jα can be written as
(17)λα=(A*A+α2I)-1A*ϖ,
where A* denotes the conjugate transposed matrix of A and I denotes the identity matrix.Next, assume that the locations{θj} of the point sources are also unknown. In such case, we will get a nonlinear system. The nonlinear system is not suitable or difficult for direct numerical computation. In order to eliminate the difficulty in implementing the numerical computation, we propose in the following to linearize the nonlinear system.For the estimate locations{θ-j}, we define the following union set
(18)Θ∶=⋃j=1MB(θ-j,ρj)∩Ω
and suppose that it contains all exact positions of the point sources with proper radius ρj. To linearize the nonlinear system, we take some additional collocation points {ξl}l=1n from the set Θ. Assume that {ξl}l=1n are uniformly distributed in Θ. On each point ξl, we put in a point source with intensity τl. Suppose that the temperature distribution generated by the n point sources {ξl}l=1n is equal to that generated by the M point sources {θj}j=1M. Subsequently, we have
(19)0Dtγu(x,t)-ℒu(x,t)=∑l=1nτlδ(x-ξl),x∈Ω,t∈(0,T),u(x,0)=0,x∈Ω,μu(x,t)+β∂u∂ν=0,x∈∂Ω,t∈(0,T),
with additional data (5), where τl is the intensity at the location ξl. By using the above proposed method to solve (19) with (5), the intensity τl of each point source ξl can then be obtained approximately. Next, we transform the intensities {τl}⊂B(θ-j,ρj)∩Ω back to a single source point as follows: the jth unknown source intensity λj associated with each ball B(θ-j,ρj)∩Ω is approximated by λ~j as
(20)λ~j∶=∑l:ξl∈B(θ-j,ρj)∩Ωτl,j=1,2,…,M.
With the approximation intensity in hand, we can start to look at how to find the locations of the point sources. For every point source θj, we use the weight sum of the location coordinate ξl in the ball B(θ-j,ρj) to approximate the exact location. More specifically, the approximation location θ~j corresponding to the intensity λj is defined by
(21)θ~j∶=1λ~j∑l:ξl∈B(θ-j,ρj)∩Ωτlξl,j=1,2,…,M.


## 4. Numerical Examples

In this section, some numerical examples are given to verify the effectiveness of the method proposed in Section3. In our computation, we use the MATLAB code developed by Hansen [35, 36] for solving the ill-conditioned system (13). To compare the accuracy of the approximation, we use the root mean square (RMS) which is defined as
(22)RMS=1M∑j=1M(λ~α,j-λj)2.
The noisy data uϵ(xi,T) at measurement points xi is obtained by adding random noise to the exact data u(xi,T) by
(23)uϵ(xi,T)=u(xi,T)+ϵrand(i)
for xi∈Ω, where rand(i) is a random number between [-1,1]. The measurement points {xi} are equally distributed in Ω. In addition, as we know, for ill-posed problem, the regularization parameter α plays an important role and hence has to be chosen appropriately. In theory, α depends on some a priori knowledge of exact solution and noise level ϵ [34]. However, in practice, the a priori knowledge and noise level may not always be known. Therefore, to compensate this lack of information for the noise level, it is necessary for us to consider some error-free parameter choice rules. Here, we adopt the L-curve criterion [35–37] to choose the regularization parameter.Example 1.
Consider the following heat conduction problem on a semi-infinite stripe domainΩ={(x,y)∣0≤x≤l,y≥0}:
(24)ut(x,y,t)-Δu=f(x,y),0<x<l,y>0,t∈(0,T),u(x,y,0)=0,0≤x≤l,y≥0,u|x=0=u|x=l=0,y≥0,u|y=0=0,0≤x≤l.
Green’s function is given by
(25)G(x,y,t;ξ,η)=2lπte-(y2+η2)/4tsinhyη2t×∑n=1∞e-n2π2t/l2sinnπξlsinnπxlH(t),
where H(·) is the Heaviside function. Without loss of generality, we take l=1.In this test, we consider the case that the source function (3) contains five source points {θj}j=15. The input source locations θ-j are randomly chosen such that
(26)|θ-j-θj|<ρ-j,forj=1,2,3,4,5,
where ρ-j are sufficient small to ensure θ-j∈Ω. The value of the parameters {ρj}j=15 in (18) used in this computations is 0.1. When the noisy data are given in the final time T=1, we demonstrate the numerical performance under two noise levels: ϵ=0.01 and ϵ=0.001. The computations are performed by using a total of Q trial centers in each ball. We report the numerical results under different Q in Tables 1 and 2. The displayed results show that the total number Q of trial centers plays no role in the convergence of the scheme. Only a small number of trial centers are sufficient to approximate the unknown source function. Therefore, we only consider the case when 100 trial centers are taken in each ball in subsequent examples.Table 1
Example1: numerical comparison for ϵ=0.01 and M=5.
Exact
λ
j
11
1
21
−5
−7
θ
j
(0.5, 1)
(0.5, 5)
(0.5, 10)
(0.5, 15)
(0.5, 25)
Q
=
25
λ
~
j
11.6978
1.0539
22.0974
−5.1806
−7.3625
θ
~
j
(0.5196, 0.9876)
(0.5349, 4.9513)
(0.5173, 10.0002)
(0.5059, 14.9989)
(0.5129, 25.0005)
Q
=
100
λ
~
j
11.7866
0.9881
22.1760
−5.1460
−7.3026
θ
~
j
(0.5226, 0.9925)
(0.5196, 4.9688)
(0.5188, 10.0002)
(0.5059, 15.0016)
(0.5102, 25.0001)
Q
=
256
λ
~
j
11.4531
1.0006
22.1702
−5.1226
−7.1781
θ
~
j
(0.5136, 0.9929)
(0.5219, 4.9689)
(0.5190, 10.0005)
(0.5027, 15.0023)
(0.5046, 24.9995)
Q
=
400
λ
~
j
11.6968
1.0962
22.1453
−5.0899
−7.1799
θ
~
j
(0.5213, 0.9977)
(0.5521, 4.9829)
(0.5186, 10.0005)
(0.5009, 15.0027)
(0.5047, 25.0008)
Q
=
625
λ
~
j
11.9726
1.0230
22.0421
−5.1289
−7.1438
θ
~
j
(0.5274, 0.9951)
(0.5309, 5.0000)
(0.5171, 10.0004)
(0.5033, 15.0033)
(0.5028, 24.9993)
Q
=
900
λ
~
j
11.5229
1.0794
22.3220
−4.9636
−7.1709
θ
~
j
(0.5163, 0.9975)
(0.5477, 4.9907)
(0.5209, 10.0004)
(0.4919, 15.0005)
(0.5035, 25.0000)Table 2
Example1: numerical comparison for ϵ=0.001 and M=5.
Exact
λ
j
5
7
2
5
7
θ
j
(0.3, 1)
(0.6, 5)
(0.5, 10)
(0.7, 15)
(0.2, 25)
Q
=
25
λ
~
j
4.9670
7.3015
1.8876
5.1820
6.9783
θ
~
j
(0.3016, 1.0113)
(0.6133, 5.0007)
(0.4816, 10.0010)
(0.7097, 14.9972)
(0.1985, 24.9974)
Q
=
100
λ
~
j
4.9553
7.1001
1.9929
5.0736
6.9830
θ
~
j
(0.3009, 1.0117)
(0.6057, 4.9999)
(0.4993, 10.0079)
(0.7061, 15.0017)
(0.1982, 24.9966)
Q
=
256
λ
~
j
4.9891
7.0950
2.0070
4.9055
6.7970
θ
~
j
(0.3041, 1.0112)
(0.6054, 4.9992)
(0.5018, 10.0075)
(0.7003, 15.0016)
(0.2031, 25.0043)
Q
=
400
λ
~
j
4.9690
7.0538
1.9133
4.9015
6.7483
θ
~
j
(0.3020, 1.0114)
(0.6041, 4.9998)
(0.4859, 10.0036)
(0.7002, 15.0014)
(0.2032, 25.0023)
Q
=
625
λ
~
j
4.9678
7.0433
1.9759
4.9636
6.6110
θ
~
j
(0.3022, 1.0115)
(0.6036, 5.0003)
(0.4966, 10.0061)
(0.7025, 15.0035)
(0.2078, 25.0015)
Q
=
900
λ
~
j
4.9803
7.0981
2.0104
5.0513
6.6630
θ
~
j
(0.3030, 1.0103)
(0.6052, 4.9991)
(0.5019, 10.0072)
(0.7047, 15.0059)
(0.2065, 25.0024)Example 2.
In this example, we consider the following inverse identification problem on a square domainΩ=[0,1]×[0,1] for 0<γ<1:
(27)0Dtγu(x,y,t)-Δu=f(x,y),(x,y)∈Ω,t∈(0,T),u(x,y,0)=0,(x,y)∈Ω,u(x,y,t)=0,(x,y)∈∂Ω,t∈(0,T).
By virtue of Laplace’s transform [2, 38], one can derive the corresponding Green’s function
(28)G(x,y,t;ξ,η)=4∑n=1∞tγEγ,γ+1(-n2π2tγ)sin(nπx)×sin(nπy)sin(nπξ)sin(nπη).Firstly, we see the robustness of the proposed algorithm about the parameterγ. For ϵ=0.01,0.001, ρj=0.1, and T=1, we report the RMS in Table 3 under different γ for three point sources located at
(29)(0.3,0.3),(0.5,0.7),(0.7,0.3)
with intensities 1, 3, and 5. The corresponding approximation locations are given in Tables 4 and 5. The displayed results show that the change of the parameter γ has little effect on the numerical computations, which reflects that the proposed method is robust about γ. On the other hand, one can see that, for smaller noise level ϵ, we obtain better numerical effect.Table 3
Example2: RMS under different ϵ and γ for M=3.
ϵ
γ
0.0909
0.1818
0.2727
0.3636
0.4545
0.5455
0.6364
0.7273
0.8182
0.9091
0.01
0.0231
0.0121
0.0388
0.0266
0.0137
0.0250
0.0184
0.0233
0.0822
0.0589
0.001
0.0083
0.0045
0.0088
0.0080
0.0089
0.0114
0.0099
0.0091
0.0096
0.0070Table 4
Example2: the approximation locations for (0.3, 0.3), (0.5, 0.7), and (0.7, 0.3).
ϵ
γ
0.0909
0.1818
0.2727
0.3636
0.4545
0.01
(0.2986, 0.2844)(0.5017, 0.7031)(0.7023, 0.2998)
(0.2983, 0.2987)(0.5013, 0.6993)(0.7018, 0.3014)
(0.3053, 0.3013)(0.5022, 0.6979)(0.7030, 0.3061)
(0.2951, 0.2875)(0.5005, 0.7029)(0.7036, 0.3016)
(0.3035, 0.3049)(0.5009, 0.6985)(0.7022, 0.3030)
0.001
(0.3002, 0.2937)(0.5010, 0.7010)(0.7024, 0.3018)
(0.3031, 0.2956)(0.5002, 0.7002)(0.7016, 0.3016)
(0.3019, 0.2951)(0.5006, 0.7003)(0.7020, 0.3017)
(0.3021, 0.2950)(0.5006, 0.7003)(0.7021, 0.3019)
(0.2996, 0.2947)(0.5008, 0.7008)(0.7024, 0.3018)Table 5
Example2: the approximation locations for (0.3, 0.3), (0.5, 0.7), and (0.7, 0.3).
ϵ
γ
0.5455
0.6364
0.7273
0.8182
0.9091
0.01
(0.3029, 0.2955)(0.5021, 0.7004)(0.7029, 0.3036)
(0.3003, 0.3007)(0.5007, 0.6983)(0.7010, 0.3029)
(0.3083, 0.3026)(0.4974, 0.6976)(0.6995, 0.3008)
(0.2957, 0.2992)(0.5050, 0.6980)(0.7018, 0.3088)
(0.2830, 0.2881)(0.4980, 0.7038)(0.7056, 0.3000)
0.001
(0.3011, 0.2940)(0.5010, 0.7006)(0.7025, 0.3022)
(0.3029, 0.2946)(0.5006, 0.7002)(0.7021, 0.3021)
(0.3017, 0.2943)(0.5009, 0.7006)(0.7021, 0.3017)
(0.2992, 0.2954)(0.5004, 0.7008)(0.7026, 0.3019)
(0.3023, 0.2972)(0.5005, 0.6996)(0.7021, 0.3024)Secondly, we also consider the effect of the final timeT on the numerical precision. Fixing γ=0.5 and choosing parameter ρj=0.1, we report the numerical results in Table 6, from which one can see that the accuracy of the approximation decreases with respect to the increase of the number of T. Such phenomenon can be explained by the nature of ill-posed inverse source identification problem.Table 6
Example2: numerical comparison for ϵ=0.01 using different T.
Exact
λ
j
1
3
5
θ
j
(0.3, 0.3)
(0.5, 0.7)
(0.7, 0.3)
T
=
0.1
λ
~
j
1.0108
3.0096
5.0584
θ
~
j
(0.2886, 0.2917)
(0.5007, 0.7032)
(0.7043, 0.3034)
T
=
1
λ
~
j
1.0039
3.0035
4.8805
θ
~
j
(0.3265, 0.3068)
(0.4966, 0.6984)
(0.6993, 0.3001)
T
=
2
λ
~
j
0.8916
2.9666
5.0708
θ
~
j
(0.2959, 0.2993)
(0.5046, 0.6980)
(0.7004, 0.3093)
T
=
5
λ
~
j
0.8560
2.9201
5.1635
θ
~
j
(0.2805, 0.3130)
(0.5134, 0.6938)
(0.7026, 0.3100)Finally, using the previous point source, we plot the exact and approximation locations of source points in Figure1 for ϵ=0.1 and T=0.1,1. The computational intensities are 0.9326, 2.9555, and 5.1065 for T=0.1 and 1.4365, 1.9797, and 5.8621 for T=1, respectively. It can be seen that even for high noise level ϵ=0.1, the proposed method produces an acceptable numerical approximation.Approximation forϵ=0.1 with T=0.1 (a), 1 (b).
(a)
(b)


## 5. Conclusion

Based on the use of Green’s function, we propose in this paper an effective numerical method to recover both the intensities and locations of point sources for a time fractional diffusion process. Some numerical results show that the proposed algorithm provides an accurate and reliable scheme.



---

*Source: 485273-2013-12-08.xml*