Inverse Problem for a Two-Dimensional Anomalous Diffusion Equation with a Fractional Derivative of the Riemann–Liouville Type

The article presents a method for solving the inverse problem of a two-dimensional anomalous diffusion equation with a Riemann–Liouville fractional-order derivative. In the first part of the present study, the authors present a numerical solution of the direct problem. For this purpose, a differential scheme was developed based on the alternating direction implicit method. The presented method was accompanied by examples illustrating its accuracy. The second part of the study concerned the inverse problem of recreating the model parameters, including the orders of the fractional derivative, in the anomalous diffusion equation. Equations of this type can be used to describe, inter alia, the heat conductivity in porous materials. The ant colony optimization algorithm was used to solve this problem. The authors investigated the impact of the distribution of measurement points, the use of different mesh sizes, and the input data errors on the obtained results.


Introduction
In recent years, derivatives of fractional order have been frequently used. Many different definitions of such derivatives are available in the literature. In applications, the Caputo derivative and the Riemann-Liouville derivative are the most often used [1,2]. Fractional calculus is used, inter alia, in mathematical modeling of various problems in mechanics, physics, and biology [3][4][5][6], in control theory [7], in image processing [8], as well as in modeling the phenomena of heat transport [9][10][11][12][13][14][15]. Many mathematical models are based on differential equations containing fractional derivatives. Hence, the development of algorithms for solving differential equations with fractional derivatives is extremely important [16][17][18][19][20]. For example, in [16], the authors presented a solution of fractional pantograph differential equations. They used the operational matrix of the derivatives of the generalized Lucas polynomials and then transformed the problem into a system of equations. The paper also included the derivation of the numerical scheme with examples. Direct problems for differential equations consist of solving a given equation (determining the state function) when all input data are known. Solving inverse problems, on the other hand, involves determining some of the input data on the basis of some information regarding the state function. In the literature, studies on inverse problems for differential equations with fractional derivatives are available. The first studies on these issues were by Murio [21,22]. In these papers, the one-dimensional problem was considered, and the mollification method was used to solve it. Zheng and Wei [23,24] presented a Fourier regularization and a convolution methods to recreate the boundary condition (function and its derivative) of a semi-infinite domain for the time-fractional inverse advection-dispersion equation. A similar issue is considered in the finite domain in [25], in this case a kernelbased meshless method was used. Xiong et al. [26] considered a two-dimensional case in the infinite domain. Yan and Yang [27] used an efficient Kansa-type method and presented numerical solution of two-dimensional time-fractional inverse diffusion problems. For the theoretical results of the stability and the uniqueness of the solution, a backward problem and inverse source problem for a multi-dimensional time-fractional diffusion equation were included in [28]. The authors of [29,30] managed the inverse problem for the two-dimensional time-fractional sideways heat equation in the infinite domain. Inverse problems for two-dimensional time-fractional diffusion equation were also considered in [29,[31][32][33]. The study [34] considered the inverse problem for two-dimensional fractional partial differential equations. The method assumes the knowledge of the state function for a certain time interval at each point of the domain under consideration. A hybrid method based on modulating functions was used for the solution. The inverse problem for the heat equation in the two-dimensional domain with the Atangana-Baleanu fractional derivative was considered in [35]. In turn, the equation with the Hilfer derivative was considered in [36].
Research on one-dimensional issues prevails in the literature. Those on two-dimensional problems mainly regard the Caputo derivative in an infinite domain. There are a few articles on the inverse problem for the two-dimensional anomalous diffusion equation with the Riemann-Liouville derivative. Furthermore, the artificial intelligence algorithm was used in only one study in relation to time-fractional Schrödinger equations. The authors of this manuscript considered this rarely described case.
One of the main problems with using models with the fractional derivative is determining the order of the derivative. There are papers in which the authors determined it based on physical experiments [37]. To determine the order of the derivative, it may be useful to solve the inverse problem, which the authors present in this paper.
The article presents the consideration of the inverse problem for a two-dimensional anomalous diffusion equation with the Riemann-Liouville fractional derivative. The inverse problem in the study was to recreate the orders of two fractional derivatives and one of the parameters in the equation. Additional information, in this case, was the solution values at selected points of the domain. To solve the inverse problem, the direct problem should be repeatedly solved for the fixed values of the searched parameters. To solve the direct problem, the alternating direction implicit method was used [38], which allows for solving of multidimensional problem by the iterative solution of a number of onedimensional problems. Using the given values of the solution and the values determined from the solution of the direct problem, a functional representing the error of the approximate solution was built. The argument for which the functional reached the minimum value was the sought solution for the inverse problem under consideration. The minimum of the functional was searched for using the ant colony optimization algorithm [39,40]. This is a probabilistic artificial intelligence algorithm inspired by the behavior of an ant colony. The authors of the study investigated the impact of the distribution of measurement points, the use of different mesh sizes, and the value of input data disturbances on the obtained results.

Space Fractional Diffusion Equation
The presented model in this section can be used to describe heat flow in a porous material with an irregular (fractal, chaotic) structure. In such a case, the heat flow can be modeled using fractional-order derivatives [10]. We considered a two-dimensional problem in which the derivative over space was a Riemann-Liouville type. Thus, the study considered the two-dimensional space fractional diffusion equation with the Riemann-Liouville derivative: specified in the following area (x, y, t) ∈ Ω × [0, T], where Ω = [0, L x ] × [0, L y ] and α, β ∈ (0, 1). Initial-boundary conditions were attached to the above equation: It was also assumed that both constants c, and functions λ x1 , λ x2 , λ y1 , λ y2 > 0 are positive for all (x, y) ∈ Ω. Equation (1) describes the phenomenon of anomalous diffusion, for example the phenomenon of heat conduction in porous materials [9][10][11]. In Equation (1), both the right-and left-hand Riemann-Liouville derivatives were used as fractional derivatives relative to space [1,2]: where α ∈ (0, 1). Similar to the formula (3), we can define the fractional derivatives for a variable y and order β ∈ (0, 1).
Considering Equation (1) with the conditions (2) as a model describing the heat conductivity in porous materials, the following notations can be assumed: c-specific heat of a medium, • -density of a medium, • u-a function describing the temperature distribution in space and time, • f -an additional heat source.

Direct Problem
In this section, we present a way of solving Equation (1) numerically, and then, we introduce an example illustrating the accuracy of the method. For the discretization of the equation, the finite difference method supplemented with the appropriate discretization of Riemann-Liouville was used.
We approximate the Riemann-Liouville derivative as follows [41] (the order of approximation O((∆x) 2 )): where: The formulas for the fractional derivative with respect to the variable y look similar. The derivative with respect to time is approximated as follows (the order O((∆t) 2 )): Taking into account the approximations (5)-(7), we can write: For simplicity, we introduce the following notation: The approximation formulas for the derivative of the variable y are analogous. Then, the difference scheme of the order O((∆t) 2 + (∆x) 2 + (∆y) 2 ) for Equation (4) takes the form: After transforming the above equation, we can write: We write the above difference scheme in matrix form. For this purpose, we define two matrices: The next step is to define block matrices S and H. We start with the matrix S (di- ), which is a diagonal block matrix, i.e., it contains matrices R x (l), l = 1, 2, . . . , M y − 1 on the main diagonal and zeros in other places.
Below is the form of matrix H (matrix H has the same dimension as matrix S). Matrix H is a block matrix where each block is a diagonal matrix of dimension We can write the difference scheme (15) in matrix form as follows: where: The matrices appearing in the difference scheme (20) are quite large; as a consequence, the resulting system of equations will be time consuming to solve, which is shown in Section 3.2. Hence, we applied the alternating direction implicit method (ADIM) to the difference scheme (15), which significantly reduces the computation time. This will be an important issue in the case of solving the inverse problem, in which direct problem should be solved repeatedly. We write the scheme (15) in the form of the directional separation product: Then, the above scheme can be broken in two parts and solved, respectively, first in direction x and, afterwards, in direction y. As a consequence, the resulting matrices for the system of equations will have significantly lower dimensions than in the case of the scheme (15). The solution algorithm consists of two successive steps: • for each fixed y j , we solve the scheme in the direction x. As a result, we obtain a temporary solution u k+1 i,j : • then, for each fixed x i , we solve the scheme in direction y: We can represent this process symbolically by means of Figure 1. For the boundary nodes and the initial condition, we assumed: Furthermore, for the ADIM method, we write the appropriate matrix equations. First, for each l = 1, 2, . . . , M x − 1, we define auxiliary vectors u * l : where . Then, the scheme (22) can be written in the following matrix form (for p = 1, 2, . . . , M y − 1): where the temporary solution takes the form u k where u k+1 l = [u k+1 l,1 , u k+1 l,2 , . . . , u k+1 l,M y −1 ] T and u * k l = [ u k l,1 , u k l,2 , . . . , u k l,M y −1 ] T . From this step of the algorithm, we solved M x − 1 systems of equations with dimensions (M y − 1) × (M y − 1) each. The Bi-CGSTAB [42,43] method was used to solve the equation systems, which also influenced the computation time. A comparison with the Gauss algorithm is presented in the next section. The proof of the convergence of the above difference schemes can be found in [38].

Numerical Results
In this section, we present the results of the numerical solution of the problem (1) and (2). The example assumes the following data: For the data defined above, the exact solution (1) and (2) is the function: where the domain under consideration is as follows: Errors are marked as follows: where i = 0, 1, . . . , M x , j = 0, 1, . . . , M y , k = 0, 1, . . . , N, and ue k i,j denotes exact values of the function u at the mesh points. First of all, we compared the results obtained with and without the ADIM for different meshes. We also provided computation times. As can be seen in Table 1, the results obtained with or without ADIM were similar; the differences were slight. However, the difference in the computation time was of considerable size; for example, for the mesh 100 × 100 × 100, in the case of the ADIM, this time was about 9 s, while for the method without ADIM, this time was about 453 s. The difference was enormous, and considering that in the inverse problem, it is necessary to solve the direct problem repeatedly, the use of ADIM was essential. Another important issue influencing the computation time was the appropriate selection of the method to solve the system of equations. In the presented algorithm, the Bi-CGSTAB method [42,43] was used to solve a system of linear equations. We compared the computation times of the algorithm in the case of using Bi-CGSTAB and the Gaussian elimination method (Table 2). In the case of mesh 200 × 200 × 200, the use of the Bi-CGSTAB caused the computation to be performed approximately 5.3 times faster than using the Gaussian elimination method.

Inverse Problem
In this section, we present a solution of the inverse problem involving the reconstruction of the selected parameters of the model described by Equation (1). The considered inverse problem consisted of selecting the unknown input parameters of the model in such a way that the model output (function u from Equation (1)) had predetermined values at selected points of the domain. A lack of information regarding the searched input parameters was compensated by the values of the function u at selected points of the domain. Due to the fact that the inverse problems are ill-conditioned, they are classified as quite difficult issues.

Formulation of the Problem
We considered the following equation: where: λ = 240, = 2100, ϕ(x, y) = u(x, y, 0) = 0, and function f (x, y, t) is given by the formula (27). We considered the differential equation where N 1 is the number of measurement points, and N 2 is the number of measurements. The values u k i,j (c, α, β) are obtained from the solution of the direct problem for fixed parameter values c, α, β.

Objective Function Minimization
In order to find minimum of the objective function (28), and then solve inverse problem, the Ant Colony Optimization (ACO) algorithm was used [39,40]. The ACO algorithm is inspired by behavior of an ant colony in nature. In the considered version of the algorithm, solutions, which in our case were vectors, are treated as the so-called pheromone spots. In the initial phase of the algorithm, these spots are randomly distributed in the considered domain. Then, they are sorted by quality, and each solution (pheromone spot) is assigned a probability of choice. The better the solution is, the higher the probability is. The probability is determined on the basis of the quality of the solution and the q parameter of the algorithm, which we can control. In the case of a low q, the best solution is preferred, and as q increases, the probability of choosing any solution becomes similar. As a result, we obtained an ordered set of solutions. Then, the iterative phase of the algorithm begins, in which each ant (which are M) selects a spot according to its probability and then transforms it to sample the neighborhood. In the considered algorithm, to transform the solution, we used the Gaussian distribution density function. As a result of the transformation, we obtained M new solutions (pheromone spots). Next, the solution set must be updated with new pheromone spots (solutions) by sorting all solutions according to quality, and then, M the worst solutions are deleted. We repeated the iterative process I times. Finally, we chose the best solution. The algorithm was adapted to parallel computations. A more detailed description of the algorithm can be found in [9].

Results
In the considered inverse problem, we determined the parameters α, β, c. The exact values of these parameters were respectively 0.8, 0.6, 900. As we mentioned earlier, to test the stability of the algorithm, we assumed different arrangements of the measurement points, mesh sizes, and input data disturbances. For example, in Figure 3, we can see a comparison of the input data not disturbed by the error and disturbed by the 10% error from the measurement point K4. The following parameter values were adopted in the ant algorithm: Due to the probabilistic nature of the ant algorithm, we assumed five runs for the given settings, and then, the best one was selected from the obtained results. For this reason, in the tables with the results, the value of the standard deviation for the obtained values of the objective function was added. As can be seen, they were inconsiderable, which proved that the algorithm gave very similar results each time. Table 3 shows the results of reconstructing the searched parameters for various meshes and the values of input data disturbances in the case of twelve measurement points. In each case, the errors of the results obtained were low and did not exceed 1%. For the input data without error, the results were the best. Even with 10% disturbance of the input data, the obtained results were satisfying, which proved the stability of the model and algorithm. Table 4 contains analogous data for the case of seven measuring points. As can be seen, the change in the arrangement of measurement points, as well as their number influenced the obtained results. In a few cases (for example, for the 100 × 100 × 200 mesh and 10% disturbance), the errors in recreating selected parameters exceeded 1%, which did not happen in the case of twelve measuring points. It can therefore be concluded that the arrangement of the measuring points and their number had a significant impact on the results. At the same time, it should be noted that also in this case, the obtained errors of the reconstruction coefficient were at a satisfactory level. Only in five cases, the recreating error exceeded 1%.
We next verified the errors of the recreating function u in the measuring points for the recreated parameters α, β, c. In all cases, these errors were at a low level. A tendency can be seen that as input data noise increased, the recreating errors of the function u increased, although this was not always the case. A similar tendency applied to the mesh density; in general, the lowest reconstruction errors were obtained for the 160 × 160 × 250 mesh and the highest for the 80 × 80 × 100 mesh. Importantly, the average errors of the recreating function u turned out to be greater for the second arrangement of measurement points (seven measurement points) in ten of twelve cases (see Tables 5 and 6). Table 3. Results of the calculations in the case of twelve measurements points P i (i = 1, 2, . . . , 12): α ireconstructed value of x-directional order derivative α, β i -reconstructed value of y-directional order derivative β, c-reconstructed value of parameter c, δ-relative error of reconstruction, J-value of the objective function, σ-standard deviation of the objective function.

Mesh Size
Noise     The greater the noise on the input data was, the greater these errors were. For exact data, the reconstruction errors were at a stable, low level. In the case of point K6, it can be seen that the reconstruction error for input data disturbed with a 5% noise error at the end of the process was greater than the corresponding error for a noise of 10%.  Now, we present a way of achieving the solution by the ACO algorithm. Figures 6 and 7 show the values of the objective function in successive iterations of the algorithm for the case of the 100 × 100 × 200 mesh, seven measurement points, and various input data noise values. The values of the objective function for the first few iterations were relatively high; therefore, due to the scale, they are not marked in the figures. It can be seen that around 15 iterations, these values stabilized, and in the last iterations, the algorithm slightly improved the searched minimum. The selection of the parameters for the ant algorithm is quite important and was obtained based on the experience of the authors [40,44].  Finally, we present the influence of the recreated parameters on the values of the function u at the measuring point K4. For this purpose, we calculated the values of the function u at the point K4 for different values of the parameters α, β, c. Firstly, we determined the values of the parameters β, c, and we changed the values of the parameter α = 0.2, 0.4, 0.6, 0.8. Then, we similarly executed with the parameter β at fixed values of the α, c. Finally, we set the α, β parameters, and the parameter c took the values 500, 700, . . . , 1500. As can be seen in Figures 8-10, a change in the value of the β parameter had a greater impact on the obtained values of the function u at point K4 than in the case of the α derivative order. Furthermore, introducing a change in the value of the parameter c had a significant impact on the function, especially in the case of lower values of the parameter c. It can also be observed that changing the α derivative parameter affected the u function to a lesser extent than the other parameters.

Conclusions
The paper presented a numerical solution of a two-dimensional differential equation with a fractional derivative with respect to a space of the Riemann-Liouville type. This type of equation can describe the phenomenon of heat conduction in porous media, which the authors will also consider in further research. The presented numerical solution was also used to solve the inverse problem consisting of recreating the orders of the derivatives α, β and the parameter c. The presented results of the numerical experiment seemed satisfactory, and the algorithm was stable even for input data burdened with noise (measurement error).
The presented algorithm was characterized by faster calculations in relation to the classic differential scheme, resistance to falling into local minima compared to other minimization methods, low requirements for the objective function (it was not necessary to require continuity, differentiability, or convexity), and the stability against input data errors. The above-described properties can be treated as its advantages. The only downside was the computation time of the ACO algorithm.
In the further part of the research, the authors plan the following: • continued tests of the proposed algorithm (reconstructing of more parameters); • application of the model described in the study (differential equation) to model the heat flow process in porous materials; • development of the algorithm to shorten the computation time; • investigation of the influence of initial conditions (due to the fact that fractional derivatives contain the memory of past events), and thus the development of the model with the Caputo derivative to the time variable.