Theoretical Analysis (Convergence and Stability) of a Difference Approximation for Multiterm Time Fractional Convection Diffusion-Wave Equations with Delay

In this paper, we introduce a high order numerical approximation method for convection diffusion wave equations armed with a multiterm time fractional Caputo operator and a nonlinear fixed time delay. A temporal second-order scheme which is behaving linearly is derived and analyzed for the problem under consideration based on a combination of the formula of L2 − 1σ and the order reduction technique. By means of the discrete energy method, convergence and stability of the proposed compact difference scheme are estimated unconditionally. A numerical example is provided to illustrate the theoretical results.


Introduction
Fractional derivatives and integrals have recently gained high interest in many fields of science. The ability of classifying and capturing the memory and hereditary properties of various materials and processes is an advantage of fractional derivatives over their integer counterparts, e.g., the modeling of anomalous diffusion by fractional differential equations gives more informative and interesting models [1]. For time-fractional differential equations, the memory feature implies that all previous information is needed to evaluate the time fractional derivative at the current time level. Accordingly, designing a numerical differentiation formula of good accuracy is as ever paramount, but especially hard. The approximation formulas based on the interpolation approximation, such as L1 [2] and L 2 − 1 σ [3], are of significance to design numerical algorithms to solve time-fractional differential equations. Demonstrated applications in numerous seemingly diverse and widespread fields of physics, such as in porous and glassy materials, in percolation clusters over fractals to semi-conductors, polymers, random media, and beyond, like geophysical and biological systems or processes (e.g., [4,5]), can be effectively modeled by time-fractional diffusion-wave equations of different types. Here, we are seeking to design a compact difference scheme that behaves linearly to numerically solve the non-linear delayed multiterm time fractional convection diffusion-wave equation (dmfCDWEs) with spatial variable coefficients. More specifically, we consider m ∑ r=0 p r ∂ α r u(x, t) ∂t α r with the following initial and boundary conditions where s > 0 is a fixed delay parameter, q 1 (x) and q 2 (x) are functions chosen, respectively, to be sufficiently and arbitrary differentiable functions. The fractional derivative is defined in Caputo sense and the fractional orders {α r | 1 ≤ r ≤ m} are specified in the manner {1 < α m ≤ α m−1 < · · · < α 0 = 2}. The existence and uniqueness of the global mild solutions for the problem of nonlinear fractional reaction-diffusion equations with delay and Caputo's fractional derivatives are addressed in [6]. This work can be considered to be an extension of our previously published work [7], in which we discussed a single term time fractional wave equation with spatial constant coefficients. The scheme was of 2 − α order in time and fourth in space. Here, we treat the multiterm time fractional order case with spatial variable coefficients and seek to have temporal second order of convergence. Accordingly, we hinged on the proposed numerical formula in [8] to approximate the multiterm Caputo fractional derivatives of order α r (0 < α r ≤ 1) at the super-convergent point. The formula L 2 − 1 σ can achieve at least second-order accuracy at this point. We rely on the compact operator proposed in [9] in order to attain a fourth order accuracy with respect to space in case of having spatial variable coefficients.
There are some results and findings available regarding the theoretical analysis and numerical computation of single term time fractional sub or super diffusion equations with delay. In [10], the authors introduced a satisfactory numerical method for time fractional diffusion equations with delay. In [11], a novel discrete Grönwall inequality is used to simplify the analysis of difference schemes for time-fractional multi-delayed sub-diffusion equations. Convergence and stability of a compact finite difference method for nonlinear time fractional reaction-diffusion equations with a fixed delay are proposed in [12] by the aid of a new discrete form of the fractional Grönwall inequalities. A numerical solution for a class of time fractional diffusion equations with delay is proposed in [13] that is based on a smooth difference approximation of specific L1 type. Additionally, there are many difference and spectral approaches proposed for multiterm or distributed order time fractional differential equations. An efficient spectral method that is based on Jacobi-Gauss-Radau collocation is applied in order to solve a system of multi-dimensional distributed-order generalized Schrödinger equations in [14]. A combined difference and Galerkin-Legendre spectral method in [15] is used to solve time fractional diffusion equations with nonlinear source term. A Legendre spectral-collocation method for the numerical solution of distributed-order fractional initial value problems is designed in [16]. In [17], the authors proposed a spectral τ-scheme to discretize the fractional diffusion equation with distributed-order fractional derivative in time and Dirichlet boundary conditions. The model solution is expanded in multi-dimensions in terms of Legendre polynomials and the discrete equations are obtained with the τ-method. The two-dimensional distributed-order time fractional cable equation is numerically solved based on the finite difference/spectral method, as clarified in [18]. Two classes of finite difference methods that are based on backward differential formula discretization in the temporal direction are proposed in [19] to efficiently solve the semilinear space fractional reaction-diffusion equation with time-delay. The coefficients of the problem are constants, a fractional centered difference approximation is employed for the space fractional derivative [20], and it gains a fourth order approximation in space due to the use of a specific compact operator [21].
The main purpose of our work is the manufacturing of a difference scheme for problems of the kind of (1). Until now, few works hace paid close attention to the multiterm fractional wave equation with variable coefficients and delay argument simultaneously. It is well known that it may be a more challenging task to solve the fractional partial differential equation with delay effectively, since the evolution of a fractional partial differential equation with delay at time t not only depends on its value at t − s, but also depends on all previous solutions due to the non-locality of the fractional operator. Higher order numerical schemes are extremely scarce and difficult in regard to the analysis and implementation for the variable coefficient multiterm time fractional convection-diffusion wave equation with delay. For a single temporal term for that kind of problem, we refer to [22]. For the multiterm problem, we design the difference scheme at a super-convergent point to gain a high order of convergence [8] which is more challenging than the single term problem at the level of theoretical analysis.
In order to simplify the problem, an exponential transformation technique [23] can be applied to system (1) to avoid the difficulties resulting from the spatial variable when constructing high order compact difference methods. That technique is used to eliminate the convection term. Assume that When substituting (2)-(4) into (1), one directly obtains with the following initial and boundary conditions In order to transform (5) to a system with zero Dirichlet boundary conditions, we define h(x, t) := with the following initial and boundary conditions According to the new transformed system (6), we analytically overcame the first degree of complexity by the elimination of the convection term. There are still two degrees of complexity to be numerically overcome; the multiterm fractional order on the one hand and the nonlinear delay on the other. Throughout this work, we assume that the function f (µ, v, x, t) and the solution ν(x, t) of (6) are sufficiently smooth in the following sense: The partial derivativesf µ (µ, ν, x, t) andf ν (µ, ν, x, t) are continuous in the 0 -neighborhood of the solution. Define The structure of this paper is arranged, as follows. First, we introduce a derivation of the compact difference scheme. Next, in the third section, convergence and stability for the compact difference scheme are carried out. Finally, the paper ends with a numerical illustration and a conclusion.

A Compact Difference Scheme
A linearized numerical method that combines the super-convergence approximation L 2 − 1 σ with the order reduction method is derived. Some further notations are fixed before we continue. Take two positive integers M and n 0 , let h = L M , τ = s n 0 and denote x i = i h for i = 0, . . . , M; t k = k τ and t k+σ = (k + σ) τ, for k = −n 0 , . . . , N, where N = T τ . Using the points x i in space and t k in time, we cover the space-time domain by Here, we give a preliminary for the Alikhanov scheme. Denote γ r = α r − 1 (0 ≤ r ≤ m) and Additionally, now we invoke the following lemmas from Alikhanov work.
For the sake of simplicity, let σ = σ * here and later. Define [3] The next two lemmas are devoted to the properties of the coefficientsĈ (k) n andb n .

Lemma 2.7 ([26]
). For any grid function w 1 , Lemma 2.10 ( [24]). Suppose that ·, · * is an inner product onŴ h and · * is a norm deduced by the inner product. For any grid functions w 0 , w 1 , . . . , w k+1 ∈Ŵ h , we have the following inequality where Additionally, it holds that Let us initiate by order reduction for system (6) by letting Subsequently, the equivalent system to (6) after order reduction can be formulated as with the following initial and boundary conditions

Compact Difference Scheme Construction
Consider (20a) at (x i , t k+σ ),, and then we get such thatf From Lemma 2.10, we conclude that and a direct expansion of Taylor type yields such thatF Acting the averaging operator A on both sides of (21), noticing Lemma 2.10 and using Taylor expansion, we arrive at Next, by using Lemma 2.5 and (23), we obtain where a constantc 0 exists in order that By considering (20b) at (x i , t 1/2 ) and (x i , t k+σ ), respectively, operating by A on both equations, we obtain by the aids of Taylor expansions and Lemmas 2.5 and 2.9, which Moreover, there exists a constantc 1 > 0, such that By omitting the small terms in (26)-(28) and noticing the initial and boundary conditions, we construct a spatial fourth order difference scheme for problem (6), as follows

The Stability and Convergence of the Constructed Difference Schemes
First, we will start with some technical lemmas, which will be helpful in the context of convergence and stability.
The nonlinear delay termf (µ, ν, x, t) is sufficiently smooth and it satisfies the Lipschitz condition where c 1 and c 2 are two positive constants.
the following estimate is satisfied Proof. RecallingF k+σ i from (24) and under the assumptions of the nonlinear delay term (33), we can deduce the following estimate For the convenience of our analysis, the following Grönwall inequality is recalled.
Proof. Starting from the l.h.s of (40), then the r.h.s of (40) is achieved directly. Additionally, starting from the l.h.s of (41), then the r.h.s of (41) is immediately held. Invoking the previous estimates and Lemma 2.10, we deduce whereẼ k+1 is defined by (44) and so (42) is achieved. The estimate (43) is simply calculated. Now, we are in a position to combine the Lemmas 3.1-3.3 to prove the convergence and stability of the proposed compact difference scheme. To that end, Let Subtracting (31) from (26)-(28), (20c)-(20e), respectively, we obtain the error equations, as follows (51f) is the smooth solution of (5) and {V k i , ν k i |0 ≤ i ≤ M, −n 0 ≤ k ≤ N} the numerical solution of the scheme (31). Subsequently, there exist positive constants h 0 and τ 0 , independent of h and τ, such that, when h ≤ h 0 and τ ≤ τ 0 , we have the error estimate Proof. The proof will be preformed in two steps. Let us tackle the first one.

Almost Unconditional Stability
To discuss the stability of the compact difference scheme (31), we also use the discrete energy method.
where k i denotes an initial perturbation term that is very small.
Proof. The perturbation equations in terms ofρ k i andē k i come by subtracting (82) from (31) and similar to the proof of Theorem 1, the conclusion of stability holds immediately.

Generalized Scheme for the Distributed Order Case
We are now in a position to consider the distributed order form of dmfCDWEs, which means that with the following initial and boundary conditions Following the same manipulations illustrated before; starting from an exponential transformation technique, then a transformation to zero Dirichlet boundary conditions, we obtain the following system with the initial and boundary conditions as A numerical quadrature rule can be adapted to transform the system (85) to dmfCDWEs. We recall Simpson's rule (also known as the three-point Newton-Cotes quadrature rule), a proof of which can be found in any descent textbook. Lemma 5.1. Consider an equidistant partition of the interval [1,2] into 2J subintervals, let ∆α = 1 2J and denote α l = 1 + l ∆α, 0 ≤ l ≤ 2J. Subsequently, the composite Simpson's rule reads where Define the function G(· ; . Suppose that G(α) ∈ C 4 ([1, 2]), then by using Lemma 5.1, we approximate the distributed derivative as , n = 0, 1, . . . , k, then the constructed compact difference scheme has the following form (88f) Remark 1. The local truncation error of the compact difference scheme (88) for the distributed order system (85) is of order O τ 2 + h 4 + (∆α) 4 . The convergence and stability estimates can be derived in the same manner as in Theorem 1 and Theorem 2.

Numerical Illustration
The purpose of the present section is to demonstrate the convergence rate of the method. We will consider the maximum absolute error between the exact solution u(x i , t k ) of the continuous problem and corresponding approximations u k i , which is given by Moreover, we define the standard rates We consider the following multiterm time fractional delay sup-diffusion problem Note that g(x, t) is defined/derived, such that u(x, t) = e x x 2 (1 − x) 2 t 3 is the exact solution. The exact solution determines the initial condition and boundary conditions. The difference scheme (31) is employed in order to obtain the numerical solution. First, the numerical accuracy of this scheme in time will be verified. Taking a sufficiently small step size h and varying step size τ, the numerical errors and numerical convergence orders are listed in the lower half of Table 1. The computational results presented in Table 1 confirm the second-order convergence of the difference scheme (31) in time. Table 1. Absolute errors and standard convergence rates in space and time when approximating the solution u of (1) with (α 1 = 1.3, α 2 = 1.5, α 3 = 1.7), while using the difference method (31). The parameters and conditions employed in this case correspond to those in Example 6.

Spatial Analysis of Convergence
(p 0 = 1, p 1 = 1.25, p 2 = 2) (p 0 = 2, p 1 = 1.75, p 2 = 1) Next, the numerical accuracy of the difference scheme in space for solving this example is examined. The numerical results of this scheme for different step sizes in space are calculated and the numerical errors, as well as the numerical convergence orders are recorded in the upper half of Table 1. Again, from which, one can find that, in this case the fourth-order convergence is achieved.

Conclusions
A linearized difference scheme for solving a class of dmfCDWEs is constructed. With the help of an easy to execute and invertible exponential transformation, the considered problem can be converted into the delay variable coefficient fractional diffusion wave equation equivalently. Subsequently, we establish a fourth-order accurate numerical scheme that is based on a variable coefficient compact operator and with a temporal second order of convergence at a super-convergent point. The convergence and stability of the current numerical scheme are proved at length and a numerical example is finally added for the sake of demonstrating the theoretical findings.
Author Contributions: All authors contributed equally. All authors have read and agreed to the published version of the manuscript.
Funding: The first author wishes to acknowledge the support of RFBR Grant 19-01-00019.