The Numerical Solutions of a Two-Dimensional Space-Time Riesz-Caputo Fractional Diffusion Equation

This paper is concerned with the numerical solutions of a two-dimensional space-time fractional differential equation used to model the dynamic properties of complex systems governed by anomalous diffusion. The space-time fractional anomalous diffusion equation is defined by replacing the second order space and the first order time derivatives with Riesz and Caputo operators, respectively. Using the Laplace and Fourier transforms, a general representation of analytical solution is obtained in terms of the Mittag-Leffler function. Grünwald-Letnikov (GL) approximation is also used to find numerical solution of the problem. Finally, simulation results for two examples illustrate the comparison of the analytical and numerical solutions and also validity of the GL approach to this problem.


Introduction
In recent years, Fractional Order Partial Differential Equations (FOPDEs) have played an important role in applied mathematics, physics, engineering, economics and other research areas.The dynamic models of a growing number of physical phenomena in elasticity theory, control theory, finance, electromagnetics, diffusion processes and in many other fields can be modeled by FOPDEs which are characterized by fractional space and/or time derivatives.Several analytical and numerical solution methods have been introduced to obtain the solutions of FOPDEs.
The space-time fractional differential equations considered in this work have had more attention from the view point of physical applications.For example, anomalous diffusion or dispersion processes can be modeled by these type of equations.Some related works in the literature are as follows.Zhang and Liu [1] studied the fundamental solutions of space-time Riesz fractional partial differential equations in terms of the Caputo fractional derivatives.Shen et al. [2] proposed the exact and numerical solutions of the Riesz advection-dispersion equation which is defined in terms of the Riemann-Liouville (RL) fractional derivatives, and they also analyzed stability and convergence of the numerical methods.Chen et al. [3] obtained the analytical and numerical solutions of a Riesz spacefractional reaction-dispersion equation by using the Laplace and Fourier transforms and also an explicit finite-difference approximation.Yang et al. [4] presented the analytical solution and numerical schemes for the time and space-symmetric fractional diffusion equation under homogeneous Dirichlet and Neumann conditions.Zhuang et al. [5] used the numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term.Meerschaert and Tadjeran [6] also developed some basic numerical methods to solve one-dimensional fractional advectiondispersion equations with variable coefficients.Yang et al. [7] analyzed the numerical solutions of a fractional partial differential equation with the Riesz space derivatives on a finite domain by using three numerical methods.Ciesielski and Leszczynski [8,9] presented a numerical method based on finite-difference method for numerical solutions of ordinary differential equations in terms of Riesz and Riesz-Feller operators, respectively.Gorenflo et al. [10] considered a Cauchy problem for a partial differential equation of fractional order in terms of the Caputo time derivative and the Riesz space Pseudo-differential operators.Saichev and Zaslavsky [11] considered a generalization of the normal diffusion equation which can be used to formulate a fractional kinetic equation.
Povstenko [12] presented a survey of nonlocal generalizations of the Fourier law and analyzed a heat conduction equation in terms of time fractional derivative and Riesz space fractional derivative.For solution of the multidimensional fractional partial differential equations, Meerschaert et al. [13] presented a practical numerical method which is improved by using a shifted Grünwald finite difference scheme.El-Sayed and Gaber [14] applied the Adomian Decomposition Method to solve some diffusion-wave equations using the connection between the Caputo and Riesz fractional derivatives.Özdemir et al. [15] researched the solution of an axial-symmetric diffusion-wave problem in polar coordinates.
Özdemir and Karadeniz [16] applied the GL approximation to axial-symmetric diffusionwave problems defined in cylindrical coordinates.There are some other papers related to the solution of diffusion-wave equations concerning analytical and numerical solutions [17,18,19,20,21].
In this work, our aim is to find the analytical and numerical solutions of an initial value problem defined in terms of space-time fractional differential equation which is used to model diffusion process.We suppose that the solution and initial condition functions belong to the Lizorkin space.The main reason for this choice is that the Lizorkin space is convenient while dealing both with the Fourier transform, and with the fractional integration and differentiation operators.Moreover, this operator is invariant with respect to the fractional integration and differentiation operators [22].
After all, we can underline this work as follows.This work is an extension of a onedimensional problem considered in [1] a two dimensional case.In addition, we not only find the analytical solution but also study the numerical solution using an approximation of the Caputo fractional derivative.Since it is not always possible to find the analytical solution of fractional differential equations, it is important to have numerical approximations of the solutions.
The paper is organized as follows.In Section 2, some basic mathematical definitions from fractional calculus which are used for formulation of the problem are given.In Section 3, the detailed description and formulation of the problem is given.To find the analytical solution Laplace and Fourier transforms are used.In Section 4, the GL approximation is applied to obtain the numerical results.The applicability of this approach is shown in two-and three-dimensional figures in Section 5.Moreover, the variation of problem parameters are analyzed in these figures.Finally, we conclude our work in Section 6.

Preliminaries
In this section, we recall some basic definitions of fractional calculus used in our problem formulation.
There are some preferences dealt with in the definition of fractional derivatives such as Riemann-Liouville, Caputo, Riesz, Weyl, etc [23,24,25,26].It is noted that each type of fractional derivative has different properties and applications.This allows someone to decide and choose the proper definition for problem construction.Here, we give the definitions of the Riesz Pseudo-differential operator and the Lizorkin space.Furthermore, we recall the Fourier transform, the Caputo fractional derivative and the Laplace transform of the Caputo definition.Definition 1.Let S be the space of rapidly decreasing test functions.Denote by V (R) the set of functions v ∈ S satisfying the conditions: The Lizorkin space Φ (R) is introduced as the Fourier pre-image of the space V (R) in the space S , i.e., According to the definition of the Lizorkin space, any function ϕ ∈ Φ (R) satisfies the orthogonality conditions For a function u of the class S of a rapidly decreasing test functions on the real axis R, the Fourier transform is defined as whereas the inverse Fourier transform has the form Definition 3. The Riesz pseudo-differential operator is defined by analytic continuation in the whole range n − 1 < α ≤ n, n ∈ Z and α = 1 as where the coefficient c = , and are the left and right-side Weyl fractional derivatives.In Eqs.( 6) and ( 7), the I β ± (β > 0) denote the Weyl fractional integrals defined as Definition 4. The Caputo fractional derivative is defined as The Laplace transform of the Caputo fractional derivative of order 0 < β ≤ 1 is where U (s) is Laplace transform of u (t) .

Analytical Solution
Let u (x, y, t) be a real-valued, multivariable function whose arguments (x, y) are the space and t ∈ R + is the time variables.We consider the anomalous diffusion problem which is described by a Riesz-Caputo fractional partial differential equation as with the following initial and boundary conditions x,y→±∞ u (x, y, t) = 0.
It should be noted that we understand the space and time fractional derivatives in a sense of the Riesz and Caputo definitions, respectively.Furthermore, we assume that the solution and initial condition functions can be expanded into the complex Fourier series.We suppose that the solution of the problem is represented by the following series where i = √ −1.Firstly, we compute ∂ α u(x,y,t) ∂|x| α by using the definitions given in Section 2. For this purpose, it is necessary to calculate the left and right side RL fractional derivatives with respect to x variable.The calculation of −∞ D α x u (x, y, t) and x D α ∞ u (x, y, t) is given as follows (see also [1]): and Therefore, we obtain In a similar way, the term ∂ α u(x,y,t) ∂|y| α can be evaluated as follows: Substituting Eqs.( 13) to (3) into Eq.(10), we obtain Next, we take the Laplace transform with respect to time.Then we have (14) where Let us consider the main value branch of A for computational purposes and using the inverse Laplace transform, Eq.( 14) yields where E β is the Mittag-Leffler function.According to our assumptions and similar to the solution function u (x, y, t), the initial condition function u 0 (x, y) can also be expanded into a Fourier series: where u 0nm (0) are the Fourier coefficients On the other hand, we take t = 0 in Eq.( 13) and obtain k nm (0) = u 0nm (0).Then the closed form of analytical solution to the problem is obtained as

GL Approximation For Numerical Solution
To find the numerical solution of the problem we use the Grünwald-Letnikov (GL) formula which is based on time discretization and is defined as a limit of a fractional-order backward difference where β > 0 is the order of fractional derivative, t is the length of time interval on which problem is formulated, h is the length of subintervals of time t.From the theoretical background of fractional calculus it is well known that when the function f (t) of positive argument has continuous derivatives of integer order 0, 1, ..., n, the RL and Grünwald-Letnikov definitions are equivalent [25].Furthermore, the relation between the Caputo and RL fractional derivatives give allows us to apply the GL formula for approximation of the Caputo fractional derivative.It was mentioned in the problem formulation that we consider the Caputo time fractional derivative.Therefore, we give firstly the relation between the RL and Caputo definitions (see [27]) Substituting Eq. ( 18) into ( 19), we obtain the numerical calculation formula for the Caputo fractional derivative in terms of Grü nwald-Letnikov definition It should be noted that our aim is to substantiate the approximate numerical solution of the problem.For this purpose we consider a finite number series to calculate the approximate value of the Caputo fractional derivative: For analyzing this approach let us consider the following problem: one-term, constant coefficient fractional differential equation under initial condition where . Now, we use the approximation of the Caputo derivative given by Eq. ( 4) and rewrite Eq. ( 21) as where M = t h and w The applicability of this numerical approximation is evident from Figures obtained by MATLAB.In the following section we present some Figures and validate our analytical and numerical solutions from different points of view.

Numerical Examples
In this section, we study numerical examples to demonstrate the applicability of the GL approximation and investigate the behavior of anomalous diffusion process for different values of parameters.For this purpose, we present two-and three-dimensional graphs.First we consider the exponential initial condition function The function u 0 (x, y) is continuated periodically on the whole axis.In the second example, the parabolic initial condition is chosen.As above, this function is continuated periodically on the whole axes.At first, we point out the overlapping of the analytical and numerical solutions for the order of time derivative β = 1 and the order of space derivative α = 0.5, the term numbers n = m = 10, the step size h = 0.01 at x = 0.2 and y = 0.5 in Figures 1a and 1b.The analytical and numerical solutions are found to be in good agreement.Figures 2a and 2b show the solution profiles for different values of β = 0.1, 0.5, 0.7, 1 for α = 0.5, x = 0.2, y = 0.5 and h = 0.01.
Similarly, Figures 3a and 3b display the variation of behavior of the solution for the values α = 0.1, 0.5, 0.7, 0.99 when β = 0.5, x = 0.2, y = 0.5 and h = 0.01.It should be recalled that α = 1 due to the definition of the Riesz operator.
In Figures 4a and 4b, we show the contribution of term numbers n = m = 5, 10, 20, 30, 35 to the solution for α = 0.99, β = 1, x = 0.2, y = 0.5 and h = 0.01.We obtain that the numerical solutions coincide for more than 35 terms.In other words, it is enough to take only 35 terms for solution series defined by Eq.( 3).Finally, we present the whole solution of the problem in Figures 5a and 5b for α = 0.99, β = 1, h = 0.01 and n = m = 35.In addition, we show the dependence of u (x, y, t) with respect on x and t variables in Figures 6a and  6b.

Conclusion
In this paper, the solutions of a twodimensional anomalous diffusion problem were researched.The solutions of the problem were assumed to have Fourier series expansion.In the problem formulation, the space derivatives were considered as the Riesz pseudodifferential operators, and the time derivative was defined as the Caputo fractional derivative.To obtain the closed form analytical solution of the problem, the Laplace transform with respect to time variable and the Fourier transform with respect to space variables were applied.The GL approximation of the Caputo fractional derivative was used to obtain the numerical solutions.The applicability of the GL numerical scheme was substantiated by two examples and is evident from the Figures obtained with MATLAB.In addition, the solutions for different values of parameters were shown in Figures.Finally, it should be pointed out that the fundamental solution of the considered anomalous diffusion problem, which behaves as subdiffusion process, can be computed easily and effectively by the GL approximation.