A Two Dimensional Discrete Mollification Operator and the Numerical Solution of an Inverse Source Problem

We consider a two-dimensional time fractional diffusion equation and address the important inverse problem consisting of the identification of an ingredient in the source term. The fractional derivative is in the sense of Caputo. The necessary regularization procedure is provided by a two-dimensional discrete mollification operator. Convergence results and illustrative numerical examples are included.


Introduction
The mollification method has been known since the 1980s as a regularization procedure for ill-posed problems [1,2].Later on, it was implemented as an accelerator of explicit schemes [2][3][4] and there are a variety of implementations and variations.Among the most recent, Reference [5] presents a new and very close method, called the convolution regularization method.
Applications of one-dimensional mollification methods are found, for instance, in References [6][7][8].Mollification in two dimensions is less common and it is worth mentioning that it serves as a smoothing operator in Reference [9], appears as an accelerator of explicit schemes in Reference [4], and is a regularization tool for an identification problem in Reference [10].
In this paper, we deal with the two-dimensional discrete mollification operator defined in Section 3.2 of Reference [4].This operator, which we call a 2DDM operator, will be defined shortly.Based on this operator, we propose the stable solution of a challenging inverse problem consisting of the identification of an ingredient in a source term for a time fractional diffusion equation, where the fractional derivative is defined in the sense of Caputo.
This inverse problem was solved in References [11][12][13].The first two papers implemented the Tikhonov method as the required regularization procedure.The third one introduced a modified quasi-boundary regularization method.Our main original contribution is the implementation of two-dimensional discrete mollification as a regularization method for the numerical solution of this inverse source problem.It is important to have alternatives and this is the case with regularization.We offer an alternative that has not been considered before for the stable solution of this inverse problem.
The solution of the inverse problem requires the solution of the corresponding direct problem.In Reference [12], there are no details on the numerical solution, although the authors provide abundant numerical results.In Reference [11], the numerical solution is obtained by a finite difference scheme, but the only examples are one-dimensional.In our numerical solution of the direct problem, we make use of the implicit scheme defined in Reference [14] for a more general equation.
Problems based on time fractional diffusion equations are useful in the description of anomalous diffusion in highly heterogeneous porous media.There are a variety of time fractional inverse source problems.According to the number of fractional differentiation orders, they are classified as single-term and multi-term problems.This paper deals with a single-term problem that is ill-posed in the Hadamard sense.However, there are other situations, for instance in Reference [15], where the problem is single-term and depends on a parameter.For some values of this parameter, it is a well-posed problem in the Hadamard sense.Finally, we mention that in Reference [16], there is a theoretical analysis and some numerical results for a fairly general multiterm time fractional inverse source problem.
The paper is organized as follows: Section 2 is devoted to the 2DDM operator.Section 3 deals with the convergence analysis for the solution of the inverse problem, and Section 4 includes the numerical experiments.After this section, there are some final comments and conclusions.

One-Dimensional Discrete Mollification
The one-dimensional mollification of a discrete function [1,2,4] is defined as follows: Let y = {y j } j∈Z be a discrete function, which can, for example, be the evaluations or cell averages of a real function y = y(x) at equidistant grid points on: The discrete mollification of y, denoted J η y, is: where η is the integer support parameter.

Mollification Weights
The weights are defined as integrals over cells.Given positive real numbers p and h, compute: This is the continuous mollification parameter.The weights are obtained by numerical integration of the truncated Gaussian kernel: where the normalization constant: is chosen in such a way that R κ pδ = 1.Thus, the mollification weights are obtained by integrations of this kernel over cells, namely, where: We usually take p = 3. Chosen in this way, the weights w i satisfy: The set of weights has a feature worth mentioning: Independence of the step size h.
then the w i s depend only on p and η.Some weights are shown in Table 1.

Two-Dimensional Discrete Mollification
Likewise [4], two-dimensional discrete mollification requires a uniform bidimensional grid of the form (x i , y j ) where x i = x 0 + ih x and y j = y 0 + jh y .Let Y be a discrete function defined by Y(x i , y j ) = Y ij .For each positive integer η, the two dimensional discrete mollification of Y is: The weights of this two dimensional convolution are denoted by w km = w k w m .If h x = h y = h and δ is given by Equation (2), the shape of a weight is shown in Figure 1.
The main relation with one dimensional mollification is [4]:

Abstract Setting
The discrete convolutions that define mollification in one and two dimensions have their counterparts in infinite dimensions.It is true that in numerical computations, one deals with discrete functions, but what about the error estimates?They relate the abstract model based on differential equations with the discretized problem.
Convolutions are troublesome near borders; thus, whenever a convolution is mentioned, some kind of domain restriction is in order.Our presentation of the abstract setting follows Reference [17].Let If s = (s 1 , s 2 ), the two-dimensional truncated Gaussian kernel is defined by multiplication of the known one dimensional kernels, that is: where . This kernel satisfies: The next lemma establishes the consistency and stability of the two-dimensional discrete mollification.
Lemma 1. Suppose g is a Lipschitz continuous function with Lipschitz constant L defined on I and G = G ij is the associated discrete function defined on a uniform grid covering I δ with step size h in both directions.Then, there exists a constant C such that: Proof.This is Lemma 3.1 Part (1) of Reference [17].
Remark 1.Notice that due to Equation (2), we may say that the right hand side bound is of the form C δ, where C depends on L, p, and η.
2. Stability: If g is a Lipschitz continuous approximation of g so that g − g ∞,I ≤ , then: Remark 2. Actually, the right hand side bound is of the form C(h + δ + ), according to the previous remark.
Proof.It follows from Part 1 that since: where µ stands for measure.Notice that this estimate is valid in any other square subset of R 2 .Only the constant may change.

The Inverse Source Problem
Fractional derivatives have been known since the seventeenth century, but only recently have they become the subject of attention due in part to a variety of interesting applications.Podlubny [18] says: "Fractional derivatives provide an excellent instrument for the description of memory and hereditary properties of various materials and processes."We mention only two examples: For materials, we chose the damage evolution of brain tissue modeled in Reference [19] by a model in which the time derivative is a Caputo-Almeida fractional derivative.For processes, we selected the anomalous diffusion of a contaminant in fractured porous rocks modeled in Reference [20].The model is based on fractional differentiation in the sense of Caputo.

Direct Problem
Our concern is a time fractional diffusion equation.More precisely, we are interested in an initial/boundary value problem, in which the time derivative is a Caputo fractional derivative of order α, 0 < α < 1.The Caputo fractional derivative of order α of a differentiable function v is: where Γ is the Gamma function given by: We are interested in the following initial/boundary value problem: where Ω ⊂ R d and: The coefficients satisfy . The theoretical analysis that follows is from References [11,13,21].The operator −L is a symmetric uniformly elliptic operator defined on D(−L) = H 2 (Ω) ∩ H 1 0 (Ω).The direct problem is to find u that safisfies Equation ( 13), with p and f known functions and α a known fractional order of differentiation.If we think of this problem as the model for a contaminant diffusion in groundwater, then f (z) plays the role of contaminant discharge intensity and p(t) is an attenuation coefficient. Let The first step is the consideration of a fractional power operator (−L) γ , γ ∈ R and the space: where (•, •) is the inner product in L 2 (Ω).It is a Hilbert space with the norm: The Mittag-Leffler function is a generalization of the exponential function and is very useful when dealing with fractional derivatives.Definition 2. A two-parameter function of the Mittag-Leffler type is defined by: where α > 0, β > 0 are real constants and Γ is the Gamma function defined by Equation (12).
The formal solution of the direct problem ( 13) can be written in terms of a Mittag-Leffler function.It is: where From Equation (18): Denote: Thus: and its Fourier coefficients are q n = f n Q n (T).We conclude f n = q n Q n (T) .

Inverse Problem
We are interested in the inverse problem of determining the unknown source ingredient f (z) in Equation ( 13) based on the overposed data u(z, T) = q(z).Time T is a final time and the data q(t) are not known exactly; what is known is a measurement q so that q − q ∞ < .
As in Reference [11], the identification of f can be formulated as the operator equation: The operator K is a linear self-adjoint compact operator and Problem ( 19) is ill-posed.
Our aim is to obtain a function f δ that minimizes the functional: where J 2 η q is the mollification of a discrete version of the noisy data q .Since K is a self-adjoint compact operator, there exists a minimizer of H, f δ , given by: Then, we have: where . Moreover, suppose q − q ∞ < and there exist m ∈ R + and E, a constant such that f ≤ E, then: where C 1 depends on α, T, λ 1 , µ(Ω). Proof. .
The inequality follows from Theorem 1.However, , E}, we obtain: Following Reference [11] and taking into consideration Lemma 1, Part 1, we conclude:

Numerical Results
There are many alternatives for the numerical solution of a time fractional diffusion problem.If we restrict our attention to methods based on finite difference discretizations, we mention Reference [22], in which Caputo fractional derivatives are computed by the so-called Caputo fractional derivative rule, which in turn depends on a modified trapezoidal rule.This method is successfully implemented in Reference [23] for plane strain/stress elasticity in the presence of nonlocal elasticity.
Our finite difference method for the solution of the time fractional diffusion equation is taken from Reference [14], which in turn is inspired by Reference [24].The method is implicit and its properties about stability and convergence are discussed in Reference [14].
We consider the following time fractional diffusion problem: where: After discretization, operator L leads to a matrix l, whose eigenvectors and eigenvalues play the roles of the eigenfunctions and eigenvalues of the differential operator.
The Q n s are approximated according to Reference [11] by: The integrals are approximated by composite trapezoidal rule.The errors are measured by the root mean square (rms) norm given by: The particular examples considered here are for: The space step size h = 1/32 is used for both coordinates, x and y.Final time is T = 1 and ∆t = 1/10.

Example 1
Let a 11 (x, y) = x 2 + 1, a 22 (x, y) = y 2 + 1.The relative l 2 norms of the identification are shown in Tables 2 and 3. Let a 11 (x, y) ≡ 1 and a 22 (x, y) ≡ 1. Relative l 2 errors for the identification of the source term are shown in Tables 4 and 5. Figure 2 gives an idea of the quality of the identification.

Discussion and Final Comments
The 2DDM operator is a useful regularization procedure when dealing with inverse source problems for time fractional diffusion equations.It is also appropriate for surface and gradient fitting as shown in Reference [17], but we do not pursue this objective here.There are more applications of this operator and we hope to master some of them; for instance, we would like to modify our identification problem and deal with the factor p(t) as unknown and the factor f (x, y) as data.