Simultaneous Inversion for Space-Dependent Di ff usion Coe ffi cient and Source Magnitude in the Time Fractional Di ff usion Equation

We deal with an inverse problem of simultaneously identifying the space-dependent diffusion coefficient and the source magnitude in the time fractional diffusion equation from viewpoint of numerics. Such simultaneous inversion problem is often of severe ill-posedness as compared with that of determining a single coefficient function. The forward problem is solved by employing an implicit finite difference scheme, and the inverse problem is solved by applying the homotopy regularization algorithm with Sigmoid-type homotopy parameter. The inversion solutions approximate to the exact solutions demonstrating that the proposed algorithm is efficient for simultaneous inversion problems in the fractional diffusion equation.


Introduction
Diffusion equations with integer-order derivatives have played important roles in modelling contaminants diffusion processes, which are called the classical diffusion models.In recent two decades, some researches indicated that the classical model was inadequate to simulate many real situations, where a particle plume spreads slower than predicted by the integer-order diffusion equation.See, e.g., Adams and Gelhar (1992), Berkowitz, Scher, and Silliman (2000), Giona, Gerbelli, and Roman (1992), Hatano, Y. and Hatano, N. (1998), Nigmatullin (1986), Xiong, Huang, and Huang (2006), Zhou and Selim (2003), and the references therein.Such slow diffusion is called anomalous subdiffusion, and one model equation is the time fractional diffusion equation given as where u = u(x, t) denotes the state variable at space point x and time t, and γ ∈ (0, 1) is called the fractional order of the derivative in time, D(x) is the space-dependent diffusion coefficient, and ω(t) > 0 is the attenuation factor, and g(x) is the source/sink magnitude, and ∂ γ u ∂t γ here means the Caputo's derivative: (2) See, e.g., Miller and Ross (1993), Podlubny (1999) and Kilbas, Srivastava and Trujillo (2006) for the definition and properties of the Caputo's derivative.
With more and more applications of the fractional derivatives in the applied sciences, there are quite a few research works on theoretical analysis and numerical methods for the forward problem of fractional differential equations.
On the other hand, there are only a few studies for the inverse problems of simultaneously identifying more than one different kinds of coefficients in the diffusion models, see (Bondarenko & Ivaschenko, 2009;Cheng, Nakagawa, Yamamoto, & Yamazaki, 2009;Rodrigues, Orlande, & Dulikravich, 2004), for instance.Rodrigues et al. (2004) considered such simultaneous inversion problem but in the 1-D integer-order diffusion equation by the conjugate gradient method; Cheng et al. (2009) studied an inverse problem of simultaneously determining the fractional order and space-dependent diffusion coefficient in the time fractional diffusion equation with the additional boundary data at one endpoint, they proved uniqueness of the simultaneous inverse problem by utilizing Gel'fand-Levitan theory and eigenfunction expansion; Bondarenko and Ivaschenko (2009) considered the same simultaneous inversion problem as in Cheng et al. (2009) but with numerical method, they gave a weighted difference scheme which is of conditional stability and convergence for solving the forward problem with zero Dirichelet boundary conditions, and presented two kinds of numerical examples by using Levenberg-Marquart algorithm.
In this article, we will consider the simultaneous inversion problem for determining the space-dependent diffusion coefficient and the source magnitude function in Equation ( 1) with final observations.The simultaneous inversion problem here could become much more difficult than those given in (Bondarenko & Ivaschenko, 2009;Cheng et al., 2009) because two different kinds of coefficient functions D(x) and g(x) in Equation ( 1) need to be determined by one kind of additional information.Another different point comes from that we will perform numerical inversions by utilizing the homotopy regularization algorithm, which is a combination of the optimal perturbation algorithm with the homotopy method.
As we know, the optimal perturbation algorithm has been verified to be efficient for solving coefficient inverse problems not only for the integer-order diffusion equations (see Li, Cheng, Yao, Liu, & Liu, 2007;Li, Yao, Jiang, & Jia, 2011;Li, Tan, Yao, Wang, & Liu, 2008;Liu et al., 2007;Su, 1995, e.g.), but also for the fractional diffusion equations (see Chi, Li, & Jia, 2011;Li, Gu, & Jia, 2012;Wei, Chen, Sun, & Li, 2010, e.g.).Wei et al. (2010) considered an inverse source problem in the spatial fractional diffusion equation, and solved the inverse problem numerically by the optimal perturbation method; Chi, Li, and Jia (2011) considered an inverse problem in the space fractional advection dispersion equation, and presented numerical inversions for determining the space-dependent source magnitude by applying the same inversion algorithm with random noisy data; Li et al. (2012) considered an inverse problem of determining the space-dependent diffusion coefficient in Equation (1), and performed numerical inversions also by applying the optimal perturbation algorithm with the regularization parameter chosen as empirical constant.However, if dealing with more complicated identification problems, such as simultaneous inversion problems, more fine regularization parameters need to be chosen to overcome the illposedness of the simultaneous inversion.In other words, we should seek to more effective algorithms to deal with simultaneous inversion problems arising in the fractional differential equations.
Fortunately, the well-known homotopy method is always utilized to approximate to a solution within a large range of convergence in many branches of mathematics, including algebra equations, mathematical programming, and differential equations, etc. Especially for solving of differential equations with analyticity, the homotopy perturbation method, which is a combination of the homotopy idea with the perturbation analysis method, has been widely studied and applied to solve nonlinear problems arising in PDE (see He, 2003He, , 2004;;Hosseinnia, Ranjbar, Ganji, Soltani, & Ghasemi, 2009, e.g.).For the considered simultaneous inversion problem here, the unknowns can not be determined if still using the optimal perturbation algorithm as done in the previous works.A possible approach for solving the inverse problem here is to combine the homotopy method with the optimal perturbation algorithm, which can be called the homotopy regularization algorithm.Recently in (Ruan, Xu, & Wang, 2011), the so-called homotopy regularization algorithm was ever utilized to solve the backward problem in the heat equation, however, we will still employ the inversion algorithm by the following two reasons: (i) The forward problem is more complicated than before.Here we will consider the time fractional diffusion equation, and the computational complexity is high; (ii) The inversion problem is different.We will simultaneously determine the space-dependent diffusion coefficient and the source magnitude with the final observations, which seems to be more difficult than those in the previous works.
This paper is organized as follows.Section 2 gives the implicit difference scheme for solving the forward problem, and in section 3 the homotopy regularization algorithm with the Sigmoid-type homotopy parameter is introduced to determine the space-dependent diffusion coefficient and the source magnitude simultaneously.In section 4, numerical inversions and simulations are performed by applying the homotopy regularization algorithm, and several concluding remarks are given in section 5.

The Difference Scheme to the Forward Problem
Consider Equation (1) with the initial and zero Neumann boundary conditions here u 0 (x) is the initial function.It is noted that an implicit finite difference scheme to Equation (1) was ever discussed in Li et al. (2012), here the discretization to the first-order derivative term in Equation ( 1) is a little different.We will utilize the forward difference method to discretize the first-order derivative instead of using the backward difference as before, and also we can get the unconditional stability and convergence of the difference scheme.
Firstly, based on Caputo's definition (2) for the fractional derivative, discretizing space domain by , and time domain by here h = 1/M is the space mesh step, and μ = 1/N is the time mesh step.
Next, discretizing the first-order derivative D (x) in Equation ( 1) by the forward difference instead of the backward difference as given in Li et al. (2012), and the second-order derivative also by the ordinary center difference, there are ∂u( and Let ), and thanks to (5), ( 6) and ( 7), we can get n k=0 where then Equation ( 8) can be re-arranged as follows: and the initial boundary conditions are discretized as Therefore, an implicit difference scheme can be given as follows.
For n = 0, there is For n > 0, there is Then, the implicit scheme ( 10)-( 11) can be rewritten in matrix form The above difference scheme is a little different from that work given in Li et al. (2012), but it is of the same stability and convergence as long as the diffusion coefficient D = D(x) being positive for x ∈ (0, 1).With a similar method as used in Li et al. (2012), it is easy to obtain the following result.
Theorem 1 The implicit difference scheme defined by ( 15) is unconditionally stable, and the difference solution is convergent to the exact solution with the order of O(h + μ) as h, μ → 0 for finite time domain.
In the following, we will employ the difference scheme (15) to solve the forward problem, and then to perform numerical inversions utilizing the homotopy regularization algorithm presented in the next section.

The Inverse Problem
In many cases for the anomalous diffusion model like Equation ( 1), the space-dependent diffusion coefficient, the source term and other physical quantities are often unknown and cannot be measured easily.However, these physical quantities can be identified from some additional information which can be measured easily.If the diffusion coefficient D(x) and the source magnitude g(x) are both unknown, we are to seek them by additional observations at the final time t = T , i.e., there is Therefore, an inverse problem is formulated by Equation (1), the initial boundary conditions (3)-( 4) with the overposed condition (16).
Suppose that the diffusion coefficient D(x) and the source magnitude g(x) are both continuous on the space domain [0, l], and D(x) > 0 for x ∈ (0, l), and Φ, Ψ ⊂ C[0, l] are the admissible sets for D(x) and g(x), respectively.A solution to the inverse problem can be represented by D(x), g(x) ∈ Φ × Ψ.Then for any prescribed D ∈ Φ and g ∈ Ψ, an unique solution of the corresponding forward problem, denoted by u(x, t; D, g), can be solved by the difference scheme (15), and then the so-called output u(x, T ; D, g) corresponding to the input D ∈ Φ and g ∈ Ψ is obtained.So, solving the inverse problem is transformed to minimize an error of the unknown between the output data and the additional data from the viewpoint of optimality, or it is equivalent to solve an error operator equation from the viewpoint of integral equation.However, both of the above direct approaches could be ill-posed, and often of nonlinearity and complexity which lead us to seek other effective and stable methods.

The Homotopy Method
For convenience of writing, denote the solution of the inverse problem as z = D, g ∈ Φ × Ψ, then solving the inverse problem is to find a solution of the operator equation: where u(x, T ; z) is the output corresponding to z ∈ Φ × Ψ, and θ(x) is the additional information given by ( 17).
Thus by the ordinary homotopy, let us define where λ ∈ [0, 1] is the homotopy parameter, and g(z) is an optional function whose solution can be obtained easily.
Then a solution of Equation ( 17) can be obtained by solving its homotopy Equation ( 18) when the homotopy parameter continuously varying from 0 to 1. On the other hand in reality, we are always to seek for a minimizer as a solution of the inverse problem, which making a norm of the error functional attain its minimum, i.e., to solve the minimization problem min where • * denotes a suitable norm on the space domain.Hence, also from the idea of homotopy, taking g(z) = z * , let us define the homotopy: here the homotopy parameter λ ∈ (0, 1) but taking values continuously decreasing from near 1 approaching to zero, which can be regarded as regularization parameter in the inversion algorithm given below.Noting to the optimal perturbation algorithm, combining with the above homotopy idea, we can deduce the homotopy regularization algorithm.

The Homotopy Regularization Algorithm
Suppose that {ϕ k (x), k = 1, 2, • • • , ∞} and {ψ s (x), s = 1, 2, • • • , ∞} are groups of basis functions of Φ and Ψ respectively, and there have approximate expansions given as follows and respectively, where D K (x), g S (x) are the K-dimensional and S -dimensional approximate solutions to D(x) and g(x) respectively, and K, S ≥ 1 are the truncated levels, and are the expansion coefficients of D(x) and g(x), respectively.It is convenient to set two limited dimensional spaces as and respectively, and , in which meaning we can say that (a, b) = D K , g S .In the following, we will set ā ) as an approximate solution to the inverse problem if there is no disorder.
Based on the above discussions, denoting u(x, t; ā) = u(x, t; D K , g S ) as the unique solution of the forward problem for any admissible diffusion coefficient D K (x) given by ( 20) and g S (x) given by ( 21), then combining with the homotopy method, a feasible way to solve the joint inversion problem numerically is to solve the following minimization problem where H(ā, λ) is the homotopy, and λ ∈ (0, 1) is the homotopy parameter, which plays the same role as regularization parameter in our inversion algorithm.
For any given ā where j ≥ 1 is the number of iterations, and δa j denotes a perturbation of ā j for each j.In the following for convenience of writing, ā j and δā j are abbreviated as ā and δā, respectively.
Taking Taylor's expansion for u(x, T ; ā + δā) at ā, and ignoring the high order terms, we get Noting to (24), let us define an error functional for the perturbation: Now, discretizing the space domain [0, l] with 0 = x 1 < x 2 < • • • < x N = l, and using the forward difference to approximate each term of the gradient vector ∇ T ā u(x, T ; ā), then we have where • 2 denotes the discrete Euclid norm, and here τ is the numerical differential step, and It is easy to verify that minimizing ( 27) can be reduced to solve the normal equation (see Kirsch, 1996, e.g.) In fact, dividing by a factor 1 − λ ∈ (0, 1) on the two side of the error functional (28), and denoting α = λ 1−λ , there is Obviously, the two error functionals ( 28) and ( 32) are equivalent, and to minimize (32) is equivalent to solve the normal equation Then, noting to α = λ 1−λ , we get the normal Equation (31).By the above discussions, we know that an optimal perturbation can be solved by the normal Equation ( 31 then an optimal solution can be approximated by the iteration procedure (25) as long as arriving at the given number of iterations, or the perturbation satisfying a prescribed convergent precision given as here eps is the given convergent precision.
The key points of implementing the above inversion algorithm lie in suitable choices of the approximate space Φ K for the diffusion coefficient, and Ψ S for the source magnitude, the homotopy parameter, the differential step, the initial iteration, and the convergent precision, etc.In the computations presented in the next section, we will choose polynomials space as the approximate space both for the diffusion coefficient and the source magnitude, and noting to positive and continuous decreasing properties of Sigmoid-type functions (see Haykin, 1994, e.g.), we will utilize a Sigmoid-type function depending on the number of iterations as the homotopy parameter given as: where j is the number of iterations, and j 0 is the preestimated number of iterations when the homotopy parameter decreasing to 0.5, and β > 0 is the adjust parameter.In all of the following numerical inversions, we will always choose j 0 = 5 and β = 0.5 in (36) if there is no specification.

Numerical Inversion
Without loss of generality, in the following computations we will take the attenuation factor in the source term of Equation ( 1) as ω(t) = exp(−t), and the final time as T = 1, and the initial function as u 0 (x) = x 2 (1 − x) 2 , and we will always set M = 100 and N = 100 in the difference scheme (15) for solving the forward problem, and take the differential step as τ = 1e − 2, the convergent precision as eps = 1e − 5, and the homotopy parameter given by ( 36), where j 0 = 5 and β = 0.5.In addition, all computations are performed in a PC of Dell Studio.
4.1 Ex. 1-Inversion for D(x) = 1 + x and g(x) = 1 − x Suppose that the exact diffusion coefficient is D(x) = 1 + x, and source magnitude is g(x) = 1 − x in this example, then by substituting them into the forward problem, the output u(x, T ; D, g) can be obtained with which the homotopy regularization algorithm is applied to reconstruct them in suitable approximate spaces.If performing the inversion algorithm in linear polynomials space Φ 2 = Ψ 2 = span{1, x}, then the exact solution of the inversion problem here is represented by ā = (1, 1, 1, −1).Firstly, taking the initial iteration as zero, the inversion results varying with the fractional orders are listed in Table 1, where γ denotes the fractional order, āinv denotes the inversion solution, and Err = ā − āinv 2 / ā 2 denotes relative error in the solutions, and j denotes the number of iterations.

Discussion and Conclusion
We give several conclusions and discussions in this section.
(i) The homotopy regularization algorithm with the homotopy parameter chosen by the Sigmoid transform function is suitable for the simultaneous inversion problem of determining the space-dependent diffusion coefficient and the source magnitude in the time fractional diffusion equation.The inversion solutions give good approximations to the exact solutions, and the inversion algorithm here is of numerical convergence with a large convergent range.
(ii) The order of the fractional derivative could indicate some ill-posedness of the related inverse problem.By the numerical computations, we find that the inversion always becomes complicated if taking larger fractional order greater than γ = 0.8, and it seems to be better when the fractional order belonging to the interval of [0.1, 0.8].Especially for Ex.3, the fractional order plays a key role in the implementation of the algorithm.
(iii) Generally speaking, the approximate space, the homotopy parameter, the differential step, the initial iteration, and computation of the forward problem are primary factors influencing the implementation of the inversion algorithm.However, by computations for the simultaneous inversion problem, we find that the initial iteration and the differential step seem to have little influences on the inversion.Theoretically speaking, the inversion coefficients should approximate to the exact solutions as the dimensions go to infinity.However, the situation here becomes complicated on concrete computations, and the inversions could become much difficult when performing the algorithm in high-dimensional approximate spaces.Maybe smooth inversions in high-dimensional approximate spaces can be obtained if choosing orthogonal basis functions instead of polynomials for the approximate space.Finally, by the inversion computations in the three examples, we find that the simultaneous inversion problem is sensitive to the choice of regularization parameter.The inversion algorithm can be realized adaptively and stably with the number of iterations increasing using the Sigmoid-type homotopy parameter given by (36), and the homotopy parameter could become very small.We will focus our attention on theoretical analysis and application of the homotopy regularization algorithm in the sequent works.
Figure 1.The exact and the inversion solutions in Ex.1

Figure 2 .
Figure 2. The inversion error with number of iterations in Ex.1

Figure 3 .
Figure 3.The exact and the inversion solutions in Ex.2 Figure 4.The exact and the inversion solutions in Ex.3

Table 3 .
The inversion results with number of iterations for γ = 0.5 in Ex.1

Table 5 .
Inversion results with initial iterations in Ex.2

Table 7 .
Inversion results with fractional orders in Ex.3