A Point Source Identification Problem for a Time Fractional Diffusion Equation

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.

1. Introduction

Let Ω be a bounded domain in ℝ² and let 𝜕Ω be the boundary of Ω. Consider the following time fractional diffusion process: ₀𝐷^𝛾_𝑡𝑢()()(𝑥)(),𝑢()()()𝑥,𝑡−ℒ𝑢𝑥,𝑡=𝑓,𝑥∈Ω,𝑡∈0,𝑇𝑥,0=0,𝑥∈Ω,ℬ𝑢𝑥,𝑡∶=𝜇𝑢𝑥,𝑡+𝛽𝜕𝑢(),𝜕𝑣=0,𝑥∈𝜕Ω,𝑡∈0,𝑇 (1) where ℒ is the uniformly elliptic operator, 𝜈 is the outward normal at the boundary 𝜕Ω, and 𝜇, 𝛽 are known constants which are not simultaneously zero. Here, ₀𝐷^𝛾_𝑡 stands for the Caputo fractional derivative operator of order 0<𝛾≤1 defined by ₀𝐷^𝛾_𝑡𝜓(𝑡)1∶=Γ1−𝛾^𝑡₀𝜓^′(𝑠)(𝑡−𝑠)^𝛾𝜓𝑑𝑠,0<𝛾<1_′(𝑡),𝛾=1, (2) 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, 𝑢(𝑥,0)=𝜑(𝑥), we are always able to simplify the system (1) into two components; that is, 𝑢=𝑣+𝑤, where 𝑣 solves the homogeneous equation with nonhomogeneous initial condition and 𝑤 satisfies the nonhomogeneous equation with homogeneous initial condition. As we know, the initial value/boundary value problem associated with 𝑣 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 𝑓 is a known function, we are asked to determine the solution function 𝑢(𝑥,𝑡) so as to satisfy (1). So posed, this is a direct problem. However, the source term 𝑓 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 𝑡=𝑇 as 𝑢(𝑥,𝑇)+𝜖∗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 𝑡 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 𝑓(𝑥,𝑡). 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 Section 2, 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 function 𝑓(𝑥) is of the form 𝑓(𝑥)=^𝑀₁𝜆_𝑗𝛿𝑥−𝜃_𝑗,𝑗=1,2,…,𝑀, (3) where 𝜃_𝑗 denotes the location of the point source and 𝜆_𝑗 is the intensity associated with each point source at 𝜃_𝑗. Thereby, the temperature distribution 𝑢=𝑢(𝑥,𝑡) inside the domain Ω is generated by 𝑓(𝑥) satisfying ₀𝐷^𝛾_𝑡𝑢(𝑥,𝑡)−ℒ𝑢(𝑥,𝑡)=^𝑀_𝑗=1𝜆_𝑗𝛿𝑥−𝜃_𝑗(),𝑢()()(),𝑥∈Ω,𝑡∈0,𝑇𝑥,0=0,𝑥∈Ω,ℬ𝑢𝑥,𝑡∶=𝜇𝑢𝑥,𝑡+𝛽𝜕𝑢(),𝜕𝜈=0,𝑥∈𝜕Ω,𝑡∈0,𝑇 (4) where 𝛿(⋅) is the Dirac delta function. Meanwhile, let 𝑁 be a natural number and let {𝑥_𝑖}^𝑁_𝑖=1 be a group of points in Ω. Here, the points 𝑥_𝑖, 𝑖=1,…,𝑁, scattered in Ω are the collocation points. Our goal is to determine the strength sources 𝜆_𝑗 and the locations 𝜃_𝑗 from user-input estimated position and the set of final measurement data 𝑢^𝜖𝑥_𝑖,𝑇∶=𝑢_𝑖,𝑇+𝜐_𝑖,𝑥_𝑖∈Ω,𝑖=1,2,…,𝑁, (5) where (𝜐_𝑖)^𝑁_𝑖=1 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 {𝜃_𝑗} are given. Under this assumption, we come to the problem of the recovery of the intensity 𝜆_𝑗 associated with the point sources 𝜃_𝑗 from the 𝑁 distinct final collocation data 𝑢^𝜖(𝑥_𝑖,𝑇). This recovery problem is ill-posed, which prompts us to use some regularization methods.

Consequently, we assume that the locations of the point sources 𝜃_𝑗 are not known but an initial guess location 𝜃_𝑗 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, 𝜃_𝑗∈𝐵𝜃_𝑗,𝜌_𝑗𝐵∩Ω,𝑗=1,2,…,𝑀,𝜃_𝑗,𝜌_𝑗∩𝐵𝜃_𝑘,𝜌_𝑘=∅,1≤𝑗<𝑘≤𝑀, (6) where 𝐵(𝜃,𝜌) 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 function 𝐺(𝑥,𝑡;𝜃) can be defined as the solution of the following equations: ₀𝐷^𝛾_𝑡𝐺()()()(),𝐺()()()𝑥,𝑡;𝜃−ℒ𝐺𝑥,𝑡;𝜃=𝛿𝑥−𝜃,𝑥,𝜃∈Ω,𝑡∈0,𝑇𝑥,0;𝜃=0,𝑥,𝜃∈Ω,ℬ𝐺𝑥,𝑡;𝜃∶=𝜇𝐺𝑥,𝑡;𝜃+𝛽𝜕𝐺().𝜕𝜈=0,𝑥,𝜃∈𝜕Ω,𝑡∈0,𝑇 (7) By applying Laplace transform technique, we have that 𝐺(𝑥,𝑡;𝜃)=^∞_𝑛=1𝑡^𝛾𝐸_𝛾,𝛾+1−𝜆²_𝑛𝑡^𝛾𝜑_𝑛(𝜃)𝜑_𝑛(𝑥), (8) where 𝜑_𝑛 is the 𝑛th orthonormal eigenfunction and 𝜆_𝑛 is the corresponding eigenvalue to the Sturm-Liouville problem ℒ𝜑+𝜆²_𝑛𝜑=0,𝜇𝜑+𝛽𝜕𝜑𝜕𝜈=0, (9) and 𝐸_𝛾,𝜁(𝑧) is the Mittag-Leffler function defined by 𝐸_𝛾,𝜁(𝑧)=^∞_𝑘=0𝑧^𝑘Γ.𝛾𝑘+𝜁 (10) 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 𝑢(𝑥,𝑡)=^𝑀_𝑗=1𝜆_𝑗𝐺𝑥,𝑡;𝜃_𝑗. (11) Therefore, when the locations {𝜃_𝑗} of point sources are known, once we obtain the final time measurement data specified in (5), we can solve the following linear algebraic equation: ^𝑀_𝑗=1𝜆_𝑗𝐺𝑥_𝑖,𝑇;𝜃_𝑗=𝑢^𝜖_𝑖,𝑇,𝑖=1,2,…,𝑁 (12) to get the unknown values of the intensities 𝜆_𝑗. Moreover, denoting 𝜆=(𝜆_𝑗) and 𝜛=(𝑢^𝜀_𝑖,𝑇), (12) can be rewritten as the following matrix form: 𝐴𝜆=𝜛, (13) where 𝐴 is an 𝑁×𝑀 matrix: 𝐴_𝑖𝑗=𝐺𝑥_𝑖,𝑇;𝜃_𝑗. (14) Taking 𝑁=𝑀, then the system of (13) contains 𝑁 linear equations with 𝑁 unknowns. Subsequently, if the matrix 𝐴 is invertible, one simply has 𝜆=𝐴⁻¹𝜛. (15) 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: min_𝜆𝐽_𝛼(𝜆)∶=min_𝜆‖𝐴𝜆−𝜛‖²+𝛼²‖𝜆‖², (16) where 𝛼>0 is the regularization parameter and ‖⋅‖ denotes the usual Euclidean norm. It is well known [34] that the minimal element of 𝐽_𝛼 can be written as 𝜆_𝛼=𝐴^∗𝐴+𝛼²𝐼⁻¹𝐴^∗𝜛, (17) where 𝐴^∗ denotes the conjugate transposed matrix of 𝐴 and 𝐼 denotes the identity matrix.

Next, assume that the locations {𝜃_𝑗} 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 {𝜃_𝑗}, we define the following union set Θ∶=^𝑀_𝑗=1𝐵𝜃_𝑗,𝜌_𝑗∩Ω (18) and suppose that it contains all exact positions of the point sources with proper radius 𝜌_𝑗. To linearize the nonlinear system, we take some additional collocation points {𝜉_𝑙}^𝑛_𝑙=1 from the set Θ. Assume that {𝜉_𝑙}^𝑛_𝑙=1 are uniformly distributed in Θ. On each point 𝜉_𝑙, we put in a point source with intensity 𝜏_𝑙. Suppose that the temperature distribution generated by the 𝑛 point sources {𝜉_𝑙}^𝑛_𝑙=1 is equal to that generated by the 𝑀 point sources {𝜃_𝑗}^𝑀_𝑗=1. Subsequently, we have ₀𝐷^𝛾_𝑡𝑢(𝑥,𝑡)−ℒ𝑢(𝑥,𝑡)=^𝑛_𝑙=1𝜏_𝑙𝛿𝑥−𝜉_𝑙,(),𝑢()()𝑥∈Ω,𝑡∈0,𝑇𝑥,0=0,𝑥∈Ω,𝜇𝑢𝑥,𝑡+𝛽𝜕𝑢(),𝜕𝜈=0,𝑥∈𝜕Ω,𝑡∈0,𝑇 (19) with additional data (5), where 𝜏_𝑙 is the intensity at the location 𝜉_𝑙. By using the above proposed method to solve (19) with (5), the intensity 𝜏_𝑙 of each point source 𝜉_𝑙 can then be obtained approximately. Next, we transform the intensities {𝜏_𝑙}⊂𝐵(𝜃_𝑗,𝜌_𝑗)∩Ω back to a single source point as follows: the 𝑗th unknown source intensity 𝜆_𝑗 associated with each ball 𝐵(𝜃_𝑗,𝜌_𝑗)∩Ω is approximated by 𝜆_𝑗 as 𝜆_𝑗∶=_{𝑙∶𝜉𝑙∈𝐵(𝜃𝑗,𝜌𝑗)∩Ω}𝜏_𝑙,𝑗=1,2,…,𝑀. (20) 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 𝜃_𝑗, we use the weight sum of the location coordinate 𝜉_𝑙 in the ball 𝐵(𝜃_𝑗,𝜌_𝑗) to approximate the exact location. More specifically, the approximation location 𝜃_𝑗 corresponding to the intensity 𝜆_𝑗 is defined by 𝜃_𝑗1∶=𝜆_{𝑗𝑙∶𝜉𝑙∈𝐵(𝜃𝑗,𝜌𝑗)∩Ω}𝜏_𝑙𝜉_𝑙,𝑗=1,2,…,𝑀. (21)

4. Numerical Examples

In this section, some numerical examples are given to verify the effectiveness of the method proposed in Section 3. 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 RMS=1𝑀^𝑀_𝑗=1𝜆_𝛼,𝑗−𝜆_𝑗². (22) The noisy data 𝑢^𝜖(𝑥_𝑖,𝑇) at measurement points 𝑥_𝑖 is obtained by adding random noise to the exact data 𝑢(𝑥_𝑖,𝑇) by 𝑢^𝜖𝑥_𝑖,𝑇=𝑢𝑥_𝑖,𝑇+𝜖rand(𝑖) (23) for 𝑥_𝑖∈Ω, where rand(𝑖) is a random number between [−1,1]. The measurement points {𝑥_𝑖} 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 𝐿-curve criterion [35–37] to choose the regularization parameter.

Example 1. Consider the following heat conduction problem on a semi-infinite stripe domain Ω={(𝑥,𝑦)∣0≤𝑥≤𝑙,𝑦≥0}: 𝑢_𝑡(),𝑢𝑢𝑥,𝑦,𝑡−Δ𝑢=𝑓𝑥,𝑦,0<𝑥<𝑙,𝑦>0,𝑡∈0,𝑇𝑥,𝑦,0=0,0≤𝑥≤𝑙,𝑦≥0,_𝑥=0=𝑢_𝑥=𝑙𝑢|=0,𝑦≥0,_𝑦=0=0,0≤𝑥≤𝑙. (24) Green’s function is given by 𝐺=2𝑥,𝑦,𝑡;𝜉,𝜂𝑙√𝑒𝜋𝑡^{−(𝑦2+𝜂2)/4𝑡}sinh𝑦𝜂×2𝑡_∞𝑛=1𝑒^{−𝑛2𝜋2𝑡/𝑙2}sin𝑛𝜋𝜉𝑙sin𝑛𝜋𝑥𝑙𝐻(𝑡), (25) where 𝐻(⋅) is the Heaviside function. Without loss of generality, we take 𝑙=1.

In this test, we consider the case that the source function (3) contains five source points {𝜃_𝑗}⁵_𝑗=1. The input source locations 𝜃_𝑗 are randomly chosen such that 𝜃_𝑗−𝜃_𝑗<𝜌_𝑗,for𝑗=1,2,3,4,5, (26) where 𝜌_𝑗 are sufficient small to ensure 𝜃_𝑗∈Ω. The value of the parameters {𝜌_𝑗}⁵_𝑗=1 in (18) used in this computations is 0.1. When the noisy data are given in the final time 𝑇=1, we demonstrate the numerical performance under two noise levels: 𝜖=0.01 and 𝜖=0.001. The computations are performed by using a total of 𝑄 trial centers in each ball. We report the numerical results under different 𝑄 in Tables 1 and 2. The displayed results show that the total number 𝑄 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 

Example 1: numerical comparison for 𝜖=0.01 and 𝑀=5.

Table 2 

Example 1: numerical comparison for 𝜖=0.001 and 𝑀=5.

Example 2. In this example, we consider the following inverse identification problem on a square domain Ω=[0,1]×[0,1] for 0<𝛾<1: ₀𝐷^𝛾_𝑡𝑢,(),𝑢𝑢().𝑥,𝑦,𝑡−Δ𝑢=𝑓𝑥,𝑦𝑥,𝑦∈Ω,𝑡∈0,𝑇𝑥,𝑦,0=0,𝑥,𝑦∈Ω,𝑥,𝑦,𝑡=0,𝑥,𝑦∈𝜕Ω,𝑡∈0,𝑇 (27) By virtue of Laplace’s transform [2, 38], one can derive the corresponding Green’s function 𝐺𝑥,𝑦,𝑡;𝜉,𝜂=4^∞_𝑛=1𝑡^𝛾𝐸_𝛾,𝛾+1−𝑛²𝜋²𝑡^𝛾()().sin𝑛𝜋𝑥×sin𝑛𝜋𝑦sin𝑛𝜋𝜉sin𝑛𝜋𝜂 (28)

Firstly, we see the robustness of the proposed algorithm about the parameter 𝛾. For 𝜖=0.01,0.001, 𝜌_𝑗=0.1, and 𝑇=1, we report the RMS in Table 3 under different 𝛾 for three point sources located at (0.3,0.3),(0.5,0.7),(0.7,0.3) (29) 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 

Example 2: RMS under different 𝜖 and 𝛾 for 𝑀=3.

Table 4 

Example 2: the approximation locations for (0.3, 0.3), (0.5, 0.7), and (0.7, 0.3).

Table 5 

Example 2: the approximation locations for (0.3, 0.3), (0.5, 0.7), and (0.7, 0.3).

Secondly, we also consider the effect of the final time 𝑇 on the numerical precision. Fixing 𝛾=0.5 and choosing parameter 𝜌_𝑗=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 𝑇. Such phenomenon can be explained by the nature of ill-posed inverse source identification problem.

Table 6 

Example 2: numerical comparison for 𝜖=0.01 using different 𝑇.

Finally, using the previous point source, we plot the exact and approximation locations of source points in Figure 1 for 𝜖=0.1 and 𝑇=0.1,1. The computational intensities are 0.9326, 2.9555, and 5.1065 for 𝑇=0.1 and 1.4365, 1.9797, and 5.8621 for 𝑇=1, respectively. It can be seen that even for high noise level 𝜖=0.1, the proposed method produces an acceptable numerical approximation.

(a)

(b)

Figure 1 

Approximation for 𝜖=0.1 with 𝑇=0.1 (a), 1 (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.

Acknowledgments

This work is supported by the Fundamental Research Funds for the Central Universities (SWJTU11BR078, ZYGX2011J104) and the NSF of China (no. 11226040).