Hexic-Chebyshev Collocation Method for Solving Distributed-Order Time-Space Fractional Diffusion Equations

Abstract

This paper presents a spectral method to solve nonlinear distributed-order diffusion equations with both time-distributed-order and two-sided space-fractional terms. These are highly challenging to solve analytically due to the interplay between nonlinearity and the fractional distributed-order nature of the time and space derivatives. For this purpose, Hexic-kind Chebyshev polynomials (HCPs) are used as the backbone of the method to transform the primary problem into a set of nonlinear algebraic equations, which can be efficiently solved using numerical solvers, such as the Newton–Raphson method. The primary reason of choosing HCPs is due to their remarkable recurrence relations, facilitating their efficient computation and manipulation in mathematical analyses. A comprehensive convergence analysis was conducted to validate the robustness of the proposed method, with an error bound derived to provide theoretical guarantees for the solution’s accuracy. The method’s effectiveness is further demonstrated through two test examples, where the numerical results are compared with existing solutions, confirming the approach’s accuracy and reliability.

1. Introduction

Since the turn of the century, fractional calculus has garnered significant interest among researchers across diverse fields of science and engineering, including physics, biology, finance, and fluid mechanics [1,2,3,4,5,6]. Fractional-order operators, characterized by their nonlocality, retain memory and hereditary properties inherent in natural phenomena. These distinctive features enable scientists to provide more accurate representations of real-life scenarios [7,8,9,10,11,12,13]. Applications in anomalous diffusion [14,15,16,17,18,19], viscoelasticity [20,21], control problems [22], and biological systems [23] have further driven the development of advanced numerical techniques that address the non-locality inherent in fractional operators. However, constant-order fractional operators present limitations in the mathematical modeling of multi-physics problems. Moreover, recent studies show that fixed-order fractional differential equations are not capable of adequately characterizing some types of complex natural phenomena, such as complex diffusion encountered in highly heterogeneous fractured or disordered porous media [24,25,26]. Thus, more robust mathematical tools are required to model increasingly complex systems effectively. The scope of fractional calculus extends beyond constant-fractional-order operators, and this flexibility allows for its application in systems with varying dynamic behavior.

The concept of variable-order differential operators was introduced to address the limitations of constant-order operators, with several definitions proposed [27]. This shift from fixed-order to variable-order operators opens new possibilities for describing complex systems mathematically. Another class of fractional differential equations, distributed-order fractional differential equations, differs fundamentally from variable-order equations in structure and complexity. Distributed-order equations do not have a fixed fractional order or rely on a single variable. Instead, the fractional order is distributed over a range of values, requiring integration across this range [15,28,29,30,31,32,33,34]. This approach captures a wider range of dynamic behaviors across different scales, making it especially useful for memory effects processes and complex temporal or spatial variations. In contrast, variable-order equations allow the fractional order to change only as a function of time or space, offering less flexibility for systems with evolving properties. Therefore, distributed-order models are better suited for scenarios where underlying processes exhibit multifaceted, non-uniform evolution over time or space, providing a more comprehensive framework for modeling real-world phenomena.

Distributed-order fractional models have been instrumental in representing materials with time-varying properties, a characteristic that traditional calculus fails to capture. The versatility of distributed-order operators makes them particularly effective for capturing the evolving nature of materials, processes, and systems, further fueling the growing interest in this area of study. Anomalous diffusion, characterized by multiple scaling exponents in the mean squared displacement, is a process effectively modeled by distributed-order fractional differential equations [30]. The ability of fractional models to couple memory effects with non-integer differentiation orders makes them particularly effective in simulating processes with long-range temporal dependencies. This is evident in applications such as heat conduction in heterogeneous materials, diffusion in porous media, and other complex systems where the standard integer-order models fall short of capturing the underlying dynamics.

Given the wide-ranging applications of distributed-order fractional differential equations in science and engineering, research in this area has grown substantially, strongly emphasizing the development of numerical methodologies for solving these equations. Over the past two decades, various models based on distributed-order fractional differential equations have been proposed [21,35,36,37,38,39,40]. Analytical and computational investigations of these models have deepened our understanding of their applicability across numerous fields [41,42]. Several efficient methods have emerged, including Petrov–Galerkin and spectral collocation methods [43], the Legendre collocation method [44], the block-centered finite difference method [45], the composite collocation method [46], the Legendre operational matrix technique [47], the shifted Legendre-series method coupled with the composite midpoint quadrature rule [48], finite difference schemes [49], space-time finite element methods [50], POD-based Crank–Nicolson/fourth-order ADI finite difference schemes [51], and the split-step method [52]. These approaches have continued to enhance the accuracy and computational feasibility of fractional differential models in real-world applications.

Classical Chebyshev polynomials are known for their orthogonality and ability to approximate functions with high accuracy using minimal terms, making them powerful tools in spectral methods for solving differential equations [53,54,55,56]. These polynomials are defined over the interval $[-1,1]$ and possess remarkable properties, such as their recurrence relations and orthogonality concerning a weight function. This property ensures that the polynomials can represent complex solutions with fewer basis functions, improving computational efficiency [57,58,59]. Hence, they are ideal for numerical methods, such as spectral collocation, where differential equations are transformed into algebraic equations to be solved at the collocation points [57].

Hexic-kind Chebyshev polynomials (HCPs) represent an extension of traditional Chebyshev polynomials and possess unique features that make them especially well suited for addressing fractional differential equations, particularly when complex boundary conditions and fractional operators are involved [60]. These polynomials demonstrate several important mathematical properties, including recurrence relations, orthogonality conditions, and generating functions, which enable the conversion of fractional differential equations into systems of nonlinear algebraic equations. When compared to conventional Chebyshev and Legendre polynomials, HCPs provide superior stability characteristics for fractional differential operators and deliver improved numerical conditioning, particularly in applications involving higher-order derivatives. The specialized recurrence relations and orthogonal structure of HCPs lead to decreased computational costs during matrix formation processes.

This study presents a Hexic-kind Chebyshev collocation (HCC) approach for solving sophisticated nonlinear time-distributed-order and two-sided space-fractional diffusion equations, building upon the spectral methodology previously established by the same authors using HCPs for variable-order fractional partial integrodifferential equations [21]. Given the inherent mathematical complexity of these equations, which renders conventional analytical methods inadequate, the development of this advanced numerical framework becomes essential. The methodology initiates by employing the Gauss quadrature rule to approximate the time-distributed-order component, subsequently utilizing HCPs within a collocation framework to convert the original problem into a tractable system of nonlinear algebraic equations. The distinctive contribution of this research lies in incorporating HCPs within an innovative spectral collocation methodology to tackle nonlinear distributed-order fractional diffusion equations that feature both time-distributed and two-sided space-fractional operators, differentiating it from conventional spectral approaches that rely on traditional Chebyshev or Legendre polynomials. The innovation centers on the effective development of differentiation matrices for two-sided Riemann–Liouville operators utilizing HCPs, combined with the integration of the Gauss quadrature for distributed-order temporal derivatives, which facilitates exponential convergence and delivers highly precise and stable solutions using a relatively small number of basis functions. This methodology proves especially advantageous for addressing problems characterized by memory-dependent behavior, irregular source terms, and complex boundary specifications. The developed method demonstrates exceptional computational efficiency and outstanding accuracy, with numerical results indicating that precise solutions can be achieved using minimal basis functions.

The manuscript is structured as follows: Section 2 establishes the fundamental definitions and essential concepts that underpin the theoretical framework, including a comprehensive introduction to HCPs, their mathematical notation, and characteristic properties. Section 3 formulates the nonlinear time-distributed-order and two-sided space-fractional diffusion equation that constitutes the primary subject of investigation. Section 4 develops the convergence analysis for the proposed methodology, establishing its theoretical validity and stability. Section 5 presents two representative numerical examples that demonstrate the accuracy and computational efficiency of the developed approach. Section 6 provides concluding observations regarding the research findings and outlines promising avenues for future investigations in this field.

2. Preliminaries

This section presents an overview of fractional derivatives and fundamental properties of Hexic-kind Chebyshev polynomials (HCPs), which are pivotal in formulating and applying the proposed numerical method. Following this, a concise review of the Gaussian–Legendre quadrature is provided.

2.1. Fractional Calculus

Definition 1

([1]). For $\beta \in (1,2)$, the β-order left- and right-sided RL derivatives are defined as

(1a)$\begin{array}{ccc}\hfill {}_{a }{}^{RL}\mathcal{D}_x^{\beta }u(x,t)& =& {\displaystyle \frac{1}{\Gamma (2-\beta )}}{\displaystyle \frac{{\partial }^2}{\partial x^2}}{\int }_a^x{(x-s)}^{1-\beta }u(s,t)\mathrm{d}s,\hfill \end{array}$

(1b)$\begin{array}{ccc}\hfill {}_{x }{}^{RL}\mathcal{D}_b^{\beta }u(x,t)& =& {\displaystyle \frac{1}{\Gamma (2-\beta )}}{\displaystyle \frac{{\partial }^2}{\partial x^2}}{\int }_x^b{(s-x)}^{1-\beta }u(s,t)\mathrm{d}s.\hfill \end{array}$

Lemma 1

([1]). The left-sided RL fractional derivative of the power function $x^r$ is given by

(2)${}_{0 }{}^{RL}\mathcal{D}_x^{\beta }{x}^r={\displaystyle \frac{\Gamma (r+1)}{\Gamma (r+1-\beta )}}{x}^{r-\beta }.$

Definition 2

([1]). For $\alpha \in (0,1)$, the α-order Caputo fractional derivative is defined as

(3)${}_0{}^C\mathcal{D}_t^{\alpha }u(x,t)={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_0^t{u}_s(x,s){(t-s)}^{-\alpha }\mathrm{d}s.$

Definition 3

([16]). The distributed-order fractional derivative ${\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}[·]$ is defined as

(4)${\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}u(x,t)={\int }_0^1\vartheta \left(\alpha \right){}_0{}^C\mathcal{D}_t^{\alpha }u(x,t)\mathrm{d}\alpha ,$

where $\alpha \in [0,1]$ and the weight function $\vartheta \left(\alpha \right)$ is a non-negative continuous function satisfying $\vartheta \left(\alpha \right)\ne 0$ and ${\int }_0^1\vartheta \left(\alpha \right)\mathrm{d}\alpha <\infty$.

2.2. Hexic-Kind Chebyshev Polynomials

Definition 4

([61]). The Hexic-kind Chebyshev polynomials are defined over the interval $[-1,1]$ and are orthogonal regarding a specific weight function. The recurrence relation characterizes them

(5)$\begin{array}{ccc}\hfill {\widehat{\psi }}_0\left(t\right)& =& 1,{\widehat{\psi }}_1\left(t\right)=t,\hfill \\ \hfill {\widehat{\psi }}_{j+1}\left(t\right)& =& t{\widehat{\psi }}_j\left(t\right)+{\overline{\Theta }}_j{\widehat{\psi }}_{j-1}\left(t\right),j=2,3,\dots ,\hfill \end{array}$

where the coefficients ${\overline{\Theta }}_j$ are given by

(6)${\overline{\Theta }}_j=-{\displaystyle \frac{(1+j+{(-1)}^{1+j})(j+2+{(-1)}^{1+j})}{4(1+j)(2+j)}}.$

Subsequently, the shifted HCPs on the interval $[0,t_{max}]$ are defined as

(7) ${\psi }_j\left(t\right)={\widehat{\psi }}_j\left({\displaystyle \frac{2}{t_{max}}}t-1\right),j=0,1,2,\dots .$

Additionally, these polynomials have an explicit analytic form as derived in [62]:

(8) ${\psi }_j\left(t\right)=\sum _{l=0}^j{b}_{l,j}{\left({\displaystyle \frac{t}{t_{max}}}\right)}^l,$

where the coefficients $b_{l,j}$ are given by:

(9) $b_{l,j}=\left\{\begin{array}{cc}{\displaystyle \frac{2^{2l-j}}{(2l+1)!}}{\sum }_{i=\lfloor {\displaystyle \frac{l+1}{2}}\rfloor }^{{\displaystyle \frac{j}{2}}}{\displaystyle \frac{{(-1)}^{\frac{j}{2}+1+i}(1+2i+l)!}{(2i-l)!}},\hfill & jiseven,\hfill \\ {\displaystyle \frac{2^{2l-j+1}}{(2l+1)!(j+1)}}{\sum }_{i=\lfloor {\displaystyle \frac{l}{2}}\rfloor }^{{\displaystyle \frac{j-1}{2}}}{\displaystyle \frac{{(-1)}^{\frac{j+1}{2}+l+i}(1+i)(l+2i+2)!}{(2i-l+1)!}},\hfill & jisodd.\hfill \end{array}\right.$

Subsequent analyses will utilize the key properties of HCPs, including their recurrence relations, orthogonality conditions, and differentiation formulas, to develop the numerical framework. In addition, their Gauss-type quadrature integration formula will be applied to approximate distributed-order terms, providing an efficient means of handling time-varying and space-fractional components in the diffusion equation.

Definition 5.

The collection of basis functions ${\left\{{\psi }_j\left(t\right)\right\}}_{j\in \mathbb{N}\cup \left\{0\right\}}$ constitutes a set of orthogonal functions associated with the weight function

(10) $\varpi \left(t\right)={(2{\displaystyle \frac{t}{t_{max}}}-1)}^2\sqrt{{\displaystyle \frac{1}{t_{max}}}\left(t-{\displaystyle \frac{t^2}{t_{max}}}\right)},$

defined on the interval $[0,t_{max}]$.

Thus, any function $u\left(t\right)\in L_{\varpi }^2[0,t_{max}]$ can be approximated by the shifted HCPs as

(11)$u\left(t\right)\simeq u_m\left(t\right)=\sum _{j=0}^m{\mathsf{u}}_j{\psi }_j\left(t\right)\triangleq {\mathbf{u}}^{\top }{\Psi }_m\left(t\right),$

where

(12)$\begin{array}{ccc}\hfill \mathbf{u}& :=& {[{\mathsf{u}}_0,\dots ,{\mathsf{u}}_j,\dots ,{\mathsf{u}}_m]}^{\top },\hfill \\ \hfill {\Psi }_m\left(t\right)& :=& {[{\psi }_0\left(t\right),\dots ,{\psi }_j\left(t\right),\dots ,{\psi }_m\left(t\right)]}^{\top },\hfill \end{array}$

in which the coefficients ${\mathsf{u}}_j$ are obtained as

(13)${\mathsf{u}}_j={\displaystyle \frac{1}{c_j}}{\int }_0^{t_{max}}u\left(t\right)\varpi \left(t\right){\psi }_j\left(t\right)\mathrm{d}t,$

with

(14)$\begin{array}{c}\hfill c_i=\left\{\begin{array}{cc}{\displaystyle \frac{\pi }{2^{2i+5}}},& ieven,\hfill \\ {\displaystyle \frac{\pi (i+3)}{2^{2i+5}(i+1)}},\hfill & iodd.\hfill \end{array}\right.\end{array}$

Theorem 1.

Suppose ${\Psi }_m\left(t\right)$ represents the vector of shifted HCPs as defined in Equation (12). Then, its fractional derivative is given by

(15) ${\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}{\Psi }_m\left(t\right)={\Phi }^{\vartheta }\left(t\right),$

where ${\Phi }^{\vartheta }\left(t\right)$ denotes the distributed-order fractional derivative of the vector ${\Psi }_m\left(t\right)$ and is characterized as

(16) ${\Phi }^{\vartheta }\left(t\right)={\left[0,\sum _{r=1}^1{\varphi }_{r,1}^{\vartheta }\left(t\right),\dots ,\sum _{r=1}^j{\varphi }_{r,j}^{\vartheta }\left(t\right),\dots ,\sum _{r=1}^m{\varphi }_{r,m}^{\vartheta }\left(t\right)\right]}^{\top },$

where

(17) ${\varphi }_{r,j}^{\vartheta }\left(t\right)=\sum _{i=0}^{\overline{M}}{\varpi }_i^{\overline{M}}\vartheta \left({\varsigma }_i\right){\displaystyle \frac{\Gamma (r+1)}{t_{max}^r\Gamma (r+1-{\varsigma }_i)}}{b}_{r,j}{t}^{r-{\varsigma }_i},$

and ${\left\{{\varpi }_i^{\overline{M}}\right\}}_{i=0}^{\overline{M}}$ represent the Legendre–Gauss weights, while the Legendre–Gauss nodes ${\left\{{\varsigma }_i\right\}}_{i=0}^{\overline{M}}$ are the zeros of ${\overline{P}}_{\overline{M}+1}\left(x\right)$.

Proof. 

Utilizing Definition 3 and the shifted HCPs, for $j=0,\dots ,m$, the distributed-order fractional derivative ${\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}{\psi }_j\left(t\right)$ is expressed as

(18)$\begin{array}{c}\hfill {\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}{\psi }_j\left(t\right)={\int }_0^1\vartheta \left(\alpha \right){}_0{}^C\mathcal{D}_t^{\alpha }{\psi }_j\left(t\right)\mathrm{d}\alpha ,\end{array}$

and by applying Lemma 2, the approximation is obtained as

(19)$\begin{array}{c}\hfill {\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}{\psi }_j\left(t\right)\simeq \sum _{i=0}^{\overline{M}}{\varpi }_i^{\overline{M}}\vartheta \left({\varsigma }_i\right){}_0{}^C\mathcal{D}_t^{\varsigma i}{\psi }_j\left(t\right),\end{array}$

where the Legendre–Gauss nodes ${\left\{{\varsigma }_i\right\}}_{i=0}^{\overline{M}}$ represent the zeros of ${\overline{P}}_{\overline{M}+1}\left(x\right)$ on $[0,1]$, and the associated Legendre–Gauss weights ${\left\{{\varpi }_i^{\overline{M}}\right\}}_{i=0}^{\overline{M}}$ are given by

(20)${\varpi }_i^{\overline{M}}={\displaystyle \frac{2}{(1-{\varsigma }_i^2){\left({{\overline{P}}^{\prime }}_{\overline{M}+1}\left({\varsigma }_i\right)\right)}^2}}.$

Considering Equation (8), the following is derived:

(21)${}_0{}^C\mathcal{D}_t^{\alpha }{\psi }_0\left(t\right)=b_{0,0}{}_0{}^C\mathcal{D}_t^{\alpha }\left(1\right)=0,$

which leads to ${\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}{\psi }_0\left(t\right)=0$, as shown by Equations (19) and (21). Furthermore, it is known that

(22)${}_0{}^C\mathcal{D}_t^{\alpha }{t}^r={\displaystyle \frac{\Gamma (r+1)}{\Gamma (r+1-\alpha )}}{t}^{r-\alpha },r\in \mathbb{N},$

thus, for $j=1,\dots ,m$, it can be deduced that

(23)$\begin{array}{c}\hfill {}_0{}^C\mathcal{D}_t^{\alpha }{\psi }_j\left(t\right)=\sum _{r=0}^j{b}_{r,j}{}_0{}^C\mathcal{D}_t^{\alpha }{\left({\displaystyle \frac{t}{t_{max}}}\right)}^r=\sum _{r=1}^j{\displaystyle \frac{\Gamma (r+1)}{t_{max}^r\Gamma (r+1-\alpha )}}{b}_{r,j}{t}^{r-\alpha },\end{array}$

and subsequently, from Equations (19) and (23), the following is derived:

(24)$\begin{array}{c}\hfill {\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}{\psi }_j\left(t\right)\simeq \sum _{r=1}^j{\varphi }_{r,j}^{\vartheta }\left(t\right),\end{array}$

where

(25)${\varphi }_{r,j}^{\vartheta }\left(t\right)=\sum _{i=0}^{\overline{M}}{\varpi }_i^{\overline{M}}\vartheta \left({\varsigma }_i\right){\displaystyle \frac{\Gamma (r+1)}{t_{max}^r\Gamma (r+1-{\varsigma }_i)}}{b}_{r,j}{t}^{r-{\varsigma }_i}.$

□

Definition 6.

Legendre polynomials are commonly utilized in numerical methods because of their advantageous properties, including orthogonal behavior and efficient convergence rates. They are defined on the interval $[-1,1]$ by the recurrence formula

(26) $\begin{array}{ccc}\hfill & & P_0\left(x\right)=1,P_1\left(x\right)=x,\hfill \\ & & P_{i+1}\left(x\right)={\displaystyle \frac{2i+1}{i+1}}x{P}_i\left(x\right)-{\displaystyle \frac{i}{i+1}}{P}_{i-1}\left(x\right),i=2,3,\dots .\hfill \end{array}$

In addition, the shifted Legendre polynomials defined on the interval $[0,\tau ]$ are expressed as

${\overline{P}}_i\left(x\right)=P_i\left({\displaystyle \frac{2}{\tau }}x-1\right).$

Lemma 2

([63]). Suppose $u\left(x\right)\in C^{2\overline{M}}[0,\tau ]$. The Legendre–Gauss quadrature formula on the interval $[0,\tau ]$ is given by

(27)${\int }_0^{\tau }u\left(x\right)\mathrm{d}x=\sum _{i=0}^{\overline{M}}{\varpi }_i^{\overline{M}}u\left({\varsigma }_i\right)+r_{LG},$

where the Legendre–Gauss nodes ${\left\{{\varsigma }_i\right\}}_{i=0}^{\overline{M}}$ are the zeros of ${\overline{P}}_{\overline{M}+1}\left(x\right)$, the Legendre–Gauss weights ${\left\{{\varpi }_i^{\overline{M}}\right\}}_{i=0}^{\overline{M}}$ are defined as

(28)${\varpi }_i^{\overline{M}}={\displaystyle \frac{2}{(1-{\varsigma }_i^2){\left({\overline{P}}_{\overline{M}+1}^{\prime }\left({\varsigma }_i\right)\right)}^2}},$

and the remainder term is given by

(29)$r_{LG}={\displaystyle \frac{{\tau }^{1+2\overline{M}}{(\overline{M}!)}^4}{(1+2\overline{M}){[\left(2\overline{M}\right)!]}^3}}{u}^{\left(2\overline{M}\right)}\left(\xi \right),\xi \in [0,\tau ].$

3. Hexic Chebyshev Collocation for Distributed-Order Fractional Diffusion

Consider the following nonlinear fractional diffusion equation with time-space distributed orders, which was previously explored by a linear source term in [64]:

(30)${\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}u(x,t)={\Theta }_l(x,t){}_{0 }{}^{RL}\mathcal{D}_x^{\beta }u(x,t)+{\Theta }_r(x,t){}_{x }{}^{RL}\mathcal{D}_{x_{\max }}^{\beta }u(x,t)+u(x,t,u(x,t)),$

within the domain $(x,t)\in \Omega \times \mathcal{I}$, accompanied by the initial and Dirichlet boundary conditions

(31a)$\begin{array}{ccc}\hfill u(x,0)& =& \rho \left(x\right),\hfill \end{array}$

(31b)$\begin{array}{ccc}\hfill u(0,t)& =& {\phi }_0\left(t\right),u(x_{max},t)={\phi }_{x_{max}}\left(t\right),\hfill \end{array}$

where the interval $\mathcal{I}:=[0,t_{max}]$ and the domain $\Omega :=[0,x_{max}]$ are specified. The functions ${\Theta }_l(x,t)$ and ${\Theta }_r(x,t)$ are known, while $u(x,t,u)$, a nonlinear source term, belongs to $C^1(\Omega \times \mathcal{I}\times \mathbb{R})$ and satisfies Lipschitz continuity with respect to u. The functions $\rho \left(x\right)$, ${\phi }_0\left(t\right)$, and ${\phi }_{x_{max}}\left(t\right)$ are also given. The operator ${\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}[·]$ denotes the time-distributed-order derivative with a non-negative continuous weight function $\vartheta \left(\alpha \right)$ over $\alpha \in [0,1]$. Additionally, ${}_{0 }{}^{RL}\mathcal{D}_x^{\beta }[·]$ and ${}_{x }{}^{RL}\mathcal{D}_{x_{\max }}^{\beta }[·]$ represent the left- and right-sided Riemann–Liouville (RL) derivatives of order $\beta \in (1,2)$.

Due to the complexity of the problem, obtaining exact solutions is generally impractical, making it necessary to develop numerical methods for approximating solutions. Despite this difficulty, few numerical schemes have been proposed for such distributed-order problems. In [64], a fast implicit difference method was introduced to solve Equations (30)–(31b) with a linear source term. Similarly, an implicit numerical method was proposed in [65] to address two-sided space-fractional and time-distributed-order advection-dispersion equations. The proposed HCP method is subsequently outlined, detailing the construction of the nonlinear algebraic equations, which can be solved using a nonlinear solver. To achieve this, the following two theorems will first be presented and proved to construct a differentiation matrix for approximating the fractional derivatives of a function using HCPs.

Theorem 2.

Consider the shifted HCPs vector ${\Psi }_n\left(x\right)$ defined as

(32) ${\Psi }_n\left(x\right):={[{\psi }_0\left(x\right),\dots ,{\psi }_i\left(x\right),\dots ,{\psi }_n\left(x\right)]}^{\top },$

where ${\psi }_i\left(x\right)={\widehat{\psi }}_i\left({\displaystyle \frac{2}{x_{max}}}x-1\right)$. Then, the left-sided RL fractional derivative of ${\Psi }_n\left(x\right)$ is denoted as ${\mathsf{D}}_{\mathrm{L}}^{RL}\left(x\right)$ and defined as

(33) $\begin{array}{ccc}\hfill & & {\mathsf{D}}_{\mathrm{L}}^{RL}\left(x\right)=\hfill \\ & & {\left[\sum _{r=0}^0{\Delta }_{\mathsf{l}\mathsf{e}\mathsf{f}\mathsf{t}}^{r,0,\beta }\left(x\right),\sum _{r=0}^1{\Delta }_{\mathsf{l}\mathsf{e}\mathsf{f}\mathsf{t}}^{r,1,\beta }\left(x\right),\dots ,\sum _{r=0}^i{\Delta }_{\mathsf{l}\mathsf{e}\mathsf{f}\mathsf{t}}^{r,i,\beta }\left(x\right),\dots ,\sum _{r=0}^n{\Delta }_{\mathsf{l}\mathsf{e}\mathsf{f}\mathsf{t}}^{r,n,\beta }\left(x\right)\right]}^{\top },\hfill \end{array}$

where

(34) ${\Delta }_{\mathsf{l}\mathsf{e}\mathsf{f}\mathsf{t}}^{r,i,\beta }\left(x\right)={\displaystyle \frac{\Gamma (r+1)}{x_{max}^r\Gamma (r+1-\beta )}}{b}_{r,i}{x}^{r-\beta }.$

Proof. 

Employing Lemma 1 and the explicit form of the shifted HCPs, for $i=0,\dots ,n$, the following is obtained

(35)$\begin{array}{c}\hfill {}_{0 }{}^{RL}\mathcal{D}_x^{\beta }{\psi }_i\left(x\right)=\sum _{r=0}^i{b}_{r,i}{}_{0 }{}^{RL}\mathcal{D}_x^{\beta }{\left({\displaystyle \frac{x}{x_{max}}}\right)}^r=\sum _{r=0}^i{\Delta }_{\mathsf{left}}^{r,i,\beta }\left(x\right),\end{array}$

where

(36)${\Delta }_{\mathsf{left}}^{r,i,\beta }\left(x\right)={\displaystyle \frac{\Gamma (r+1)}{x_{max}^r\Gamma (r+1-\beta )}}{b}_{r,i}{x}^{r-\beta }.$

□

Theorem 3.

Let ${\Psi }_n\left(x\right)$ represent the shifted HCPs vector as defined in (32). Then, the right-sided RL fractional derivative of ${\Psi }_n\left(x\right)$, denoted as ${\mathsf{D}}_{\mathrm{R}}^{RL}\left(x\right)$, is given by:

(37) ${\mathsf{D}}_{\mathrm{R}}^{RL}\left(x\right)={\left[\sum _{r=0}^0{\Delta }_{\mathsf{r}\mathsf{i}\mathsf{g}\mathsf{h}\mathsf{t}}^{r,0,\beta }\left(x\right),\sum _{r=0}^1{\Delta }_{\mathsf{r}\mathsf{i}\mathsf{g}\mathsf{h}\mathsf{t}}^{r,1,\beta }\left(x\right),\dots ,\sum _{r=0}^i{\Delta }_{\mathsf{r}\mathsf{i}\mathsf{g}\mathsf{h}\mathsf{t}}^{r,i,\beta }\left(x\right),\dots ,\sum _{r=0}^n{\Delta }_{\mathsf{r}\mathsf{i}\mathsf{g}\mathsf{h}\mathsf{t}}^{r,n,\beta }\left(x\right)\right]}^{\top },$

where

(38) $\begin{array}{cc}\hfill & {\Delta }_{\mathsf{r}\mathsf{i}\mathsf{g}\mathsf{h}\mathsf{t}}^{r,i,\beta }\left(x\right)=\\ & {\sum }_{r=0}^i{\displaystyle \frac{b_{r,i}}{x_{max}^r\Gamma (2-\beta )}}{\displaystyle \frac{{\partial }^2}{\partial x^2}}\left({(x_{max}-x)}^{2-\beta }{\sum }_{l=0}^r{\displaystyle \frac{r!x_{max}^l\Gamma (l-\beta +2)}{l!\Gamma (r-\beta +3)}}{x}^{r-l}\right).\end{array}$

Proof. 

Employing the explicit form of the shifted HCPs for $i=0,\dots ,n$, the following is derived:

(39)$\begin{array}{ccc}\hfill {}_{x }{}^{RL}\mathcal{D}_{x_{\max }}^{\beta }{\psi }_i\left(x\right)& =& \sum _{r=0}^i{b}_{r,i}{}_{x }{}^{RL}\mathcal{D}_{x_{\max }}^{\beta }{\left({\displaystyle \frac{x}{x_{max}}}\right)}^r\hfill \\ \hfill & =& \sum _{r=0}^i{\displaystyle \frac{b_{r,i}}{x_{max}^r\Gamma (2-\beta )}}{\displaystyle \frac{{\partial }^2}{\partial x^2}}\left({\int }_x^{x_{max}}{\displaystyle \frac{s^r}{{(s-x)}^{\beta -1}}}\mathrm{d}s\right).\hfill \end{array}$

Next, the integral expression in Equation (39) is evaluated. Introducing the substitution $s=(x_{max}-x)\lambda +x$, where $\lambda \in [0,1]$, yields:

(40)$\begin{array}{c}\hfill {\int }_x^{x_{max}}{\displaystyle \frac{s^r}{{(s-x)}^{\beta -1}}}\mathrm{d}s={\int }_0^1{\left((x_{max}-x)\lambda \right)}^{1-\beta }{\left((x_{max}-x)\lambda +x\right)}^r(x_{max}-x)\mathrm{d}\lambda \\ \hfill ={(x_{max}-x)}^{2-\beta }{\int }_0^1{\lambda }^{1-\beta }{\left(x_{max}\lambda +(1-\lambda )x\right)}^r\mathrm{d}\lambda .\end{array}$

It follows that

(41)${\left(x_{max}\lambda +(1-\lambda )x\right)}^r=\sum _{l=0}^r{\displaystyle \frac{r!}{l!(r-l)!}}{\left(x_{max}\lambda \right)}^l{\left((1-\lambda )x\right)}^{r-l}.$

Using the identity in Equation (41), the integral becomes

(42)$\begin{array}{ccc}\hfill {\int }_x^{x_{max}}{\displaystyle \frac{s^r}{{(s-x)}^{\beta -1}}}\mathrm{d}s& =& {(x_{max}-x)}^{2-\beta }\sum _{l=0}^r{\displaystyle \frac{r!x_{max}^l}{l!(r-l)!}}{x}^{r-l}{\int }_0^1{\lambda }^{l-\beta +1}{(1-\lambda )}^{r-l}\mathrm{d}\lambda \hfill \\ \hfill & =& {(x_{max}-x)}^{2-\beta }\sum _{l=0}^r{\displaystyle \frac{r!x_{max}^l}{l!(r-l)!}}{x}^{r-l}\mathrm{B}(l-\beta +2,r-l+1)\hfill \\ \hfill & =& {(x_{max}-x)}^{2-\beta }\sum _{l=0}^r{\displaystyle \frac{r!x_{max}^l\Gamma (l-\beta +2)\Gamma (r-l+1)}{l!(r-l)!\Gamma (r-\beta +3)}}{x}^{r-l}\hfill \\ & =& {(x_{max}-x)}^{2-\beta }\sum _{l=0}^r{\displaystyle \frac{r!x_{max}^l\Gamma (l-\beta +2)}{l!\Gamma (r-\beta +3)}}{x}^{r-l},\hfill \end{array}$

where $\mathrm{B}(·,·)$ is the Beta function. Therefore, in light of Equation (39), the following is derived:

(43)$\begin{array}{c}\hfill {}_{x }{}^{RL}\mathcal{D}_{x_{\max }}^{\beta }{\psi }_i\left(x\right)={\Delta }_{\mathsf{right}}^{r,i,\beta }\left(x\right),\end{array}$

with

(44)$\begin{array}{cc}\hfill & {\Delta }_{\mathsf{right}}^{r,i,\beta }\left(x\right)=\\ & {\sum }_{r=0}^i{\displaystyle \frac{b_{r,i}}{x_{max}^r\Gamma (2-\beta )}}{\displaystyle \frac{{\partial }^2}{\partial x^2}}\left({(x_{max}-x)}^{2-\beta }{\sum }_{l=0}^r{\displaystyle \frac{r!x_{max}^l\Gamma (l-\beta +2)}{l!\Gamma (r-\beta +3)}}{x}^{r-l}\right).\end{array}$

□

An approximate solution for Equation (30) is

(45)$u(x,t)\approx u_{n,m}(x,t)=\sum _{i=0}^n\sum _{j=0}^m{\delta }_{i,j}{\psi }_i\left(x\right){\psi }_j\left(t\right)={\Psi }_n{\left(x\right)}^{\top }\mathbf{D}{\Psi }_m\left(t\right),$

where ${\Psi }_m\left(t\right)$ and ${\Psi }_n\left(x\right)$ are defined in (12) and (32), respectively. The coefficients matrix $\mathbf{D}$ given by

(46)$\mathbf{D}=\left(\begin{array}{ccc}{\delta }_{0,0}& \dots & {\delta }_{0,m}\\ ⋮& \ddots & ⋮\\ {\delta }_{n,0}& \dots & {\delta }_{n,m}\end{array}\right),$

is unknown and necessitates determination. By considering Equations (30) and (45), and employing Theorems 1–3, the following is derived

(47)$\begin{array}{c}\hfill \Phi (x,t)\triangleq {\Psi }_n{\left(x\right)}^{\top }\mathbf{D}{\Phi }^{\vartheta }\left(t\right)-{\Theta }_l(x,t){\mathsf{D}}_{\mathrm{L}}^{RL}{\left(x\right)}^{\top }\mathbf{D}{\Psi }_m\left(t\right)\\ \hfill -{\Theta }_r(x,t){\mathsf{D}}_{\mathrm{R}}^{RL}{\left(x\right)}^{\top }\mathbf{D}{\Psi }_m\left(t\right)-u(x,t,{\Psi }_n{\left(x\right)}^{\top }\mathbf{D}{\Psi }_m\left(t\right))\approx 0,\end{array}$

and from the initial and boundary conditions (31a) and (31b) alongside Equation (45), the following is obtained

(48a)$\begin{array}{ccc}\hfill {\Lambda }_1\left(t\right)& \triangleq & {\Psi }_n{\left(0\right)}^{\top }\mathbf{D}{\Psi }_m\left(t\right)-{\phi }_0\left(t\right)\approx 0,\hfill \end{array}$

(48b)$\begin{array}{ccc}\hfill {\Lambda }_2\left(t\right)& \triangleq & {\Psi }_n{\left(x_{max}\right)}^{\top }\mathbf{D}{\Psi }_m\left(t\right)-{\phi }_{x_{max}}\left(t\right)\approx 0,\hfill \end{array}$

(48c)$\begin{array}{ccc}\hfill \Pi \left(x\right)& \triangleq & {\Psi }_n{\left(x\right)}^{\top }\mathbf{D}{\Psi }_m\left(0\right)-\rho \left(x\right)\approx 0.\hfill \end{array}$

The discrete shifted Hexic-kind Chebyshev basis yields quadrature nodes as

(49a)$\begin{array}{c}\hfill x_0=0,x_n=x_{max},\forall i=1,\dots ,n-1,x_i\mathrm{are}\mathrm{roots}\mathrm{of}{\psi }_{n-1}\left(x\right),\end{array}$

(49b)$\begin{array}{c}\hfill \forall j=1,\dots ,m,t_j\mathrm{are}\mathrm{roots}\mathrm{of}{\psi }_m\left(t\right).\end{array}$

By evaluating Equations (47)–(48c) at collocation points $(x_i,t_j)$ defined in Equations (49a) and (49b), a system of $(n+1)\times (m+1)$ nonlinear algebraic equations can be formulated as

(50)$\left\{\begin{array}{c}\Phi (x_i,t_j)=0,i=1,\dots ,n-1,j=1,\dots ,m,\hfill \\ {\Lambda }_l\left(t_j\right)=0,l=1,2,j=1,\dots ,m,\hfill \\ \Pi \left(x_i\right)=0,i=0,\dots ,n.\hfill \end{array}\right.$

Hence, the system (50), comprising $(n+1)\times (m+1)$ nonlinear algebraic equations, facilitates the determination of the unknown coefficients ${\delta }_{i,j},i=0,1,\dots ,n,j=0,1,\dots ,m$. Solving this system enables the approximation of the solution $u_{n,m}(x,t)$ for the problem (30)–(31b) in the form (45).

When addressing highly nonlinear problems, the Newton–Raphson algorithm can demonstrate dependency on the selection of the initial conditions. In practical applications, we employ solutions derived from linearized versions of the original problem to provide initial approximations. Damping strategies are implemented as needed to guarantee convergence. Additionally, it should be noted that the Gauss–Legendre quadrature operates under the assumption of smooth weight functions $\vartheta \left(\alpha \right)$. When dealing with discontinuous or piecewise functions, precision may degrade in regions near discontinuities. Under these circumstances, the composite Gauss–Legendre quadrature incorporating subintervals around discontinuity locations or adaptive quadrature approaches to optimize node distribution [63] can be utilized. These adaptations ensure accurate approximation of the integral in Equation (19) while preserving the collocation structure. It is important to recognize that solving the nonlinear system (50) through the Newton–Raphson approach requires resolving a linear system of dimension $(n+1)(m+1)$ at each iteration. Considering the computation of fractional derivatives and quadrature terms at collocation nodes, the asymptotic computational complexity for Jacobian evaluation reaches $\mathcal{O}\left((n+1)(m+1)\overline{M}\right)$. Linear system resolution presents a worst-case complexity of $\mathcal{O}\left({(n+1)}^3{(m+1)}^3\right)$; however, the Jacobian’s structural characteristics, featuring sparse components from the HCPs basis, reduce this to approximately $\mathcal{O}\left({(n+1)}^2{(m+1)}^2\right)$ when utilizing sparse solvers in MATLAB (2019a). To reduce computational expenses for large values of n and m, a preconditioned conjugate gradient approach may be implemented to decrease iteration requirements. Enhancing $\overline{M}$, representing the quantity of Gauss–Legendre quadrature nodes used for approximating the distributed-order integral in Equation (19), enhances the precision of integral approximation, since the remainder term $r_{LG}$ in Equation (29) diminishes with increasing $\overline{M}$. For smooth $\vartheta \left(\alpha \right)$ functions (such as $\Gamma (3-\alpha )$ in Example 1), numerical experiments demonstrate that $\overline{M}=6$ produces errors below ${10}^{-8}$ when $n=m=10$. For non-smooth $\vartheta \left(\alpha \right)$ functions, including step functions, $\overline{M}\ge 8$ per subinterval within composite quadrature becomes necessary to preserve low error levels. Beyond $\overline{M}=10$, computational expenses increase substantially while error reduction becomes negligible due to the predominance of spectral approximation errors. Therefore, $\overline{M}$ = 7–8 typically represents the optimal choice for smooth cases, whereas composite quadrature with $\overline{M}$ = 8–10 per subinterval is advised for non-smooth scenarios.

4. Convergence Analysis

This section delves into investigating the convergence behavior exhibited by the numerical solution (45) towards the exact solution of the primary problem (30).

Theorem 4

([18]). Consider a function $u(x,t)\in L_w^2\left([0,x_{max}]\times [0,t_{max}]\right)$ with a weight function $w(x,t)=\varpi \left(x\right)\varpi \left(t\right)$, admitting the expansion $u(x,t)={\displaystyle \sum _{i=0}^{\infty }}{\displaystyle \sum _{j=0}^{\infty }}{\mathsf{u}}_{i,j}{\psi }_i\left(x\right){\psi }_j\left(t\right)$. Assume that ${∥{\displaystyle \frac{{\partial }^{6}f(x,t)}{\partial x^3\partial t^3}}∥}_2\le \overline{c}$ for a positive constant $\overline{c}$. Then, for expansion coefficients ${\mathsf{u}}_{i,j}$ with $i,j>3$, the inequality $|{\mathsf{u}}_{i,j}|<{\displaystyle \frac{\overline{c}}{i^3{j}^3}}$ holds. Furthermore, if

(51)$u_{n,m}(x,t)=\sum _{i=0}^n\sum _{j=0}^m{\mathsf{u}}_{i,j}{\psi }_i\left(x\right){\psi }_j\left(t\right),$

is an approximation of $u(x,t)$, then the following estimates hold:

(52a)$\begin{array}{ccc}\hfill |u(x,t)-f_{n,m}(x,t)|& <& {\displaystyle \frac{\overline{c}}{2^{n+m}}},\hfill \end{array}$

(52b)$\begin{array}{ccc}\hfill |u_x(x,t)-{\displaystyle \frac{\partial f_{n,m}}{\partial x}}(x,t)|& <& \kappa {\displaystyle \frac{n}{2^{n+m-2}}},\hfill \end{array}$

where κ is a positive constant.

Theorem 5.

Let $u(x,t)$ denote the exact solution of the differential Equation (30), and $u_{n,m}(x,t)$ represent its numerical approximation. Let $r_{n,m}(x,t)$ denote the residual error associated with this approximation. Then, it holds that:

(53) $|r_{n,m}(x,t)|<{\displaystyle \frac{\overline{\sigma }m+\mu }{2^{n+m-2}}}+r_{LG},$

where $\overline{\sigma }$ and μ are positive constants.

Proof. 

Based on the assumption, $u_{n,m}(x,t)$ satisfies the following equation:

(54)$\begin{array}{c}\hfill {\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}{u}_{n,m}(x,t)={\Theta }_l(x,t){}_{0 }{}^{RL}\mathcal{D}_x^{\beta }{u}_{n,m}(x,t)\\ \hfill +{\Theta }_r(x,t){}_{x }{}^{RL}\mathcal{D}_{x_{\max }}^{\beta }{u}_{n,m}(x,t)+u(x,t,u_{n,m}(x,t))+r_{n,m}(x,t).\end{array}$

The term $e_{n,m}(x,t):=u(x,t)-u_{n,m}(x,t)$ is considered as the error function, and hence combining Equations (30) and (54) results in

(55)$\begin{array}{c}\hfill |r_{n,m}(x,t)|\le |{\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}\left(e_{n,m}(x,t)\right)|+\left|{\Theta }_l(x,t)\right||{}_{0 }{}^{RL}\mathcal{D}_x^{\beta }\left(e_{n,m}(x,t)\right)|\\ \hfill +|{\Theta }_r(x,t)\left|\right|{}_{x }{}^{RL}\mathcal{D}_{x_{\max }}^{\beta }\left(e_{n,m}(x,t)\right)|+|u(x,t,u(x,t))-u(x,t,u_{n,m}(x,t))|.\end{array}$

Employing Definition 3 and Lemma 4 yields

(56)$\begin{array}{c}\hfill {\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}\left(e_{n,m}(x,t)\right)={\int }_0^1\vartheta \left(\alpha \right){}_0{}^C\mathcal{D}_t^{\alpha }\left(e_{n,m}(x,t)\right)\mathrm{d}\alpha \\ \hfill =\sum _{i=0}^{\overline{M}}{\varpi }_i^{\overline{M}}\vartheta \left({\varsigma }_i\right){}_0{}^C\mathcal{D}_t^{\varsigma i}\left(e_{n,m}(x,t)\right)+r_{LG}.\end{array}$

Additionally, using Definition 2 and Theorem 4 gives

(57)$\begin{array}{ccc}\hfill \left|{}_0{}^C\mathcal{D}_t^{\varsigma i}\left(e_{n,m}(x,t)\right)\right|& \le & {\int }_0^t{\displaystyle \frac{{|(t-s)}^{-{\varsigma }_i}|}{\Gamma (1-{\varsigma }_i)}}\left|{\displaystyle \frac{\partial }{\partial s}}{e}_{n,m}(x,s)\right|\mathrm{d}s\hfill \\ & <& {\displaystyle \frac{\kappa m}{\Gamma (1-{\varsigma }_i)2^{n+m-2}}}{\int }_0^t\left|{(t-s)}^{-{\varsigma }_i}\right|\mathrm{d}s,\hfill \end{array}$

where $\kappa$ is a positive constant. Given that $0<s<t\le t_{max}$, it follows that

(58)$\begin{array}{c}\hfill \left|{}_0{}^C\mathcal{D}_t^{\varsigma i}\left(e_{n,m}(x,t)\right)\right|<{\displaystyle \frac{\kappa t_{max}^{1-{\varsigma }_i}m}{\Gamma (1-{\varsigma }_i)2^{n+m-2}}},\end{array}$

Thus, combining Equations (56) and (58) results in

(59)$\begin{array}{c}\hfill \left|{\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}\left(e_{n,m}(x,t)\right)\right|<\overline{\sigma }{\displaystyle \frac{m}{2^{n+m-2}}}+r_{LG},\end{array}$

where

(60)$\overline{\sigma }=\sum _{i=0}^{\overline{M}}{\varpi }_i^{\overline{M}}\vartheta \left({\varsigma }_i\right){\displaystyle \frac{\kappa t_{max}^{1-{\varsigma }_i}}{\Gamma (1-{\varsigma }_i)}}.$

Utilizing Definition 1 and Theorem 4 obtains

(61)$\begin{array}{ccc}\hfill \left|{}_{0 }{}^{RL}\mathcal{D}_x^{\beta }\left(e_{n,m}(x,t)\right)\right|& \le & {\displaystyle \frac{{\partial }^2}{\partial x^2}}{\int }_0^x{\displaystyle \frac{{(x-s)}^{1-\beta }}{\Gamma (2-\beta )}}\left|e_{n,m}(s,t)\right|\mathrm{d}s\hfill \\ \hfill & <& {\displaystyle \frac{\overline{c}}{\Gamma (2-\beta )2^{n+m}}}{\displaystyle \frac{{\partial }^2}{\partial x^2}}{\int }_0^x{(x-s)}^{1-\beta }\mathrm{d}s\hfill \\ \hfill & =& {\displaystyle \frac{\overline{c}}{\Gamma (2-\beta )2^{m+n}}}(1-\beta )x^{-\beta }\hfill \\ & \le & {\displaystyle \frac{\overline{c}{\widehat{\theta }}_1}{\Gamma (1-\beta )2^{m+n}}},\hfill \end{array}$

where $\overline{c}$ is a positive constant and ${\widehat{\theta }}_1=\underset{x\in (0,x_{max}]}{max}\left\{x^{-\beta }\right\}$. Moreover, employing Definition 1 and Theorem 4 yields

(62)$\begin{array}{ccc}\hfill \left|{}_{x }{}^{RL}\mathcal{D}_{x_{\max }}^{\beta }\left(e_{n,m}(x,t)\right)\right|& \le & {\displaystyle \frac{{\partial }^2}{\partial x^2}}{\int }_x^{x_{max}}{\displaystyle \frac{{(s-x)}^{1-\beta }}{\Gamma (2-\beta )}}\left|e_{n,m}(s,t)\right|\mathrm{d}s\hfill \\ \hfill & <& {\displaystyle \frac{\overline{c}}{\Gamma (2-\beta )2^{n+m}}}{\displaystyle \frac{{\partial }^2}{\partial x^2}}{\int }_x^{x_{max}}{(s-x)}^{1-\beta }\mathrm{d}s\hfill \\ \hfill & =& {\displaystyle \frac{\overline{c}}{\Gamma (2-\beta )2^{n+m}}}(1-\beta ){(x_{max}-x)}^{-\beta }\hfill \\ \hfill & \le & {\displaystyle \frac{\overline{c}{\widehat{\theta }}_2}{\Gamma (1-\beta )2^{n+m}}},\hfill \end{array}$

where ${\widehat{\theta }}_2=\underset{x\in (0,x_{max}]}{max}\left\{{(x_{max}-x)}^{-\beta }\right\}$. The function F adheres to the Lipschitz condition concerning u, implying

(63)$|u(x,t,u(x,t))-u(x,t,u_{n,m}(x,t))|\le {\eta }_F\left|e_{n,m}(x,t)\right|,$

where ${\eta }_F$ is a positive real constant. Thus, using Theorem 4 gives

(64)$\begin{array}{c}\hfill |u(x,t,u(x,t))-u(x,t,u_{n,m}(x,t))|<{\displaystyle \frac{\overline{c}{\eta }_F}{2^{n+m}}}.\end{array}$

Let

(65)${\zeta }_l=\underset{(x,t)\in \Omega \times \mathcal{I}}{max}\left\{|{\Theta }_l(x,t)|\right\},{\zeta }_r=\underset{(x,t)\in \Omega \times \mathcal{I}}{max}\left\{|{\Theta }_r(x,t)|\right\},$

hence, considering the relations (55) and (59)–(64), it follows that

(66)$\begin{array}{ccc}\hfill |r_{n,m}(x,t)|& <& \overline{\sigma }{\displaystyle \frac{m}{2^{n+m-2}}}+\mu {\displaystyle \frac{1}{2^{n+m}}}+r_{LG}\hfill \\ & <& {\displaystyle \frac{\overline{\sigma }m+\mu }{2^{n+m-2}}}+r_{LG},\hfill \end{array}$

where

(67)$\mu ={\zeta }_l{\displaystyle \frac{\overline{c}{\widehat{\theta }}_1}{\Gamma (1-\beta )}}+{\zeta }_r{\displaystyle \frac{\overline{c}{\widehat{\theta }}_2}{\Gamma (1-\beta )}}+\overline{c}{\eta }_F.$

□

The differentiation matrices ${\mathsf{D}}_{\mathrm{L}}^{RL}\left(x\right)$ and ${\mathsf{D}}_{\mathrm{R}}^{RL}$ (Theorems 2 and 3) corresponding to the left- and right-sided Riemann–Liouville fractional derivatives are formulated using the explicit representations of HCPs (Equations (34) and (36)). The condition number of these matrices, expressed as

(68)$\kappa \left(\mathsf{D}\right)={\parallel \mathsf{D}\parallel }_2{\parallel {\mathsf{D}}^{-1}\parallel }_2,$

depends on the parameters n and $\beta$. For values of $n\le 15$, numerical experiments demonstrate that $\kappa \left({\mathsf{D}}_{\mathrm{L}}^{RL}\right)\approx {10}^3-{10}^4$ and $\kappa \left({\mathsf{D}}_{\mathrm{R}}^{RL}\right)\approx {10}^4-{10}^5$ when $\beta \in [1.3,1.8]$, suggesting reasonable conditioning attributed to the stable recurrence relationships inherent in HCPs [60].

5. Numerical Examples

This section applies the developed approach to solve the nonlinear two-sided spatial-fractional diffusion equation incorporating time-distributed order. To evaluate the method’s precision, the maximum error criterion is utilized, defined as

(69)$e_n=\underset{1\le i,j\le n}{max}|u(x_i,t_j)-u_n(x_i,t_j)|,$

where $u_n(x,t)$ denotes the approximate solution derived from Equation (45). The convergence rate (${\mathtt{CO}}_n$) is subsequently calculated as

(70)${\mathtt{CO}}_n={log}_{{\displaystyle \frac{n_1}{n_2}}}{\displaystyle \frac{e_{n_1}}{e_{n_2}}}.$

For nonlinear source terms, we compute the Jacobian matrix symbolically within the Newton–Raphson algorithm. For extended systems, symbolic sparsity structures are employed to enhance convergence speed.

All computational procedures are executed using MATLAB on an Intel (R) Core (TM) i5-4210U CPU @ 1.70 GHz processor. This configuration guarantees efficient and precise calculation of the outcomes.

The selection of parameters n, m, and $\overline{M}$ in Examples 1 and 2 is made to optimize the trade-off between precision and computational expense. These selections are determined through preliminary convergence evaluations to deliver accurate outcomes with manageable computational requirements. Uniform spatial and temporal discretization is adopted unless otherwise stated, meaning $n=m$. This approach ensures equivalent resolution in both spatial and temporal dimensions, which is suitable for the isotropic diffusion scenario. These parameter values are refined through error analysis across a parameter grid of $n,\overline{M}$ values, choosing the minimal parameters that satisfy the desired accuracy threshold.

Example 1.

Consider the following nonlinear two-sided spatial-fractional diffusion equation with time-distributed order [64]:

(71) $\begin{array}{ccc}\hfill {\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}u(x,t)& =& x^{0.6}(1+t){}_0{}^{RL}\mathcal{D}_x^{\beta }u(x,t)\hfill \\ \hfill & & +{(1-x)}^{0.6}(1+t){}_x{}^{RL}\mathcal{D}_1^{\beta }u(x,t)+\overline{\eta }{u}^2(x,t)+u(x,t),\hfill \end{array}$

over the domain $(x,t)\in [0,1]\times [0,1]$, where $\overline{\eta }$ represents a real constant, and $\vartheta \left(\alpha \right)=\Gamma (3-\alpha )$. The function $u(x,t)$ is defined as

(72) $\begin{array}{ccc}\hfill u(x,t)& =& {\displaystyle \frac{-2(t^2-t)}{ln\left(t\right)}}{x}^2(1-x^2)+\hfill \\ \hfill & & {\displaystyle \frac{2}{\Gamma (7-\beta )}}((t-1){(1+t)}^2{x}^{2.6-\beta }({\beta }^2-11\beta +3)\times \hfill \\ \hfill & & (12+12x^2+6x(\beta -4)-7\beta +{\beta }^2))\hfill \\ \hfill & & +{\displaystyle \frac{2(t-1)}{\Gamma (5-\beta )}}({(1+t)}^2{(1-x)}^{2.6-\beta }\hfill \\ & & \times (12x^2-6x\beta +(-1+\beta )\beta )-\overline{\eta }{x}^4{(1-x)}^4{(1-t^2)}^2.\hfill \end{array}$

Notably, the exact solution of this example is given by:

(73) $u(x,t)=x^2{(1-x)}^2(1-t^2).$

Figure 1 depicts the comparison between the exact and numerical solutions across various time levels, considering parameters $\beta =1.6$, $\overline{\eta }=2$, $t_{max}=1.5$, $n=8$, and $\overline{M}=7$.

Figure 2 illustrates the surface plot of absolute error (top) and the corresponding contour plot (bottom) for two scenarios with $n=10$ and $n=15$. These plots are generated under the conditions of $\beta =1.5$, $\overline{\eta }=1$, $t_{max}=1$, and $\overline{M}=7$. These visualizations offer insights into the accuracy and behavior of the numerical solution under varying conditions.

A comparative study is also performed against a fast implicit difference (FID) method outlined in [64], utilizing parameters such as $h={\displaystyle \frac{x_{max}}{M}}$, $\tau ={\displaystyle \frac{t_{max}}{N}}$, ${\mathtt{CO}}_{\tau }={log}_2{\displaystyle \frac{e(h,\tau ,\sigma )}{e(h,\tau /2,\sigma )}}$, $\sigma =0.2$, and $M=2000$. Table 1 provides a comprehensive analysis, presenting the maximum error, CPU time (in seconds), and convergence order for different values of n when $\beta =1.3$ and $1.8$. These evaluations are conducted under conditions where $\overline{\eta }=0$, $t_{max}=1.5$, and $\overline{M}=8$.

When evaluating computational efficiency by normalizing accuracy against processing time, the proposed approach demonstrates significant advantages over the FID method referenced in the HYJian citation. Using $\beta =1.3$, the new method ($n=12$) produces an error of approximately $1.3\times {10}^{-12}$ within 13.012 s, whereas the FID approach ($N=20$, $M=2000$) generates an error of roughly $8.7\times {10}^{-5}$ in an estimated 15–20 s based on complexity analysis from the HYJian study. This represents a remarkable improvement of about seven orders of magnitude in accuracy for equivalent computational investment. Similar performance patterns emerge with $\beta =1.8$, where the proposed method achieves $3.3\times {10}^{-12}$ error ($n=12$, 12.842 s) compared to the FID method’s $4.9\times {10}^{-5}$ error ($N=20$, approximately 15–20 s). These results demonstrate the proposed method’s exceptional precision-to-cost ratio, attributed to its exponential convergence characteristics.

Higher-order implementations (such as $n=20$) present conditioning challenges, with the condition number $\kappa$ from Equation (68) escalating to ${10}^5$–${10}^6$ ranges due to polynomial basis ill-conditioning at elevated degrees. Additionally, as $\beta$ approaches 2, mathematical expressions including $\Gamma (r+1-\beta )$ terms in Equation 34 encounter near-singular behavior, amplifying $\kappa$ by factors of 2–5 when comparing $\beta =1.95$ to $\beta =1.5$. To address these stability issues, diagonal preconditioning techniques are implemented, achieving approximately 30% reduction in effective condition numbers for cases with $n=15$ and $\beta =1.8$, thereby maintaining numerical solution reliability as documented in Table 1.

Example 2.

Consider the following nonlinear two-sided spatial-fractional diffusion equation with time-distributed order:

(74) ${\mathcal{D}}_{0,t}^{\vartheta \left(\alpha \right)}u(x,t)=x^2{e}^{-t}{}_0{}^{RL}\mathcal{D}_x^{\beta }u(x,t)+x^2{e}^{\top }{}_x{}^{RL}\mathcal{D}_1^{\beta }u(x,t)+sin\left(u(x,t)\right)+u(x,t),$

over the domain $(x,t)\in [0,1]\times [0,1]$, where $\vartheta \left(\alpha \right)=\Gamma (4-\alpha )sinh\left(\alpha \right)$ and

(75) $\begin{array}{c}\hfill u(x,t)={\displaystyle \frac{1}{2\Gamma (2-\beta )}}{e}^{t-x-1}t(t-1)x^2{(-{(x-1)}^2)}^{-\beta }\left\{e{(-{(x-1)}^2)}^{\beta }\right((-1+{(-1)}^{\beta }{e}^{2x})\\ \hfill \times \Gamma (2-\beta )+\Gamma (2-\beta ,1-x))-e^x({(-1)}^{\beta }{e}^2{(1-x)}^{\beta }(\beta -x)-(x-2+\beta )\\ \hfill \times {(x-1)}^{\beta }+{(-1)}^{\beta }{e}^{x+1}{(-{(x-1)}^2)}^{\beta }\Gamma (2-\beta ,x-1)\left)\right\}+{\displaystyle \frac{1}{\Gamma (8-\beta )}}\{e^{-t}t(t-1)\\ \hfill \times x^{3-\beta }({\mathcal{F}}_2^1(\left[1\right];[2-{\displaystyle \frac{\beta }{2}},{\displaystyle \frac{5}{2}}-{\displaystyle \frac{\beta }{2}}];{\displaystyle \frac{x^2}{4}})(5040-8028\beta +5104{\beta }^2-1665{\beta }^3+295{\beta }^4\\ \hfill -27{\beta }^5+{\beta }^6)+2{\mathcal{F}}_2^1(\left[2\right];[3-{\displaystyle \frac{\beta }{2}},{\displaystyle \frac{7}{2}}-{\displaystyle \frac{\beta }{2}}];{\displaystyle \frac{x^2}{4}})x^2(294-175\beta +33{\beta }^2-2{\beta }^3)\\ \hfill +8x^4{\mathcal{F}}_2^1(\left[3\right];[4-{\displaystyle \frac{\beta }{2}},{\displaystyle \frac{9}{2}}-{\displaystyle \frac{\beta }{2}}];{\displaystyle \frac{x^2}{4}})\left)\right\}-sin(t(1-t)sinh\left(x\right))+{\displaystyle \frac{1}{{(-1+ln{\left(t\right)}^2)}^3}}\\ \hfill \times sinh\left(x\right)\{-6t^2-4cosh\left(1\right)(e^{-1}+2t(5-2t)+e(4t-7))ln\left(t\right)+2(6t^2\\ \hfill {-cosh\left(1\right)-t(3+4cosh\left(1\right)))ln\left(t\right)}^2-3sinh\left(1\right)+2sinh\left(1\right)(2t-1)ln{\left(t\right)}^5+2t\\ \hfill \times (4+2cosh\left(1\right)+sinh\left(1\right))+\left((6+{\displaystyle \frac{3}{e}}+e)t-6t^2-2cosh\left(1\right)+3sinh\left(1\right)\right)ln{\left(t\right)}^4\\ \hfill +2ln{\left(t\right)}^3(2t^2+3cosh\left(1\right)+sinh\left(1\right))-t(5+2cosh\left(1\right)+4sinh\left(1\right))\}sinh\left(x\right),\end{array}$

where ${\mathcal{F}}_q^p([a_1,\dots ,a_p];[b_1,\dots ,b_q];t)$ denotes the generalized hypergeometric function. It is important to observe that the exact solution for this example is given by

(76) $u(x,t)=t(1-t)sinh\left(x\right).$

Figure 3 exhibits the exact and numerical solutions of $u(x,t)$ for the scenario where $beta=1.85$, $n=10$, and $barM=6$. It visually compares the computed and exact solutions, offering a comprehensive understanding of the numerical accuracy achieved.

Table 2 provides a deeper analysis of detailed information regarding the maximum error, convergence order, and CPU-time (s) for different values of n when $\beta =1.25$ and $\beta =1.75$. This table serves as a valuable resource for evaluating the performance of the numerical method under various parameter settings, shedding light on its efficiency and accuracy across different scenarios.

Furthermore, Figure 4 presents a comprehensive analysis of the surface of absolute error and the contour plot for specific cases with $n=8$ and $n=12$. By visualizing the error distribution and contour patterns, this figure provides additional insights into the behavior of the numerical solution, particularly in scenarios where $\beta =1.5$ and $\overline{M}=8$. This extended analysis enhances understanding of the method’s performance and sensitivity to different parameter combinations.

6. Conclusions

This paper presented a collocation method based on Hexic-kind Chebyshev polynomials (HCPs) for solving nonlinear diffusion equations involving time-distributed-order and two-sided space-fractional terms. The inherent complexity of such equations, which arise from the nonlinearity and fractional nature of the derivatives, makes traditional analytical methods impractical. A numerical approach was developed to address this, transforming the problem into a set of nonlinear algebraic equations that can be efficiently solved using the Newton–Raphson method. HCPs in the collocation method proved highly effective due to their strong recurrence relations and differentiation properties, allowing for accurate approximations with minimal basis functions. A comprehensive convergence analysis was conducted, and the derived error bounds provided theoretical guarantees for the method’s robustness and accuracy. The approach was further validated through two test examples, where the numerical results closely aligned with existing solutions, demonstrating the method’s accuracy and reliability. For future research, extending the method to handle more complex equations with varying coefficients or mixed boundary conditions, exploring its application in multi-dimensional settings, and investigating alternative numerical schemes for comparison could offer valuable insights for further advancements in this field.

Author Contributions

A.B. led the investigation, providing expertise in synthesizing and analyzing the literature; S.B. contributed to software development and methodology, creating tools for data analysis and establishing the systematic review protocol; B.P.M. contributed to conceptualization, offering valuable insights into framing the review’s scope and objectives; A.D. was responsible for checking for technical errors, for software development, for methodology, for completing a comprehensive review, and for editing to verify the manuscript’s clarity and cohesion. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study. Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 The comparison between the exact (solid line) and numerical (dashed line) solutions of Example 1 across various time levels.

Figure 2 Surface (a,c) and contour (b,d) plots of absolute error for Example 1 with $\beta =1.5$, $\overline{\eta }=1$, $t_{max}=1$, $\overline{M}=7$. (a,b) $n=10$, error range: ${10}^{-7}$ to ${10}^{-10}$; (c,d) $n=15$, error range: ${10}^{-9}$ to ${10}^{-12}$.

Figure 3 Comparison of exact (solid line) and numerical (dashed line) solutions at different time levels for Example 2.

Figure 4 Surface (a,c) and contour (b,d) plots of absolute error for Example 2 with $\beta =1.5$, $\overline{\eta }=1$, $t_{max}=1$, $\overline{M}=8$. (a,b) $n=8$, error range: ${10}^{-6}$ to ${10}^{-9}$; (c,d) $n=12$, error range: ${10}^{-8}$ to ${10}^{-11}$.

Table 1

Comparison of maximum errors ($e(h,\tau ,\sigma )$ for FID [64], $e_n$ for proposed method), convergence orders (CO_τ, CO_n), and CPU times (seconds, proposed method only) for Example 1. Parameters: $\beta =1.3,1.8$, $\overline{\eta }=0$, $t_{\max }=1.5$, $\overline{M}=8$. FID uses $h=x_{\max }/M$, $\tau =t_{\max }/N$, $\sigma =0.2$, $M=2000$. FID CPU times estimated from [64].

	FID Method [64] with $\beta =1.3$		Proposed method with $\beta =1.3$
N	$e(h,\tau ,\sigma )$	${\mathtt{CO}}_{\tau }$	n	$e_n$	${\mathtt{CO}}_n$	CPU-time
5	$6.897311\times {10}^{-4}$	−	6	$2.3274\times {10}^{-6}$	−	$5.272$
10	$2.640520\times {10}^{-4}$	$1.3852$	9	$4.1124\times {10}^{-9}$	$15.6327$	$9.241$
20	$8.702462\times {10}^{-5}$	$1.6013$	12	$1.2569\times {10}^{-12}$	$28.1319$	$13.012$
40	$4.570569\times {10}^{-5}$	$0.9290$	15	$3.0292\times {10}^{-16}$	$37.3336$	$17.411$
	FID Method [64] with $\beta =1.8$		Proposed method with $\beta =1.8$
N	$e(h,\tau ,\sigma )$	${\mathtt{CO}}_{\tau }$	n	$e_n$	${\mathtt{CO}}_n$	CPU-time
5	$2.559967\times {10}^{-4}$	−	6	$4.1421\times {10}^{-6}$	−	$6.062$
10	$1.099115\times {10}^{-4}$	$1.2198$	9	$6.7791\times {10}^{-9}$	$15.8216$	$9.421$
20	$4.935518\times {10}^{-5}$	$1.1551$	12	$3.3011\times {10}^{-12}$	$26.5130$	$12.842$
40	$2.399388\times {10}^{-5}$	$1.0405$	15	$2.4761\times {10}^{-15}$	$32.2453$	$18.332$

Table 2

Impact of varying $\beta$ and n on maximum error, convergence order, and CPU time (in seconds) in Example 2.

	$\beta =1.25$			$\beta =1.75$
n	$e_n$	${\mathtt{CO}}_n$	CPU-time	$e_n$	${\mathtt{CO}}_n$	CPU-time
3	$1.3411\times {10}^{-2}$	−	$4.137$	$2.1578\times {10}^{-4}$	−	$4.056$
6	$6.2078\times {10}^{-7}$	$14.399$	$7.556$	$1.1431\times {10}^{-7}$	$10.8823$	$7.379$
9	$8.9377\times {10}^{-11}$	$21.8166$	$11.561$	$7.1387\times {10}^{-11}$	$18.1979$	$11.467$
12	$1.0362\times {10}^{-13}$	$23.4976$	$24.006$	$1.2120\times {10}^{-13}$	$22.1717$	$23.973$

References

1. Podlubny, I. Fractional Differential Equations; Academic Press: San Diego, CA, USA, 1999.

2. Metzler, R.; Krafter, J. The random walk’s guide to anomalous diffusion: A fractional dynamic approach. Phys. Rep.; 2000; 339, pp. 1-72. [DOI: https://dx.doi.org/10.1016/S0370-1573(00)00070-3]

3. Mathieu, B.; Melchior, P.; Oustaloup, A.; Ceyral, C. Fractional differentiation for edge detection. Signal Process.; 2003; 83, pp. 2421-2432. [DOI: https://dx.doi.org/10.1016/S0165-1684(03)00194-4]

4. Baleanu, D.; Diethelm, K.; Scalas, E.; Trujillo, J.J. Fractional Calculus Models and Numerical Methods; World Scientific: Singapore, 2017.

5. Ionescu, C.; Lopes, A. The role of fractional calculus in modeling biological phenomena: A review. Commun. Nonlinear Sci. Numer. Simul.; 2017; 51, pp. 141-159. [DOI: https://dx.doi.org/10.1016/j.cnsns.2017.04.001]

6. Amin, A.; Lopes, A.M. Fractional-order shifted Legendre collocation method for solving non-linear variable-order fractional Fredholm integro-differential equations. Comput. Appl. Math.; 2022; 41, 2.

7. Owolabi, K.; Atangana, A. Numerical approximation of nonlinear fractional parabolic differential equations with Caputo–Fabrizio derivative in Riemann-Liouville sense. Chaos Soliton. Fract.; 2017; 99, pp. 171-179. [DOI: https://dx.doi.org/10.1016/j.chaos.2017.04.008]

8. Babaei, A.; Banihashemi, S. A stable numerical approach to solve a time-fractional inverse heat conduction problem. Iran. J. Sci. Technol. Trans. A.; 2017; 42, pp. 2225-2236. [DOI: https://dx.doi.org/10.1007/s40995-017-0360-4]

9. Babaei, A.; Banihashemi, S. Reconstructing unknown nonlinear boundary conditions in a time-fractional inverse reaction-diffusion-convection problem. Numer. Methods Part. Differ. Equ.; 2018; 35, pp. 976-992. [DOI: https://dx.doi.org/10.1002/num.22334]

10. Babaei, A.; Jafari, H.; Banihashemi, S. A collocation approach for solving time-fractional stochastic heat equation driven by an additive noise. Symmetry; 2020; 12, 904. [DOI: https://dx.doi.org/10.3390/sym12060904]

11. El-Ajou, A.; Moa’ath, N.; Al-Zhour, Z.; Momani, S. A class of linear non-homogenous higher order matrix fractional differential equations: Analytical solutions and new technique. Fract. Calc. Appl. Anal.; 2020; 23, pp. 356-377. [DOI: https://dx.doi.org/10.1515/fca-2020-0017]

12. Bologna, M.; Tsallis, C.; Grigolini, P. Anomalous diffusion associated with nonlinear fractional derivative Fokker-Planck-like equation: Exact time-dependent solutions. Phys. Rev. E; 2000; 62, 2213. [DOI: https://dx.doi.org/10.1103/PhysRevE.62.2213]

13. Pramukkul, P.; Svenkeson, A.; Grigolini, P.; Bologna, M.; West, B. Complexity and the fractional calculus. Adv. Math. Phys.; 2013; 2013, 498789. [DOI: https://dx.doi.org/10.1155/2013/498789]

14. Moghaddam, B.P.; Zaky, M.A.; Lopes, A.M.; Galhano, A. A Fractional Time-Space Stochastic Advection-Diffusion Equation for Modeling Atmospheric Moisture Transport at Ocean-Atmosphere Interfaces. Fractal Fract.; 2025; 9, 211. [DOI: https://dx.doi.org/10.3390/fractalfract9040211]

15. Lorenzo, C.F.; Hartley, T.T. Variable order and distributed order fractional operators. Nonlinear Dyn.; 2002; 29, pp. 57-98. [DOI: https://dx.doi.org/10.1023/A:1016586905654]

16. Tavares, D.; Almeida, R.; Torres, D.F.M. Caputo derivatives of fractional variable order: Numerical approximations. Commun. Nonlinear Sci. Numer. Simul.; 2016; 35, pp. 69-87. [DOI: https://dx.doi.org/10.1016/j.cnsns.2015.10.027]

17. Sweilam, N.; AL-Mekhlafi, S.; Shatta, S.; Baleanu, D. Numerical Study for a Novel Variable-Order Multiple Time Delay Awareness Programs Mathematical Model. Appl. Numer. Math.; 2020; 158, pp. 212-235. [DOI: https://dx.doi.org/10.1016/j.apnum.2020.07.016]

18. Babaei, A.; Jafari, H.; Banihashemi, S. Numerical solution of variable order fractional nonlinear quadratic integro-differential equations based on the sixth-kind Chebyshev collocation method. J. Comput. Appl. Math.; 2020; 377, 112908. [DOI: https://dx.doi.org/10.1016/j.cam.2020.112908]

19. Moniri, Z.; Babaei, A.; Moghaddam, B.P. Robust numerical framework for simulating 2D fractional time-space stochastic diffusion equation driven by spatio-temporal noise: L1-FFT hybrid approach. Commun. Nonlinear Sci. Numer. Simul.; 2025; 146, 108791. [DOI: https://dx.doi.org/10.1016/j.cnsns.2025.108791]

20. Moghaddam, B.P.; Dabiri, A.; Machado, J.A.T. Application of variable-order fractional calculus in solid mechanics. Appl. Eng. Life Soc. Sci. Part A; 2019; 7, pp. 207-224.

21. Moghaddam, B.P.; Babaei, A.; Dabiri, A.; Galhano, A. Fractional Stochastic Partial Differential Equations: Numerical Advances and Practical Applications—A State of the Art Review. Symmetry; 2024; 16, 563. [DOI: https://dx.doi.org/10.3390/sym16050563]

22. Machado, J.T. Optimal variable-order fractional PID controllers for dynamical systems. J. Comput. Appl. Math.; 2018; 339, pp. 40-48.

23. Ding, W.; Patnaik, S.; Sidhardh, S.; Semperlotti, F. Applications of distributed-order fractional operators: A review. Entropy; 2021; 23, 110. [DOI: https://dx.doi.org/10.3390/e23010110] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33467618]

24. Yang, X.J.; Machado, J.A.T.; Srivastava, H.M. A new numerical technique for solving the local fractional diffusion equation: Two-dimensional extended differential transform approach. Appl. Math. Comput.; 2016; 274, pp. 143-151. [DOI: https://dx.doi.org/10.1016/j.amc.2015.10.072]

25. Sun, H.G.; Chen, W.; Chen, Y. Variable-order fractional differential operators in anomalous diffusion modeling. Phys. A Stat. Mech. Its Appl.; 2009; 388, pp. 4586-4592. [DOI: https://dx.doi.org/10.1016/j.physa.2009.07.024]

26. Obembe, A.D.; Hossain, M.E.; Abu-Khamsin, S.A. Variable-order derivative time fractional diffusion model for heterogeneous porous media. J. Pet. Sci. Eng.; 2017; 152, pp. 391-405. [DOI: https://dx.doi.org/10.1016/j.petrol.2017.03.015]

27. Samko, S.G.; Ross, B. Integration and differentiation to a variable fractional order. Integral Transform. Spec. Funct.; 1993; 1, pp. 277-300. [DOI: https://dx.doi.org/10.1080/10652469308819027]

28. Caputo, M. Distributed order differential equations modelling dielectric induction and diffusion. Fract. Calc. Appl. Anal.; 2001; 4, pp. 421-442.

29. Caputo, M. Diffusion with space memory modelled with distributed order space fractional differential equations. Ann. Geophys.; 2003; 46, [DOI: https://dx.doi.org/10.4401/ag-3395]

30. Chechkin, A.; Gorenflo, R.; Sokolov, I.M. Retarding subdiffusion and accelerating superdiffusion governed by distributed-order fractional diffusion equations. Phys. Rev. E; 2002; 66, 046129. [DOI: https://dx.doi.org/10.1103/PhysRevE.66.046129]

31. Duan, J.; Baleanu, D. Steady periodic response for a vibration system with distributed order derivatives to periodic excitation. J. Vib. Control.; 2018; 24, pp. 3124-3131. [DOI: https://dx.doi.org/10.1177/1077546317700989]

32. Srokowski, T. Levy flights in nonhomogeneous media: Distributed-order fractional equation approach. Phys. Rev. E; 2008; 78, 031135. [DOI: https://dx.doi.org/10.1103/PhysRevE.78.031135]

33. Atanackovic, T.M.; Pilipovic, S.; Zorica, D. Time distributed-order diffusion-wave equation. I. Volterra-type equation. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci.; 2009; 465, pp. 1869-1891. [DOI: https://dx.doi.org/10.1098/rspa.2008.0445]

34. Solís-Pérez, J.; Gómez-Aguilar, J.; Atangana, A. Novel numerical method for solving variable-order fractional differential equations with power, exponential and Mittag-Leffler laws. Chaos Solitons Fractals; 2018; 114, pp. 175-185. [DOI: https://dx.doi.org/10.1016/j.chaos.2018.06.032]

35. Diethelm, K.; Ford, N.J. Numerical solution methods for distributed order differential equations. Fract. Calc. Appl. Anal.; 2001; 4, pp. 531-542.

36. Mainardi, F.; Pagnini, G.; Gorenflo, R. Some aspects of fractional diffusion equations of single and distributed order. Appl. Math. Comput.; 2007; 187, pp. 295-305. [DOI: https://dx.doi.org/10.1016/j.amc.2006.08.126]

37. Diethelm, K.; Ford, N.J. Numerical analysis for distributed-order differential equations. J. Comput. Appl. Math.; 2009; 225, pp. 96-104. [DOI: https://dx.doi.org/10.1016/j.cam.2008.07.018]

38. Luchko, Y. Boundary value problems for the generalized time-fractional diffusion equation of distributed order. Fract. Calc. Appl. Anal.; 2009; 12, pp. 409-422.

39. Gorenflo, R.; Luchko, Y.; Stojanovic, M. Fundamental solution of a distributed order time-fractional diffusion-wave equation as probability density. Fract. Calc. Appl. Anal.; 2013; 16, pp. 297-316. [DOI: https://dx.doi.org/10.2478/s13540-013-0019-6]

40. Duong, P.L.T.; Kwok, E.; Lee, M. Deterministic analysis of distributed order systems using operational matrix. Appl. Math. Model.; 2016; 40, pp. 1929-1940. [DOI: https://dx.doi.org/10.1016/j.apm.2015.09.035]

41. Naber, M. Distributed order fractional sub-diffusion. Fractals; 2004; 12, pp. 23-32. [DOI: https://dx.doi.org/10.1142/S0218348X04002410]

42. Sabermahani, S.; Ordokhani, Y.; Rahimkhani, P. Application of generalized Lucas wavelet method for solving nonlinear fractal-fractional optimal control problems. Chaos Solitons Fractals; 2023; 170, 113348. [DOI: https://dx.doi.org/10.1016/j.chaos.2023.113348]

43. Kharazmi, E.; Zayernouri, M.; Karniadakis, G.E. Petrov–Galerkin and spectral collocation methods for distributed order differential equations. SIAM J. Sci. Comput.; 2017; 39, pp. 1003-1037. [DOI: https://dx.doi.org/10.1137/16M1073121]

44. Zaky, M.A. A Legendre collocation method for distributed-order fractional optimal control problems. Nonlinear Dyn.; 2018; 91, pp. 2667-2681. [DOI: https://dx.doi.org/10.1007/s11071-017-4038-4]

45. Li, X.; Rui, H. A block-centered finite difference method for the distributed-order time-fractional diffusion-wave equation. Appl. Numer. Math.; 2018; 131, pp. 123-139. [DOI: https://dx.doi.org/10.1016/j.apnum.2018.04.013]

46. Rahimkhani, P.; Ordokhani, Y.; Lima, P.M. An improved composite collocation method for distributed-order fractional differential equations based on fractional Chelyshkov wavelets. Appl. Numer. Math.; 2019; 145, pp. 1-27. [DOI: https://dx.doi.org/10.1016/j.apnum.2019.05.023]

47. Pourbabaee, M.; Saadatmandi, A. A novel Legendre operational matrix for distributed order fractional differential equations. Appl. Math. Comput.; 2019; 361, pp. 215-231. [DOI: https://dx.doi.org/10.1016/j.amc.2019.05.030]

48. Arianfar, M.; Moghaddam, B.P.; Babaei, A. Computational technique for a class of nonlinear distributed-order fractional boundary value problems with singular coefficients. Comput. Appl. Math.; 2021; 40, 190. [DOI: https://dx.doi.org/10.1007/s40314-021-01576-6]

49. Moghaddam, B.P.; Tenreiro Machado, J.A.; Morgado, M.L. Numerical approach for a class of distributed order time fractional partial differential equations. Appl. Numer. Math.; 2019; 136, pp. 152-162. [DOI: https://dx.doi.org/10.1016/j.apnum.2018.09.019]

50. Bu, W.; Ji, L.; Tang, Y.; Zhou, J. Space-time finite element method for the distributed-order time fractional reaction diffusion equations. Appl. Numer. Math.; 2020; 152, pp. 446-465. [DOI: https://dx.doi.org/10.1016/j.apnum.2019.11.010]

51. Abbaszadeh, M.; Dehghan, M. A POD-based reduced-order Crank-Nicolson/fourth-order alternating direction implicit (ADI) finite difference scheme for solving the two-dimensional distributed-order Riesz space-fractional diffusion equation. Appl. Numer. Math.; 2020; 158, pp. 271-291. [DOI: https://dx.doi.org/10.1016/j.apnum.2020.07.020]

52. Kazmi, K.; Khaliq, A.Q. An efficient split-step method for distributed-order space-fractional reaction-diffusion equations with time-dependent boundary conditions. Appl. Numer. Math.; 2020; 147, pp. 142-160. [DOI: https://dx.doi.org/10.1016/j.apnum.2019.08.019]

53. A Chebyshev spectral method based on operational matrix for initial and boundary value problems of fractional order. Comput. Math. Appl.; 2011; 62, pp. 2364-2373. [DOI: https://dx.doi.org/10.1016/j.camwa.2011.07.024]

54. Dabiri, A.; Butcher, E.A.; Nazari, M. Coefficient of restitution in fractional viscoelastic compliant impacts using fractional Chebyshev collocation. J. Sound Vib.; 2017; 388, pp. 230-244. [DOI: https://dx.doi.org/10.1016/j.jsv.2016.10.013]

55. Khader, M.; Saad, K. A numerical approach for solving the fractional Fisher equation using Chebyshev spectral collocation method. Chaos Solitons Fractals; 2018; 110, pp. 169-177. [DOI: https://dx.doi.org/10.1016/j.chaos.2018.03.018]

56. Srivastava, H.; Saad, K.M.; Khader, M. An efficient spectral collocation method for the dynamic simulation of the fractional epidemiological model of the Ebola virus. Chaos Solitons Fractals; 2020; 140, 110174. [DOI: https://dx.doi.org/10.1016/j.chaos.2020.110174]

57. Dabiri, A.; Butcher, E.A. Efficient modified Chebyshev differentiation matrices for fractional differential equations. Commun. Nonlinear Sci. Numer. Simul.; 2017; 50, pp. 284-310. [DOI: https://dx.doi.org/10.1016/j.cnsns.2017.02.009]

58. Dabiri, A.; Butcher, E.A. Stable fractional Chebyshev differentiation matrix for the numerical solution of multi-order fractional differential equations. Nonlinear Dyn.; 2017; 90, pp. 185-201. [DOI: https://dx.doi.org/10.1007/s11071-017-3654-3]

59. Ezz-Eldien, S. Efficient Chebyshev spectral methods for solving multi-term fractional orders differential equations. Appl. Math. Model.; 2011; 35, pp. 5662-5672.

60. Babaei, A.; Banihashemi, S.; Parsa Moghaddam, B.; Dabiri, A. An algorithmic exploration of variable order fractional partial integro-differential equations via hexic shifted Chebyshev polynomials. Comput. Appl. Math.; 2024; 43, pp. 1-18. [DOI: https://dx.doi.org/10.1007/s40314-024-02703-9]

61. Masjed-Jamei, M. A basic class of symmetric orthogonal polynomials using the extended Sturm–Liouville theorem for symmetric functions. J. Math. Anal. Appl.; 2007; 325, pp. 753-775. [DOI: https://dx.doi.org/10.1016/j.jmaa.2006.02.007]

62. Abd-Elhameed, W.M.; Youssri, Y.H. Sixth-kind Chebyshev spectral approach for solving fractional differential equations. J. Nonlinear Sci. Numer. Simul.; 2019; 20, pp. 191-203. [DOI: https://dx.doi.org/10.1515/ijnsns-2018-0118]

63. Kahaner, D.; Moler, C.; Nash, S. Numerical Methods and Software; Prentice-Hall Series in Computational Mathematics; Prentice-Hall: Englewood Cliffs, NJ, USA, 1989.

64. Jian, H.Y.; Huang, T.Z.; Zhao, X.L.; Zhao, Y.L. A fast implicit difference scheme for a new class of time distributed-order and space fractional diffusion equations with variable coefficients. Adv. Differ. Equ.; 2018; 2018, 205. [DOI: https://dx.doi.org/10.1186/s13662-018-1655-2]

65. Hu, X.; Liu, F.; Turner, I.; Anh, V. An implicit numerical method of a new time distributed-order and two-sided space-fractional advection-dispersion equation. Numer. Algorithms; 2016; 72, pp. 393-407. [DOI: https://dx.doi.org/10.1007/s11075-015-0051-1]