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.
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 in space is replaced by a fractional derivative of order , denoted as , we get the space FADE given as where denotes concentration of the solute in space point and time , denotes the order of the fractional derivative, is the dispersion coefficient, denotes the average velocity, and is a source term. The derivative here denotes a fractional derivative at space point from the left hand-side of the domain in the meaning of Grünwald-Letnikov or Riemann-Liouville (see , 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.  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.  considered an inverse problem of determining the space-dependent source function in (1.1) in the case of 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.  considered an inverse problem of simultaneously determining the diffusion coefficient and the source term in 1D integer-order 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.
2. Numerical Solution to the Forward Problem
2.1. The Implicit Difference Scheme
For and , 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 is assumed to be continuous on and .
Let the space step be and the time step be , with and for and . According to standard Grünwald formula (see , 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 , for and .
Denote and and , then this can be rearranged to yield for and . Let , , and is the coefficient matrix given by (2.6), where the coefficients are defined by for and , and for , respectively. Thus, we get the implicit finite difference scheme in matrix form given as follows:
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 for finite time domain.
2.2. 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 , , and take the initial function as , the source term as and utilize the definition of Grünwald-Letnikov fractional derivative; the -order fractional derivatives of the functions and are given as respectively, then the exact solution of the froward problem can be obtained which is given as .
In the numerical simulation, we will take as example, and set the diffusion coefficient , the velocity ; numerical results are listed in Tables 1 and 2, respectively, where denotes the number of space and time grids, respectively, and Err denotes relative error in the exact and numerical solutions at the given time. Moreover, the exact and numerical solutions at and at for are plotted in Figure 1, respectively.
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.
3. 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 , and the flow velocity 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 , , and 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 integer-order 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 the three parameters , , and simultaneously in the FADE. The inversion algorithm is given as follows.
Denote , and suppose that is the unique solution of the forward problem for any prescribed ; then a feasible way to solve the inverse problem here is to solve the minimization problem where is the regularization parameter, is the computational output obtained by the solution taking values at , and is the final observation referred to (3.1).
Suppose that here is called a perturbation for given . Thus, in order to get from the given , we only need to get an optimal perturbation . In what follows, for convenience of writing, and are abbreviated as and , respectively. Then, we only need to work out an perturbation vector
Taking Taylor’s expansion for at and ignoring higher order terms, we can get Noting to (3.2), and with the help of (3.5), we define an error functional given as follows: Now, discretizing the space domain with , where denotes the number of grids, the above norm can then be reduced to the discrete Euclidean norm given as where and is the numerical differentiation step, and , , and , and
It is not difficult to testify that minimizing functional (3.7) is equivalent to the solving of the following normal equation (see , e.g.) 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.
4. Numerical Inversion
Taking the space domain as , the final time as , the source term as , and the initial function as , we will carry out numerical simulations in this section. For the given exact parameters , , and in (1.1), the forward problem is solved by the difference scheme (2.9) with , , 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.
4.1. 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 , the diffusion coefficient is , and the average velocity is , that is, is regarded as the exact solution in this subsection. In the concrete implementation of the inversion algorithm, set initial iteration as , numerical differential step as , and convergent criterion as . The computational results are listed in Table 3, where denotes the regularization parameter, denotes the inversion solution, is the CPU time of each iteration, is the total CPU time for each inversion, denotes the number of iterations, and 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 , and the solutions error still remain even with very small regularization parameter . The inversion errors varying with regularization parameters in this case are plotted in Figure 2(a). Furthermore, if taking numerical differential step as , and other inversion parameters unchanged, the computational results are listed in Table 4, where , , Err, and 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.
4.2. 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 , and take the regularization parameter as , 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 , Err, and all denote the same meanings as referred in the above. Moreover, the errors in the solutions varying with the numerical differential steps for 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 in this example.
4.3. 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 . By the above computations, the regularization parameter here is , the differential step is , and the initial iteration is as before. The computational results are listed in Table 6, where denotes the exact solution, denotes the fractional order, and , Err, and 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 , which demonstrates that ill-posedness of the inversion problem could become severe if coping with small fractional orders.
4.4. 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 is the accurate additional information given by (3.1), is noise level, and is a random vector ranged in . In the following, we will take the fractional order as example, that is, the exact solution in this subsection is given as .
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 and and we choose regularization parameter as , and , 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, denotes the average inversion solution of the ten-time computations, denotes the average relative error in the solutions, and with 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 with different noise levels, where also denotes the noise level and , , and with 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.
(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 , the solutions error becomes small and the computational complexity is reduced; however, if is less than , 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 in our sequent works.
The authors thank for the referees and the editor for their fruitful suggestions and comments. This work is supported by the National Natural Science Foundation of China (no. 11071148) and the Natural Science Foundation of Shandong Province, China (no. ZR2011AQ014).
E. E. Adams and L. W. Gelhar, “Field study of dispersion in a heterogeneous aquifer: 2,” Spatial Moments Analysis, Water Resources Research, vol. 28, no. 12, pp. 3293–3307, 1992.View at: Google Scholar
D. A. Benson, The fractional advection-dispersion equation: development and application [Ph.D. thesis], University of Nevada, Reno, Nev, USA, 1998.
Y. Hatano and N. Hatano, “Dispersive transport of ions in column experiments: an explanation of long-tailed profiles,” Water Resources Research, vol. 34, no. 5, pp. 1027–1033, 1998.View at: Google Scholar
Y. Pachepsky, D. Benson, and W. Rawls, “Simulating scale-dependent solute transport in soils with the fractional advective-dispersive equation,” Soil Science Society of American Journal, vol. 64, no. 4, pp. 1234–1243, 2000.View at: Google Scholar
B. Berkowitz, H. Scher, and S. E. Silliman, “Anomalous transport in laboratory-scale heterogeneous porous media,” Water Resources Research, vol. 36, no. 1, pp. 149–158, 2000.View at: Google Scholar
L. Z. Zhou and H. M. Selim, “Application of the fractional advection-dispersion equations in porous media,” Soil Science Society of American Journal, vol. 67, no. 4, pp. 1079–1084, 2003.View at: Google Scholar
Y. W. Xiong, G. H. Huang, and Q. Z. Huang, “Modeling solute transport in one-dimensional homogeneous and heterogeneous soil columns with continuous time random walk,” Journal of Contaminant Hydrology, vol. 86, no. 3-4, pp. 163–175, 2006.View at: Google Scholar
I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, Calif, USA, 1999.
F. Liu, V. Anh, and I. Turner, “Numerical solution of the space fractional Fokker-Planck equation,” Journal of Computational and Applied Mathematics, vol. 166, no. 1, pp. 209–219, 2004.View at: Google Scholar
Q. Z. Huang, G. H. Huang, and H. B. Zhan, “A finite element solution for the fractional advection dispersion equation,” Advances in Water Resources, vol. 31, no. 12, pp. 1578–1589, 2008.View at: Google Scholar
B. T. Jin and W. Rundell, “An inverse problem for a one-dimensional time-fractional diffusion problem,” Inverse Problems, vol. 28, no. 7, Article ID 075010, 2012.View at: Google Scholar
G. S. Li, W. J. Gu, and X. Z. Jia, “Numerical inversions for space-dependent diffusion coefficient in the time fractional diffusion equation,” Journal of Inverse and Ill-Posed Problems, vol. 20, no. 3, pp. 339–366, 2012.View at: Google Scholar
F. A. Rodrigues, H. R. B. Orlande, and G. S. Dulikravich, “Simultaneous estimation of spatially dependent diffusion coefficient and source term in a nonlinear 1D diffusion problem,” Mathematics and Computers in Simulation, vol. 66, no. 4-5, pp. 409–424, 2004.View at: Publisher Site | Google Scholar | Zentralblatt MATH
G. S. Li, J. Cheng, D. Yao, H. L. Liu, and J. J. Liu, “One-dimensional equilibrium model and source parameter determination for soil-column experiment,” Applied Mathematics and Computation, vol. 190, pp. 1365–1374, 2007.View at: Google Scholar
G. S. Li, Y. J. Tan, D. Yao, X. Q. Wang, and H. L. Liu, “A nonlinear mathematical model for undisturbed soil-column experiment and source parameter identification,” Inverse Problems in Science and Engineering, vol. 16, no. 7, pp. 885–901, 2008.View at: Google Scholar
G. S. Li, D. Yao, and Y. Z. Wang, “Data reconstruction for a disturbed soil-column experimentusing an optimal perturbation regularization algorithm,” Journal of Applied Mathematics, vol. 2012, Article ID Article ID 732791, 16 pages, 2012.View at: Google Scholar
C. W. Su, Numerical Methods and Applications of Inverse Problems in PDE, Northwestern Polytechnical University Press, Xi’an, China, 1995, (in Chinese).
A. Kirsch, An Introduction to the Mathematical Theory of Inverse Problems, Springer, New York, NY, USA, 1996.View at: Publisher Site