Two unconditionally stable difference schemes for time distributed-order differential equation based on Caputo–Fabrizio fractional derivative

We consider distributed-order partial differential equations with time fractional derivative proposed by Caputo and Fabrizio in a one-dimensional space. Two finite difference schemes are established via Grünwald formula. We show that these two schemes are unconditionally stable with convergence rates O(τ2+h2+Δα2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$O(\tau ^{2}+h^{2}+ \Delta \alpha ^{2})$\end{document} and O(τ2+h4+Δα4)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$O(\tau ^{2}+h^{4}+\Delta \alpha ^{4})$\end{document} in discrete L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$L^{2}$\end{document}, respectively, where Δα, h, and τ are step sizes for distributed-order, space, and time variables, respectively. Finally, the performance of difference schemes is illustrated via numerical examples.


Introduction
The concept of fractional derivative can be traced back to the letter of Leibniz to L'Hospital in 1695. Mathematician Euler discovered the Γ function in 1729, which established the foundation of fractional derivative. Compared with the classical integer-order models, fractional derivatives provide a more profound and comprehensive explanation about memory, heredity, non-locality of complex phenomena and processes. At present, the most commonly used definitions of fractional derivatives are the Grünwald-Letnikov derivative, the Riemann-Liouville derivative, the Riesz derivative, and the Caputo derivative [1]. Some new fractional derivatives, such as the Caputo-Fabrizio derivative and the Atangana-Baleanu fractional derivative, have tremendously promoted the capability of modeling complex physical phenomena and processes. The connatural behavior has been analyzed for evolution equations generated by three fractional derivatives namely the Riemann-Liouville, Caputo-Fabrizio, and Atangana-Baleanu fractional derivatives [2][3][4]. The Caputo-Fabrizio derivative is defined via an integral operator without singular kernel [5,6], and it has supplementary motivating properties compared to the others. For example, it can describe substance heterogeneities and configurations with different scales, which noticeably cannot be characterized by the local theories. The properties of numerical approximation are discussed for the initial value problems with Caputo-Fabrizio fractional derivative operator [7,8]. Several valuable applications to the real world data, such as blood ethanol concentration system, chickenpox disease model, and dengue fever outbreak model, have been discussed; the discussions were based on the fractional derivatives in the Caputo sense, Caputo-Fabrizio sense, and Atangana-Baleanu-Caputo sense, respectively [9][10][11][12][13][14][15].
Usually, it is difficult to obtain analytical solutions of these models, so different numerical schemes for fractional partial differential equations(FPDE), which can exhibit history dependence and nonlocal features, have been developed. The common numerical methods for fractional partial differential equations include finite-difference [37][38][39][40][41], finiteelement methods [42,43], finite volume schemes [44][45][46], mixed finite element schemes [47], spectral/spectral-element schemes [48,49], the decomposition method [50,51], and others [52]. In [53][54][55][56], numerical analysis for distributed-order FPDEs was extensively investigated. More recently, Tomovski and Sandev [57] expressed the solutions of generalized distributed-order diffusion equations equipped with fractional time derivatives by use of the Fourier-Laplace transform. The authors [58] obtained the existence about Cauchy problems for the diffusion equations with time distributed-order derivative, they also computed it by Laplace transform and Fourier transform. Luchko [56] discussed the dependence of the uniqueness and continuity of the generalized time fractional diffusion equations on initial conditions. However, there are few numerical studies about the distributed-order equations, especially the study about high-order schemes is almost blank. It has motivated us to find efficient numerical schemes for distributed-order FPDEs.
Consider the following distributed-order fractional diffusion equation equipped with the fractional derivative developed by Caputo and Fabrizio [5]: where the distributed-order derivative in the Caputo-Fabrizio sense is defined as follows: and 0 < v 1 < v 2 < 1, Ω = (a, b), and f (x, t) represents the source term. This paper has the following organization. Section 2 provides some preliminary work, including the space partition, the inner products, and norms. Section 3 is devoted to the finite difference scheme, which is proved to be second-order with respect to distributedorder, space, and time respectively, and the stability and convergence rates are also studied. In Sect. 4, a compact finite difference method is proposed, the stability and convergence are rigorously proved. Section 5 provides numerical tests that illustrate the reliability of theoretical analysis. In Sect. 6 we summarize this paper and indicate the possible work in the future.

Preliminary
In order to establish finite difference methods, the space interval [a, b] is partitioned by and the average operator It is easy to see where I denotes the identical operator. We also denote Av = (Av 1 , Av 2 , . . . , Av M ) for v = (v 1 , v 2 , . . . , v M ), and A(u, v) = (Au, v) with (·, ·) the discrete L 2 inner product defined below.
, the discrete inner products and norms are defined as follows: and u 2 2 = (u, u), By summation by parts and the boundary conditions, it is easy to get that For the average operator A, we also define } the space of all grid functions with respect to Ω τ . Given a grid function v ∈ V τ , we denote the difference quotient in time direction as

Lemma 1 ([59]) The following inequalities hold for any grid function
Throughout this paper, the symbol C will indicate a genetic constant that depends on the function u and may assume different values at different occurrences, and the symbol O(S) indicates a quantity bounded above by S with a constant.

A second-order scheme for space and distributed order
Differential equation (1) has the form at the node (x i , t n ) where Noting that u(x, 0) = 0, we have and Substituting (4) and (5) into (2), we get and there exists a constant C satisfying The boundary and initial conditions show that Canceling the infinitesimal part and replacing the true solution u k i with approximate solution U k i , we get the finite difference method as follows:

Stability analysis
We give a lemma about G α l k which is useful in the stability analysis.

Lemma 3 ([61]) From the definition of G
Rearranging equation (8), we get Denote Multiplying both sides of equation (9) by hU n i and making summation with respect to i from 1 to M -1, we obtain Noting Lemma 1 and Young's inequality, we have Using the triangle inequality, we obtain Consequently, we get Namely, Then it follows that Making summation with respect to i from 1 to N shows that Observing that U 0 = 0 implies Q(U 0 ) = 0, we have Theorem 1 For scheme (8), stability inequality is as follows:
Making summation with respect to n from 1 to N shows that Theorem 2 For scheme (8), the following stability inequality holds:

A fourth-order method for space and distributed order
The following lemma is necessary for establishing a scheme with fourth-order accuracy in spatial variable, The differential equation at the node (x i , t n ) has the form
Noting the boundary and initial conditions, we have Canceling the infinitesimal part (r 2 ) n i and replacing the exact solution u k i with the approximate solution U k i , a compact difference method is obtained as follows:

Stability analysis
Sorting equation (16), we get the following equation: We set η 2 = α 2J l=0 d l w(α l ) 1 α l τ . Multiplying both sides of equation (17) by hU n i and making summation with respect to i from 1 to M -1 leads to Noting Lemma 1 and Young's inequality, we have Using the triangle inequality, we obtain It follows that Transposition leads to Making summation from n = 1 to n = N shows that Observing that u 0 = 0 implies Q(u 0 ) = 0, then we have Theorem 3 For scheme (16), the following stability inequality holds:

Convergence analysis
Combining equations (14), (15) with (16), an error equation can be obtained as follows:  Using the triangle inequality, we obtain Rearranging the above equation leads to Finally we get Denote Q e n = η 2 2 n k=1 G α l n-k e k 2 A .
Make summation with respect to n from 1 to N to get Observing that e 0 = 0 implies Q(e 0 ) = 0, it follows that Theorem 4 For scheme (16), the following stability inequality holds:

Numerical tests
We perform numerical experiments to verify the proposed numerical format. Here, we consider the convergence speed of the numerical solution format. For the problem to be solved, we take the domain Ω = [0, π], T = 0.5. The numerical simulations are performed with MATLAB2015 on a 8GB memory computer. Consider time distributed-order model as follows.
Firstly, the numerical convergence orders of the two difference formats are tested by Example 1. Let the spatial step size h and the distributed-order subinterval size α be fixed small enough. We take h = π/900, α = 1/300, and τ = 1/4, 1/8, 1/16, 1/32 to observe the convergence rates about time step τ for both of the two difference formats. Table 1 Table 2. From the table we can see that, when the spatial step size h, the discrete step of the distributed order α, and the time step τ are halved simultaneously, the errors between the numerical solution and the true solution will be reduced to quarter in both of maximum norm and the discrete L 2 norm. These results show that the numerical convergence orders for difference format (8) are approximately 2 about time variable, space variable, and the distributed order. Table 1 The discrete L 2 error and convergence rates in time for schemes (8) and (16) Table 4 The discrete L 2 errors, maximum errors, and their convergence rates for scheme (8)   In addition, we choose a group of mesh steps satisfying 8N = J 2 and M = floor(πJ), which means that τ varies synchronously with α 2 and h 2 . The simulation results of numerical format (16) for Example 1 are given in Table 3. We observe that when the spatial step h and the distributed-order step size α are halved, the corresponding numerical errors in discrete norm are reduced to 1/16 times. This means both of the convergence rates for spatial variable and distributed order are of fourth order and thus the time convergence rate is of second order. (1), if we give the true solution u(x, t) = 2t sin(x) with w(α) = α 2 , accordingly, f (x, t) can be derived as follows:

Example 2 For equation
For Example 2, we set the distributed-order interval [v 1 , v 2 ] = [0.1, 0.9] and the numbers of discrete points are required to satisfy M = N = J for scheme (8). The simulation results in Table 4 show that scheme (8) has convergence rates of second order for all of the spatial step h, the temporal step τ , and the distributed-order step α.
In addition, we choose a group of mesh steps satisfying 8N = J 2 and M = J, which means that τ varies synchronously with α 2 and h 2 . The simulation results of numerical format (16) for Example 2 are given in Table 5. We observe that when the spatial step h and the distributed-order step size α are halved, the corresponding numerical errors in discrete norm are reduced to 1/16 times. This means both of the convergence rates for spatial variable and distributed order are fourth order, and accordingly the convergence rate for time variable is second order.

Conclusion
In this paper, two effective finite difference schemes are developed for solving time distributed-order partial differential equations equipped with the new Caputo-Fabrizio fractional derivative in a one-dimensional space. The first scheme (8) is based on secondorder central difference quotient in space direction and the other one (16) is based on compact difference. We have proved that both of the two finite difference schemes are unconditionally stable. Theoretical analysis shows that the first scheme has convergence rates of second order for all of the spatial step h, the temporal step τ , and the distributed-order step α; and the second scheme (compact difference) has convergence rates of fourth order for both of the spatial step h and the distributed-order step α, while second order for the temporal step τ . Finally, some numerical experiments, which demonstrate the performance of the schemes and verify the correctness of the theoretical results, are carried out.
The basic model is in a one-dimensional space, an extension to multi-dimensional models is possible and feasible. With the development of the theory about fractional calculus, fractional partial differential equations with variable order become the potentially powerful tool for describing complex phenomena [63]. We hope to extend the idea of this paper to this kind of fractional equations in the future. Moreover, because of the nonlocal property of fractional derivatives, the requirement for memory and computational complexity will increase rapidly along with the size of discrete systems. So we plan to develop some techniques of fast solution in the future work.