Improved Gegenbauer spectral tau algorithms for distributed-order time-fractional telegraph models in multi-dimensions

The distributed-order fractional telegraph models are commonly used to describe the phenomenas of diffusion and wave-like anomalous, which can model processes without a power-law scale across the entire temporal domain. To increase the range of implementation of distributed-order fractional telegraph models, there is a need to present effective and accurate numerical algorithms to solve these models, as these models are hard to solve analytically. In this work, a novel matrix representation of the distributed-order fractional derivative based on shifted Gegenbauer (SG) polynomials has been derived. Also, two efficient algorithms based on the aforementioned operatonal matrix and the spectral tau method have been constructed for solving the one- and two-dimensional (1D and 2D) distributed-order time-fractional telegraph models with spatial variable coefficients. Also, a new operational matrix of the multiplication of space vectors has been built to have the ability in applying the tau method in the 2D case. The convergence and error bound analysis of the presented techniques are investigated. Moreover, the presented algorithms are applied on four miscellaneous test examples to illustrate the robustness and effectiveness of these algorithms.


Introduction
Single or multi-term fractional differential equations (FDEs) provided by spaceand/or time-fractional derivatives have attracted a lot of attention due to there efficiency in modelling many real-world phenomena arising in several fields (see [1][2][3][4][5][6]). Finding analytical or exact solutions of FDEs is difficult to obtain, so several numerical methods are proposed for solving different kinds of FDEs based on the finite difference, B-linear spline, second-order Newton interpolations, spectral method and Bernstein polynomials (see [7][8][9][10][11]). The distributed-order FDEs are a generalisation of the aforementioned two classes of FDEs, in which the fractional derivative order is integrated within a specified range. Subsequently, the distributeorder aggregates the contributions of the model order over the physical domain. As a result, the solutions of the distributed-order FDEs are not mathematically so direct. One of the main reasons for the interest in distributed-order operator is the relationship of this operator to physical processes that lack temporal scaling. Accelerating superdiffusion, random processes subject to Wiener's processes and decelerating subdiffusion are examples of processes that cannot be handled with fractional operator of specified order and must be modelled in terms of distributed operators [12][13][14].
The earliest results on the general solution for linear differential equations of distributed-order were published by Caputo [15]. The distributed-order timefractional diffusion equation (DOTFDE) explains a sub-diffusion random mechanism subjugated with the Wiener's process and has an exponent of decreasing diffusion over time [16]. Atanackovic et al. [17] modeled and analysed waves in a finite length viscoelastic rod by using restricted domain fractional diffusion-wave equations of distributed-order. They acquired stress and displacement in a stress relaxation test by prescribing boundary conditions on displacement. Eab et al. [18] presented the distributed-order Langevin-like equations of fractional order for characterising anomalous diffusion in the absence of scaling exponent or a unique diffusion and applied it for modelling some phenomena. Luchko [19] established the uniqueness of the generalised bounded domains DOTFDE and its continuous sensitivity on initial conditions. Meerschaert et al. [20] introduced stochastic analogues and explicit strong solutions for bounded domains DOTFDE. The fractional telegraph equation has significant for modelling wave-like anomalous and diffusive phenomena. In addition to describe the models of the blood flow under magnetic and vibration mode in a porous tube with thermochemical properties [21], enhancing denoising of image structure of elastic manifolds [22], bimolecular reactions in homogeneous media [23] and anomalous diffusion [24]. So it has many applications in various scientific fields. The single distributed-order time-fractional diffusion-wave equation (DOTFDWE) is generalised to the distributed-order time-fractional telegraph equation (DOTFTE) when the two distinct distributed-order differentiation on the intervals 1 2 and 0 1 are presented simultaneously. The DOTFTE can be associated with super/subdiffusive processes for specific density function choices. Vieira et al. [25] presented the fundamental solution of the DOTFTE in the n-dimensional space by using Fox H-functions via a combination of the Mellin, Fourier and Laplace transforms. Also, the provided techniques are generalisation of approaches for DOTFDWE that have been employed by numerous authors (see [26,27]).
The approximate solutions of distributed-order FDEs have a growing interest to study the behaviour of these equations while finding the exact analytical solutions of such equations are almost impossible [28][29][30]. In [31], Ye et al. presented a compact difference scheme for the DOTFDWE and proved the convergency and unconditional stability of the presented method. Abbaszadeh and Dehghan [32] developed a meshless method combined the alternative direction implicit procedure with the interpolating element-free Galerkin method for the DOTFDWE. Gao et al. [33] constructed two difference schemes for the 1D time-fractional wave equations of distributed-order where the time fractional derivatives approximated by using a weighted shifted Grünwald formula. Shi et al. [34] developed a completely discrete technique combined a finite element method on an unstructured mesh and the Crank-Nicolson procedure for the multi-term Riesz space and time fractional wave equation of distributed-order. Arianfar et al. [35] established the shifted Legendre collocation approach for the numerical solution of a class of nonlinear fractional boundary value problems of distributed-order with singularity in coefficients where the composite midpoint quadrature rule has been used to discretise the distributed-order time-fractional derivative.
The approximation methods based on spectral techniques and orthogonal functions have become a rigorous tool for solving various fractional dynamical systems [36,37]. These nonlocal methods are efficient numerical schemes for discretising the nonlocal fractional differential operators [38,39]. The main benefit of the spectral techniques is their ability to generate accurate approximate results with a few degrees of freedom. The orthogonal functions are often classified into three categories based on their structure: sine-cosine functions; piecewise orthogonal constant functions (block-pulse, hybrid, etc.); and orthogonal polynomials (Gegenbauer, Legendre, Laguerre, etc.). Using orthogonal polynomials as a basis function for the spectral methods gives preferable approximate results [40][41][42]. The main advantage of applying orthogonal functions is that they can be used to convert problems into systems of algebraic equations, substantially simplifying them. As far as we know, the spectral solutions for distributed-order FDEs are rarely discussed. Dehghan and Abbaszadeh [43] proposed the Legendre spectral element technique for simulating the neutral delay distributed-order fractional diffusion-wave equation with damping. Pourbabaee and Saadatmandi [44] presented a spectral method based on operational matrix of shifted Legendre for distributed-order time-FDEs and diffusion equations. Zaky and Machado [45] developed spectral schemes for discretising the multi-dimensional distributed-order time-fractional diffusion equations. Zhang et al. [46] presented a Crank-Nicolson alternating direction implicit (ADI) Galerkin-Legendre spectral technique for the 2D space-Riesz advection-diffusion equation of distributed-order. Zhang et al. [47] proposed a spectral scheme based on ADI Legendre-Laguerre for the 2D-DOTFDWE on a semi-infinite domain. Moreover, the definite integrals of the distributed-order fractional derivatives in aforementioned articles have been discretised by using Gauss quadrature formula or composite Simpson formula or composite trapezoid formula. The Gauss quadrature formula is more stable than the other mentioned formulas and has a higher computational accuracy of order , where denotes the used polynomials number and is a suitable constant.
In the current work, we consider the following 1D and 2D DOTFTE with spatial variable coefficients: (1.4) It is worth to point out that if 2 0 and the density function 1 , where 1 is the Dirac delta function, then (1.1) reduces to the distributed-order time-fractional diffusion-wave equation with damping (DOTFDWED). If 1 2 0, then (1.1) reduces to the DOTFDWE. In this paper, we present two efficient spectral algorithms for solving 1D and 2D distributed-order time-fractional telegraph models. We construct a novel SG operational matrix for approximating the fractional derivatives of distributed-order included in these models. Also, a new operational matrix of the multiplication of space vectors has been built to have the ability for completely applying the tau method in the 2D case. The presented spectral algorithms convert the given problem into an algebraic equations system, which is easy to solve. Also, we discuss the convergence and error bound analysis of the presented methods. Some numerical applications are presented to demonstrate the the effectiveness of the presented method.
The paper is arranged as follows: Section 2 introduces mathematical preliminaries of SG polynomials and its using for function approximation. In Section 3, the operational matrix of distributed-order fractional differentiation of the SG polynomials is derived. Section 4 develops two spectral algorithms for solving the 1D and 2D distributed-order models (4.1), (4.2) and (4.13), (4.14). In Section 5, the convergence and error bound analysis of the presented techniques are demonstrated. In Section 6, miscellaneous test examples on the multi-dimensional model (1.1), (1.2) are presented. Section 7 ends the paper with a conclusion.

SG polynomials
Gegenbauer polynomials with degree and parameter 1 2 arises as the eigensolution on the interval 1 1 for the following singular Sturm-Liouville problem The Gegenbauer polynomials are standardised in [48] such that 1 1 for all non-negative integer . Moreover, the following recurrence relation defines the standardised Gegenbauer polynomials : The Caputo fractional derivative of the SG polynomials has the following form [37] 1 1 0 . (2.5)

Functions approximation and fractional derivative operational matrix
The SG polynomials can be used to approximate a function as follows: where denotes the transpose, the SG expansion coefficients can be calculated from (2.7) and The SG polynomials can be used to approximate a function as follows: here A is an unknown matrix whose entries can be calculated from The Kronecker product of a matrix of order and a matrix of order is a matrix of order which has the following block form 11  where the symbol is the Kronecker tensor product. Also, the Kronecker product satisfies the following properties for the matrices and .
(2.11) Accordingly, the function 2 can be written with the following structure: where the coefficients matrix A is written in the block form as follows: The following theorem gives the matrix representation of the fractional/integer order derivatives for the SG vector .

The proposed methodologies
In this section, the procedures solution of the presented strategies for solving the multi-dimensional DOTFTE models are discussed in details.

1D-case
Consider the following 1D DOTFTE: We begin the approach by approximate the known functions for 1 2 3 and for 0 1 as follows: where B 1 B 2 and B 3 are 1 column matrices whose entries can be calculated similarly as (2.7). Also, F F 0 F 1 G 0 and G 1 are 1 1 known matrices whose entries 0 1 0 and 1 can be calculated similarly as (2.10).
According to [49], consider that B is the vector form approximation of the smooth and continuous function with B 0 1 ... , then the following relation holds: where W is an 1 square matrix defined as follows: Now, by employing (2.9), (2.15), (3.5) and (4.4), the following approximations are attained 2 1 D 1 2 A (4.6) (4.9) Substituting (4.3), (4.6), (4.7), (4.8) and (4.9) into (4.1), we get D 1 2 A D 0 1 AW 1 AW 2 AD 2 W 3 F 0 (4.10) by applying the orthogonality relation (2.3) and by using the typical tau method, we can extract an 1 1 system of linear algebraic equations in terms of unknown coefficients as follows: 1 1 (4.12) By solving the algebraic system of (4.11) and (4.12), the unknown coefficients are computed and hence the approximate solution is determined. The main steps of the suggested methodology for solving the 1D DOTFTE are implemented in the following algorithm.

2D-case
Consider the following 2D DOTFTE: where is the finite rectangular region 0 0 in 2 , and 0 represents the time interval.
Firstly, the functions for 1 2 3 4 can be approximated by the SG polynomials as follows: where F F 0 F 1 G 0 G 1 H 0 and H 1 are known block form matrices written likewise (2.12) and their elements can be calculated likewise (2.14).
The following theorem is substantial to remove the obstacle that facing the use of the tau method for solving the variable coefficients 2D DOTFTE model (4.13):    where By applying the orthogonality relation (2.3) for (4.30) and by using the typical tau method, we can extract an 1 1 1 algebraic system of linear equations in terms of unknown coefficients as follows: By solving the algebraic system of (4.32) and (4.33) for the unknowns coefficients. As a result, the approximate solution is determined easily. The main steps of the suggested methodology for solving the 2D DOTFTE are implemented in the following algorithm.

Convergence and error estimation
This section focuses on the convergence and error bound analysis of the proposed numerical technique in the weighted 2 -norm. Algorithm 2 2D DOTFTE.

1D-case
Consider the polynomials space of -degree span 0 . Let 0 0 . We define the weighted norm and inner product as follows: In this subsection, we will consider the following notations before introducing the main theorems:

2D-case
Consider the polynomials space of -degree span 0 . Let 0 0 0 . We define the weighted norm and inner product as follows: where . The generalised Taylor's expansion [50] of finite terms can be defined as follows: In this subsection, we will consider the following notations before introducing the main theorems: (5.14) The triple integration of (5.14) can be estimated as follows: Now, by substituting from (5.15) in (5.14), then the square root of the resulting equation completes the proof.
The following corollaries 5.4, 5.5, 5.6 and 5.7 can be proved by using the same procedure as in Theorem 5.3. The weight function of example 2 has the form where the parameter 0 physically represents the relaxation time. Figure 2 demonstrates the effect of on the behaviour of the wave. It is noted that, with gradually increase of , the wave becomes higher. The graphical representation of our obtained results is in good agreement with the known results.  1.
The MAEs at various time levels of cases 1 and 2 of example 4 are listed in Table 3  with 4, and 0.5. Figure 3 shows the approximate solution and  absolute error surfaces of case 3 at time level 0.5 with 4, and 0.5. In Table 4

Conclusion
In this paper, the DOTFTEs in multi-dimensions were considered where the distributed-order over the intervals 1 2 and 0 1 are existed simultaneously. The certain choices of the density function and the coefficients can give the so-called DOTFDWED and DOTFDWE which are relevant to sub/super-diffusive processes; therefore, the existence of the accurate solutions of these equations is required. Fig. 3 The approximate solution and absolute error with =4, = =10 at 0.5 for case 3 of example 4 The proposed method 1 Method in [32] 0.49, The main objective of this paper was to offer improved and accurate spectral tau algorithms for numerical solutions of the aforementioned equations in multidimensions with constant and variable space coefficients based on a novel operational matrix of distributed-order fractional derivative of SG polynomials. A new operational matrix of the multiplication of space vectors has been built to have the ability in applying the tau method in the 2D DOTFTE. The convergence analysis of the presented algorithms has been discussed. The numerical solutions obtained by the proposed spectral tau algorithms have high accuracy with a small number of basis functions for smooth and non-smooth solutions as shown from the test examples. For non-smooth solutions, we used a smooth transformation by using SG polynomials to avoid the convergence deteriorated.