Semianalytic Solution of Space-Time Fractional Diffusion Equation

We study the space-time fractional diffusion equation with spatial Riesz-Feller fractional derivative and Caputo fractional time derivative.The continuation of the solution of this fractional equation to the solution of the corresponding integer order equation is proved.The series solution of this problem is obtained via the optimal homotopy analysis method (OHAM). Numerical simulations are presented to validate themethod and to show the effect of changing the fractional derivative parameters on the solution behavior.


Introduction
Fractional derivatives have been utilized for describing numerous models in different branches of applied science.Anomalous diffusion is one process that has been successfully modeled by fractional derivatives.This type of diffusion is characterized by the nonlinear dependence of the mean square displacement () of a diffusing particle over time :  2 () ∝     , and it can be interpreted as the Lévy stable densities.On the other hand, in the case of classical diffusion, linear dependence  2 () ∝  occurs and it follows Gaussian statistics and Fick's second law for running processes at time  [1].
Anomalous diffusion is described by fractional partial differential equations (FPDEs) by replacing classical derivatives with derivatives of fractional order.A type of anomalous diffusion that received great attention is the fractional diffusion equation with spatial Riesz and Riesz-Feller fractional derivatives [2,3].These fractional derivative operators are used to model phenomena as anomalous diffusion [1], discrete random walk models [3,4], and continuous random walk models [5] in the form of FPDEs.
General solutions are obtained for some classes of linear FPDEs, using integral transforms.For example, in [6], a general solution is given for a time fractional diffusion-wave equation defined in a bounded space domain.The fractional time derivative is described in the Caputo sense.The finite sine Laplace transform technique is used to obtain the solution.In [7], analytic techniques are introduced to solve certain types of multiterm time-space Caputo-Riesz fractional advection-diffusion equations on a finite domain using the equivalent relationship between the Laplacian operator and the Riesz fractional derivative.In [8], a fundamental solution is investigated for the space-time fractional diffusion-wave equation, which is obtained from the standard diffusion equation by replacing the second-order space derivative with a Riesz-Feller derivative and the first-order time derivative with a Caputo derivative.The paper [9] presents a solution of a unified reaction-diffusion equation with the Riemann-Liouville fractional derivative as the time derivative and Riesz-Feller derivative as the space derivative in terms of -function.In [10], a closed form solution in terms of Mittag-Leffler function is obtained for a space-time fractional telegraph equation with Riesz-Feller derivative as the space derivative.In [11], a solution of a space-time fractional advection-dispersion equation is derived in terms of the Green functions.
The proposed abovementioned methods work in case of linear FPDEs.However, for the nonlinear case, the numerical methods and the semianalytic approaches have the advantage.Several numerical techniques are introduced in the literature, and we state some of the known methods here.In [12], the matrix transform method with the method of lines is used to obtain the numerical solution of a fractional partial differential equation with Riesz space fractional derivatives on a finite domain.An explicit finite-difference approximation is used to solve the variable-order nonlinear space fractional diffusion equation [13].Crank-Nicolson method for the fractional diffusion equation with the Riesz fractional derivative is used in [14], while the improved matrix transform method with Pade approximation is used in [15].Yet, semianalytic solutions for this type of problems include only the series solution by the variational iteration method approach (see [16,17]) and by the homotopy analysis method (HAM) (see, e.g., [2,18]).In [19], the homotopy analysis method with the Fourier transform is used to obtain and analytically approximate solutions of Riesz fractional diffusion equation and Riesz fractional advection-dispersion equation.
As seen in the literature review, only few articles dealt with applying iterative techniques to these types of FPDEs due to the difficulty of repeated application of Riesz or Riesz-Feller fractional derivatives to the solution components.In this work, we consider the space-time fractional diffusion equation of the following form: where     denotes the Caputo fractional derivative (in time) of order  and  ,  denotes the Riesz-Feller fractional derivative (in space) of order  and skewness .Parameters , , and  are real numbers restricted to and the two functions  and  are continuous functions in .We aim to establish the continuation of the solution of (1) to the exact solution of the corresponding equation in Riesz fractional derivative as the skewness parameter approaches zero.This objective is carried out theoretically then via approximate series solution obtained iteratively by applying the OHAM.This paper is organized as follows.In Section 2, basic definitions of fractional derivative operators involved are presented.Proof of solution continuation is presented in Section 3. The OHAM is illustrated in Section 4. The results of numerical experiments are presented in Section 5. Section 6 contains the conclusion of this work.

Fractional Derivatives and Integrals
Definition 1.A real function (),  > 0, is said to be in the space   ,  ∈ R, if there exists a real number  > , such that () =    1 (), where  1 () ∈ (0, ∞), and it is said to be in the space    if   ∈   ,  ∈ N.
Definition 5.The Riesz-Feller partial fractional derivative  ,  is defined as [3] where   ± () are the Weyl fractional derivatives defined as International Journal of Differential Equations 3 and   ± denote the Weyl fractional integrals of order  > 0, given by The coefficients  + and  − have the following forms: where || ≤ min (, 2 − ).When  = 0, the Weyl fractional derivative degenerates into the identity operator: For continuity, we get Evidently, in the case of  = 2 and  = 0, we define For the case of  = 1 and  = 0, we have where  is the Hilbert transform and the integral is understood in the Cauchy principal value sense.For the cases  = 1 and  = ±1, we obtain the first derivative: The special case  = 0 is known as Riesz fractional derivative.

Optimal Homotopy Analysis Method (OHAM)
In this section, we begin by illustrating the basic steps of the classical homotopy analysis method (HAM).Consider the following nonlinear equation: where  is a nonlinear operator, (, ) is the unknown function, and  and  denote spatial and temporal independent variables, respectively.Liao [20] generalized the traditional homotopy method to construct the so-called zeroorder deformation equation given by where  ∈ [0, 1] is an embedding parameter, ℏ is a nonzero auxiliary parameter, (, ) is an auxiliary function,  is an auxiliary linear operator,  0 (, ) is an initial guess of (, ), and (, ; ) is an unknown function.Obviously, when  = 0 and  = 1, we have (, ; 0) =  0 (, ) and (, ; 1) = (, ), respectively.Thus, as  increases from 0 to 1, the solution (, ; ) varies from the initial guess  0 (, ) to the required solution (, ).By expanding (, ; ) in Taylor series with respect to , we have where If the auxiliary linear operator, the initial guess, and the auxiliary parameter ℏ and the auxiliary function are so properly chosen, then, as proved by Liao [20], series (17) converges at  = 1 and one has which, as proven by Liao [20], must be one of the solutions of the original nonlinear equation.Using definition (18), the governing equation of the HAM can be deduced from the zero-order deformation equation ( 16) as follows.Define the vector ⃗   = { 0 (, ) ,  1 (, ) ,  2 (, ) , . . .,   (, )} . (20) From ( 16), the so-called th-order deformation equation is given by where ( Applying the inverse operator  −1 to both sides of (21),   (, ) can be calculated by symbolic computations software.The HAM has been successfully applied to solve various classes of equations and applied problems [21][22][23][24][25].
In the classical HAM, choosing the value of parameter ℏ depends on inspecting the graph of the quantity of interest: the solution or one of its derivatives.Yet, when (, ) is fixed, it is obvious that   (, ) contains only one control parameter ℏ.Thus, by constructing a formula for the residual error, the OHAM solution is obtained by choosing the value International Journal of Differential Equations for parameter ℏ that minimizes the error.Here, the averaged residual error defined for ordinary differential equations in [26] is generalized to the case of two variable partial differential equations in the following form: which is a nonlinear algebraic equation of one unknown: the convergence-control parameter ℏ.Thus, the optimal value of ℏ is determined by the minimum of the averaged residual error   to ensure the fast convergence of the homotopy series.

Continuation of the Solution
In this section, we define Riesz-Feller fractional derivative on unbounded domain in Caputo sense.This definition is considered as a generalization of the definition given in [27] for Riesz fractional derivative on bounded domains.
Then, for  ∈ (0, 2),  ̸ = 1, Proof.Consider the case  ∈ (0, 1).Riesz-Feller fractional derivative can be written in the following form: which can be written as For   ,  , we have which by substitution yields which is the same as (29).A similar argument holds for the case  ∈ (1, 2).
As the equivalence between the Riemann and Caputo definitions of Riesz-Feller fractional derivative is deduced, we discuss the continuation of the solution obtained.
Proof.Consider the set of functions: which is the exact solution of the Riesz fractional diffusion equation (35).
Similar arguments hold to prove the continuation of the solution when both parameters  and  tend to their integer limits.

Numerical Simulation
In this section, we consider linear and nonlinear problems to illustrate the efficiency of the OHAM to this type of problems and to illustrate the continuation of the solution we proved in Section 4. Yet, to apply the recursive technique of OHAM to problem (1), a repeated evaluation of Riesz-Feller fractional derivative to solution components is needed.This is accomplished by using a property of Riesz-Feller fractional derivative that is proved in the following lemma.Lemma 10.Let  ∈ (0, 2),  ̸ = 1, || ≤ min (, 2 − ), and  > 0.Then, Proof.
Case  ∈ (0, 1).From the definition of Riesz-Feller fractional derivative and by substitution, we can write which reduces to Substituting for  + and  − , we have Utilizing Euler form and trigonometric identities, we have which reduces to Carrying out the same procedure, the lemma follows.
The solution (the first five terms in the OHAM series) behavior as the parameters , , and  change at a fixed time  = 5.0 and in the interval 0 ≤  ≤  is shown in Figures 1-3.The amplitude of the solution decreases  and/or  increases.As  reaches 2 and/or  increases up to unity, the series solution approximately coincides with the exact solution of the corresponding integer order equation  =  −( 2 / 2 −) sin (/), represented by the solid line in the graph (refer to Figures 2 and 3    in plots is the partial sum of the first five terms:  = 5 (summing  0 to  4 ).The optimal convergence parameter ℏ in each case is obtained by minimizing the residual error   displayed in (23) in the intervals 0 ≤  ≤  and 0 ≤  ≤ 6.0, and the obtained results are shown in Table 1.
Example 12. Consider the nonlinear fractional problem (1) with () = , () = 4 2 − , and () = sin (2): Here, we choose the auxiliary linear operator as and the nonlinear operator  is chosen as (61) We choose (, ) = 1 and  0 = sin (2).Following the same procedure in the previous example to obtain the series solution terms, and the solution is thus obtained as Figures 4-6 show the solution (the first five terms in the OHAM series) behavior as the parameters , , and  change at a fixed time  = 0.2 and in the interval 0 ≤  ≤ /2.As  and/or  increases, the amplitude of the solution decreases.As  tends to 2 and/or  goes up to unity, the series solution approximately coincides with the exact solution of the corresponding integer order equation  =  − sin(2), represented by the solid line in the graph (see Figures 5  and 6), and this confirms the continuation of the solution numerically.The series displayed in plots is the partial sum of the first five terms:  = 5 (summing  0 to  4 ).The optimal convergence parameter ℏ in each case is obtained by minimizing the residual error   displayed in (23) in the intervals

Table 1 :
The estimated optimal convergence parameter ℏ and the corresponding residual error   at different fractional parameters used in Example 11.