Numerical Identification of Multiparameters in the Space Fractional Advection Dispersion Equation by Final Observations

This paper deals with an inverse problem for identifying multiparameters in 1D space fractional advection dispersion equation FADE on a finite domain with final observations. The parameters to be identified are the fractional order, the diffusion coefficient, and the average velocity in the FADE. The forward problem is solved by a finite difference scheme, and then an optimal perturbation regularization algorithm is introduced to determine the three parameters simultaneously. Numerical inversions are performed both with the accurate data and noisy data, and several factors having influences on realization of the algorithm are discussed. The inversion solutions are in good approximations to the exact solutions demonstrating the efficiency of the proposed algorithm.


Introduction
Many diffusion processes in nature and engineering, such as contaminants transport in the soil, oil flow in porous media, long distance transport of pollutants in the groundwater, and so forth, are referred to as anomalous diffusion, where the particle plume spreads faster or slower than predicted by the classical diffusion model.In recent two decades, the space FADE has been found to be an efficient mathematical model to describe such anomalous diffusion phenomena see 1-8 , e.g. .If the usual second-order derivative ∂ 2 /∂x 2 in space is replaced by a fractional derivative of order 1 < α < 2, denoted as ∂ α /∂x α , we get the space FADE given as where c c x, t denotes concentration of the solute in space point x and time t, 1 < α < 2 denotes the order of the fractional derivative, D > 0 is the dispersion coefficient, v > 0 denotes the average velocity, and S x, t is a source term.The derivative ∂ α c/∂x α here denotes a fractional derivative at space point x from the left hand-side of the domain in the meaning of Gr ünwald-Letnikov or Riemann-Liouville see 9 , e.g. .
There are quite a few research works on the fractional differential equations, refer to 10-16 for numerical methods and model simulations on 1.1 , and there are still a few studies in the time fractional diffusion equations recently see 17-22 , e.g. .However, research on inverse problems for the space FADE has not been paid much attention to our knowledge.Wei et al. 23 studied an inverse source problem in the spatial fractional diffusion equation and solved the inverse problem numerically based on the best perturbation method with Tikhonov regularization.Chi et al. 24 considered an inverse problem of determining the space-dependent source function f x in 1.1 in the case of S x, t e −t f x with final observations and presented numerical inversions by applying the optimal perturbation regularization algorithm.
In this paper, we will deal with an inverse problem of simultaneously determining the fractional order, the diffusion coefficient, and the average velocity in 1.1 on a finite domain also with the final observations.Rodrigues et al. 25 considered an inverse problem of simultaneously determining the diffusion coefficient and the source term in 1D integerorder diffusion equation by the conjugate gradient method with aids of an adjoint problem, and there seem to have few similar works reported in the known literatures especially for the space fractional diffusion equations.It is noticeable that the inversion problem here for simultaneously identifying these three kinds of model parameters becomes difficult as compared with those problems of determining a single unknown, and regularization strategies have to be utilized to overcome the ill-posedness.The forward problem is solved numerically by utilizing an implicit difference scheme based on the shifted Gr ünwall formula which was discussed in 10, 12, 24 .In regard of simultaneous inversion for the three kinds of parameters, we will apply the optimal perturbation regularization algorithm which was developed in 26-29 .However, since the inverse problem is different, the inversion algorithm has to be changed in detail so as to get effective solutions.The inversion algorithm is performed not only in the case of using the accurate data, but also with the random noisy data.The inversion solutions give good approximations to the exact solutions showing that the inversion algorithm is stable and suitable for the inverse problem of determining multiparameters in the FADE.The paper is arranged as follows.
In Section 2, an implicit finite difference scheme for solving the forward problem is given, and a numerical example is presented to support the numerical method.Section 3 gives the inverse problem of determining the parameters of the fractional order, the diffusion coefficient, and the flow velocity simultaneously by final observations, and the optimal perturbation regularization algorithm is introduced.In Section 4, numerical inversions both with the accurate data and noisy data are implemented, and the factors influencing realization of the inversion algorithm are discussed.Finally, several concluding remarks are given in Section 5.

The Implicit Difference Scheme
For 0 < x < l and 0 < t < T, consider an initial boundary value problem for the FADE given by 1.1 , where the initial boundary value conditions are given as where the function g x is assumed to be continuous on 0, l and g 0 g l 0. Let the space step be h l/M and the time step be r T/N, with x i ih and t n nr for i 0, 1, . . ., M and n 0, 1, . . ., N. According to standard Gr ünwald formula see 9 , e.g., , a discrete approximation to the fractional derivative is given as follows: with which a modified approximation which is called the shifted Gr ünwald formula is obtained with the help of Gamma function see 10, 12 , e.g.
Then we have a difference scheme by 1.1 with the initial boundary value conditions 2.1 where c n i ≈ c x i , t n , S n i ≈ S x i , t n for i 0, 1, . . ., M and n 0, 1, . . ., N − 1. Denote and p vr/h and q Dr/h α , then this can be rearranged to yield for i 1, 2, . . ., M − 2 and j 1, 2, . . ., M − 1, and for i M − 1, respectively.Thus, we get the implicit finite difference scheme in matrix form given as follows: Ac n 1 c n rS n 1 , n 0, 1, . . ., N − 1.

2.9
For the above difference scheme, similarly as done in 10, 12 , it is not difficult to get its convergence and stability with the help of spectrum analysis to the coefficient matrix.Lemma 2.1 see 10, 12 .The implicit difference scheme defined by 2.9 is unconditionally stable, and the difference solution is convergent to the exact solution of the forward problem as h, r → 0 for finite time domain.

Numerical Testification
In this subsection, an example is presented to show effectiveness of the finite difference scheme 2.9 for solving the forward problem numerically.For the forward problem 1.1 -2.1 , set l π, T 1, and take the initial function as g x x π − x , the source term as and utilize the definition of Gr ünwald-Letnikov fractional derivative; the α-order fractional derivatives of the functions x and x 2 are given as respectively, then the exact solution of the froward problem can be obtained which is given as c x, t x π − x e −t .In the numerical simulation, we will take α 1.6 as example, and set the diffusion coefficient D 0.1, the velocity v 1; numerical results are listed in Tables 1 and 2, respectively, where M N denotes the number of space and time grids, respectively, and E rr denotes relative error in the exact and numerical solutions at the given time.Moreover, the exact and numerical solutions at t 0.4 and at t 0.8 for M N 40 are plotted in Figure 1, respectively.

2.11
From Tables 1 and 2 and Figure 1, we can see that the difference scheme given by 2.9 is of numerical convergence; and the numerical solutions are in good agreement with the exact solution.In what follows, numerical inversions can be performed based on the the above numerical method given by 2.9 by applying the optimal perturbation regularization algorithm.

The Inverse Problem and the Inversion Algorithm
In many cases for the anomalous diffusion model given by 1.1 , the fractional order, the diffusion coefficient, and other physical quantities are often unknown and cannot be measured easily.However, these physical quantities can be identified from some known data which can be observed or measured easily based on the forward problem.The method of obtaining these quantities is considered as an inverse problem.Suppose that the fractional order α, the diffusion coefficient D, and the flow velocity v in 1.1 are all unknown, and we have final observation as the additional information given as Thus, an inverse problem of simultaneously determining these three parameters α, D, and v is formulated by 1.1 and the initial boundary value conditions 2.1 with the additional condition 3.1 .In what follows, we are to deal with the above multiparameters identification problem from numerics by utilizing an optimal perturbation algorithm.As we know, most of inversion algorithms are based on regularization strategies so as to overcome ill-posedness of the inverse problem, and different kinds of inverse problems could need different approximate methods on the basis of conditional well-posedness analysis see 25, 30-33 , e.g. .It is noted that the optimal perturbation algorithm has been attested to be effective for identifying model parameters at least for the ordinary integerorder diffusion equations see 26-29, 34 , e.g. , and it is also efficient for determining source functions in the FADE see 23, 24 , e.g. .In this section, we will also employ it to determine Denote a α, D, v a 1 , a 2 , a 3 , and suppose that c x, t; a is the unique solution of the forward problem for any prescribed a ∈ R 3 ; then a feasible way to solve the inverse problem here is to solve the minimization problem min where μ > 0 is the regularization parameter, c x, T ; a is the computational output obtained by the solution taking values at t T , and θ x is the final observation referred to 3.1 .Suppose that a j 1 a j δa j , j 0, 1, . . ., 3.3 here δa j is called a perturbation for given a j .Thus, in order to get a j 1 from the given a j , we only need to get an optimal perturbation δa j .In what follows, for convenience of writing, a j and δa j are abbreviated as a and δa, respectively.Then, we only need to work out an perturbation vector δa δa 1 , δa 2 , δa 3 .

3.4
Taking Taylor's expansion for c x, T ; a δa at a and ignoring higher order terms, we can get c x, T ; a δa ≈ c x, T ; a ∇ T c x, T ; a • δa.

3.5
Journal of Applied Mathematics 7 Noting to 3.2 , and with the help of 3.5 , we define an error functional given as follows: μ δa 2 2 .

3.6
Now, discretizing the space domain 0, l with 0 , where K denotes the number of grids, the above L 2 norm can then be reduced to the discrete Euclidean norm given as where and τ is the numerical differentiation step, and e 1 1, 0, 0 , e 2 0, 1, 0 , and e 3 0, 0, 1 , and 3.9 It is not difficult to testify that minimizing functional 3.7 is equivalent to the solving of the following normal equation see 35 , e.g.μI B T B δa B T η − ξ .

3.10
So, an optimal perturbation can be solved via the formula: with which an optimal inversion solution can be approximated by the iteration procedure 3.3 as long as the perturbation is satisfying a given precision.Obviously, it lies in suitable choices of the regularization parameter, the numerical differential step, the initial iteration, and the convergent precision to perform the above inversion algorithm.In the next section, numerical inversions will be performed by employing the above inversion algorithm.

Numerical Inversion
Taking the space domain as 0, π , the final time as T 1, the source term as S x, t xe −t , and the initial function as g x x 2 π − x , we will carry out numerical simulations in this section.For the given exact parameters α, D, and v in 1.1 , the forward problem is solved by the difference scheme 2.9 with M 100, N 20, and then the additional data are obtained with which the inversion algorithm is performed to reconstruct the exact parameters.All computations are finished in a PC of Tsinghua Tongfang.

Influence of Regularization Parameter on the Algorithm
It is important to choose an optimal regularization parameter by theoretical analysis in solving inverse problems, where regularization methods are needed.However, it is still a feasible approach to the choice of regularization parameter by trial and error especially for moderate ill-posed inverse problems.The regularization parameter here is selected by numerical testification, and its influence on the inversion algorithm for given differential step is discussed in this subsection.Suppose that the fractional order is α 1.6, the diffusion coefficient is D 1, and the average velocity is v 1, that is, a 1.6, 1, 1 is regarded as the exact solution in this subsection.In the concrete implementation of the inversion algorithm, set initial iteration as a 0  1.1, 0, 0 , numerical differential step as τ 0.1, and convergent criterion as δa 2 ≤ 1e−4.The computational results are listed in Table 3, where μ denotes the regularization parameter, a inv denotes the inversion solution, T cpu /I is the CPU time of each iteration, T cpu is the total CPU time for each inversion, j denotes the number of iterations, and E rr a inv − a 2 / a 2 is the relative error in the solutions.
It is noticeable that the inversion results have little changes if utilizing any positive regularization parameters smaller than 1e − 10, and the solutions error still remain E rr 3.0466225e −5 even with very small regularization parameter μ 1e −36.The inversion errors varying with regularization parameters in this case are plotted in Figure 2 a .Furthermore, if taking numerical differential step as τ 0.01, and other inversion parameters unchanged, the computational results are listed in Table 4, where μ, a inv , E rr , and T cpu /j all denote the same meanings as used in Table 3.
From Tables 3 and 4 and Figure 2 a , we can find that regularization parameter has important impact on the inversion algorithm for given numerical differential step, and the regularization parameter must be chosen large for small numerical differential step.However, the inversion algorithm can be performed smoothly with any small regularization parameters strictly greater than zero if choosing large numerical differential steps.

Influence of Numerical Differential Step on the Algorithm
By the above computations together with 3.8 , we find that numerical differentiation is another key point in the realization of the inversion algorithm, and suitable differential steps are necessary in order to realize a successful inversion.As done for the regularization parameter, we will investigate the influence of the differential step on the algorithm by testification.Also set the exact solution to be a 1.6, 1, 1 , and take the regularization parameter as μ 0.1, and other inversion parameters also unchanged as used in the above.The inversion results are listed in Table 5, where τ denotes the numerical differential step, and a inv , E rr , and T cpu /j all denote the same meanings as referred in the above.Moreover, the errors in the solutions varying with the numerical differential steps for μ 0.1 are plotted in Figure 2 b .
From Table 5 and Figure 2 b , we can see that the numerical differential step plays a similar role as to that of the regularization parameter in the realization of the algorithm.The inversion algorithm could be a failure in case of using large or small differential steps.The differential steps should fall in the interval of 0.03, 0.9 in this example.

Influence of Fractional Order on the Algorithm
In this subsection, we will investigate influence of the fractional order on the inversion algorithm with the convergent criterion given as δa 2 ≤ 10 −6 .By the above computations, the regularization parameter here is μ 0.1, the differential step is τ 0.1, and the initial iteration is a 0 1.1, 0, 0 as before.The computational results are listed in Table 6, where a denotes the exact solution, α denotes the fractional order, and a inv , E rr , and T cpu /j also refer to the same notations as used in the above.
From Table 6, we find that the inversion results are not so good when the fractional order is smaller than 1.2, which demonstrates that ill-posedness of the inversion problem could become severe if coping with small fractional orders.

Inversion with Noisy Data
It is difficult to perform an inversion algorithm in the case of using random noisy data, especially for inverse problems arising from the fractional diffusion.Noting computational errors and data noises, the additional information utilized for real inverse problems is often given as where θ x is the accurate additional information given by 3.1 , ε > 0 is noise level, and ζ is a random vector ranged in −1, 1 .In the following, we will take the fractional order α 1.5 as example, that is, the exact solution in this subsection is given as a 1.5, 1, 1 .Generally speaking, in case of using regularization strategy to damp the noises, an optimal regularization parameter should be chosen according to the noise level.However, the situation becomes simple for the inverse problem considered here.Without loss of generality, we set the noise levels be ε 1% and ε 5% and we choose regularization parameter as μ 0.5, μ 0.1 and μ 0.01, respectively, and other inversion parameters also unchanged as in the above.The average inversion results by ten-time computations are listed in Tables 7  and 8 respectively, where ε denotes the noise level, a inv denotes the average inversion solution of the ten-time computations, E rr a inv − a 2 / a 2 denotes the average relative error in the solutions, and T cpu with j denote the average CPU time and iteration number of the inversion, respectively.
From Tables 7 and 8, we can see clearly that for the given noise level, the inversion solutions are almost unchanged, although the number of iterations and the corresponding CPU time increase to some extent when using large regularization parameters.Furthermore, Table 9 gives inversion results by using regularization parameter as μ 0.1 with different noise levels, where ε also denotes the noise level and a inv , E rr , and T cpu with j denote the same meanings as used in the above.
From Table 9, it can be seen that the inversion results are satisfactory in the case of using noisy data, and the error in the solutions becomes small and approaches to zero with the noise level goes to zero demonstrating that the inversion algorithm is of numerical stability.

Conclusion
i An inverse problem for identifying multiparameters of the fractional order, the diffusion coefficient, and the average velocity simultaneously in the FADE with final observations is investigated.An implicit finite difference scheme is employed to solve the forward problem, and an optimal perturbation regularization algorithm is applied to reconstruct the three parameters.By the numerical simulations, we conclude that the optimal perturbation regularization algorithm is of numerical stability, and it is efficient for multiparameters identification problem arising from the FADE.
ii There are few factors that have influences on the inversion algorithm, but the regularization parameter and the numerical differential step here seem to play a much more important role in the realization of the algorithm.The inversion algorithm can be performed successfully for the regularization parameter and the numerical differential step belonging to suitable intervals, respectively.However, the inversion algorithm may be a failure if the numerical differential step is taking too small values.In addition, it seems to be insensitive to the choice of the regularization parameter in concrete realization of the algorithm showing that the inverse problem discussed here is of moderate ill-posedness.
iii By the inversion computations, we find that the fractional order in the FADE has some influence on the forward problem and the inverse problem.If the fractional order α goes to 2, the solutions error becomes small and the computational complexity is reduced; however, if α is less than 1.2, the ill-posedness of the inverse problem becomes severe, and the inversion results could be unacceptable.We will deal with the forward problem and the corresponding inverse problem with the fractional order smaller than 1.2 in our sequent works.

Figure 1 :
Figure 1: The exact and numerical solutions with M N 40.

Figure 2 :
Figure 2: The solutions errors with regularization parameter and differential step.

Table 1 :
Errors in the solutions at t 0.4 with number of grids.

Table 2 :
Errors in the solutions at t 0.8 with number of grids.

Table 3 :
Influence of regularization parameter on the algorithm with τ 0.1.

Table 4 :
Influence of regularization parameter on the algorithm with τ 0.01.

Table 5 :
Influence of numerical differential step on the algorithm with μ 0.1.

Table 6 :
Influence of fractional order on the inversion algorithm.