Finite Difference/Fractional Pertrov–Galerkin Spectral Method for Linear Time-Space Fractional Reaction–Diffusion Equation

Abstract

Achieving high-order accuracy in finite difference/spectral methods for space-time fractional differential equations often relies on very restrictive and usually unrealistic smoothness assumptions in the spatial and/or temporal domains. For spatial discretization, spectral methods using smooth basis functions are commonly employed. However, spatial–fractional derivatives pose challenges, as they often lack guaranteed spatial smoothness, requiring non-smooth basis functions. In the temporal domain, finite difference schemes on uniformly graded meshes are commonly employed; however, achieving accuracy remains challenging for non-smooth solutions. In this paper, an efficient algorithm is adopted to improve the accuracy of finite difference/Pertrov–Galerkin spectral schemes for a time-space fractional reaction–diffusion equation, with a hyper-singular integral fractional Laplacian and non-smooth solutions in both time and space domains. The Pertrov–Galerkin spectral method is adapted using non-smooth generalized basis functions to discretize the spatial variable, and the L1 scheme on a non-uniform graded mesh is used to approximate the Caputo fractional derivative. The unconditional stability and convergence are established. The rate of convergence is $O\left(N^{\mu -\gamma }+K^{-min\{\rho \beta ,2-\beta \}}\right),$ achieved without requiring additional regularity assumptions on the solution. Finally, numerical results are provided to validate our theoretical findings.

1. Introduction

Reaction–diffusion models have widespread applicability across numerous fields, including chemistry, biology, physics, and engineering [1]. In these models, the Laplace operator typically represents the spatial diffusion that characterizes transport mechanisms arising from Brownian motion. However, several experiments have shown that many complex systems and phenomena involve anomalous diffusion, where the underlying stochastic processes are non-Brownian and are typically characterized by Lévy motion [2,3]. This distinction highlights the limitations of classical integer-order models in capturing the dynamics of such systems. To address these limitations, fractional models have been proposed, where the diffusion process is governed by the fractional Laplacian [2]. These fractional models are particularly well-suited for representing long-range interactions and power law invasion profiles observed in numerous applications. The fractional Laplacian, from a probabilistic perspective, corresponds to the infinitesimal generator of a stable Lévy process [4].

The d-dimensional fractional Laplacian exhibits rotational invariance, a critical property used in accurately modeling isotropic anomalous diffusion in various scientific and engineering applications [5]. This invariance ensures that the diffusion process behaves uniformly in all directions, making the fractional Laplacian particularly suited for systems where isotropy is essential. This distinguishes the fractional Laplacian from other types of fractional derivatives, which often lack rotational invariance and may not be appropriate for capturing the behavior of isotropic processes [6]. The fractional diffusion equation with both time and space fractional derivatives has been widely applied to study the anomalous diffusive processes associated with sub-diffusion (fractional in time) and super-diffusion (fractional in space) across many fields. Despite their advantages, the nonlocal nature of fractional operators poses substantial challenges in both theoretical analysis and numerical computation [7].

In this paper, we consider the following time-fractional reaction–diffusion equation with a hyper-singular integral fractional Laplacian:

(1)$\begin{array}{cc}\hfill {}_0{}^{C}D_t^{\beta }\varphi (x,t)+{(-\Delta )}^{\mu /2}\varphi (x,t)+\lambda \varphi (x,t)=& f(x,t),x\in I,t\in (0,T],\hfill \\ \hfill \varphi (x,t)=& 0,x\in I^c,t\in (0,T],\hfill \\ \hfill \varphi (x,0)=& g_0\left(x\right),x\in I,\hfill \end{array}$

where $I=(-1,1)$, $I^c=\mathbb{R}\setminus I$, and ${}_0{}^{C}D_t^{\beta }$ is the Caputo derivative of order $\beta \in (0,1)$ defined as follows:

(2)${}_0{}^{C}D_t^{\beta }\varphi \left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }(1-\beta )}}{\int }_0^t{(t-s)}^{-\beta }{\displaystyle \frac{\partial }{\partial s}}\varphi \left(s\right)ds.$

The the fractional Laplacian ${(-\Delta )}^{\mu /2}$ is defined by the following:

$\begin{array}{c}\hfill {(-\Delta )}^{\mu /2}\varphi \left(x\right)=c_{1,\mu }\mathrm{P}.\mathrm{V}.{\int }_{\mathbb{R}}{\displaystyle \frac{\varphi \left(x\right)-\varphi \left(t\right)}{{|x-t|}^{1+\mu }}}dt,\mu \in (1,2),\mathrm{for}\mathrm{all}\varphi \in \mathcal{S},\end{array}$

where $\mathcal{S}$ denotes the Schwartz space of rapidly decreasing functions on $\mathbb{R}$, and $c_{1,\mu }$ is the normalization constant, which can be expressed as follows:

$\begin{array}{c}\hfill c_{1,\mu }={\displaystyle \frac{2^{\mu -1}\mu \mathsf{\Gamma }\left(\frac{\mu +1}{2}\right)}{\sqrt{\pi }\mathsf{\Gamma }(1-\mu /2)}}.\end{array}$

The principal value integral (P.V.) is defined as in the following equation:

(3)$\mathrm{P}.\mathrm{V}.{\int }_{\mathbb{R}}{\displaystyle \frac{\varphi \left(x\right)-\varphi \left(t\right)}{{\left|x-t\right|}^{1+\mu }}}dt=\underset{ϵ\to 0}{lim}{\int }_{\mathbb{R}\setminus B_{ϵ}\left(x\right)}{\displaystyle \frac{\varphi \left(x\right)-\varphi \left(t\right)}{{\left|x-t\right|}^{1+\mu }}}dt.$

The spatial fractional Laplacian operator ${(-\Delta )}^{\mu /2}$ governs the anomalous diffusion mechanism, capturing long-range interactions through its nonlocal formulation, where particle displacements follow Lévy flights rather than Brownian motion. The temporal Caputo derivative ${}_0{}^{C}D_t^{\beta }$ introduces memory effects into the system, modeling sub-diffusive processes where particles experience waiting times between jumps with heavy-tailed distributions. The source term $f(x,t)$ represents external influences on the system, such as the localized injection/removal of material, thermal energy sources in heat transfer problems, or biological growth/death rates in population dynamics. The linear reaction term $\lambda \varphi (x,t)$ characterizes the intrinsic system behavior; positive $\lambda$ values model the autocatalytic growth or amplification of feedback mechanisms, while negative values describe decay processes like radioactive disintegration or damping effects. This combination of fractional operators extends classical reaction–diffusion frameworks to systems exhibiting simultaneous non-Markovian temporal evolution and nonlocal spatial interactions, making it particularly relevant for modeling complex phenomena in viscoelastic materials, heterogeneous media, and biological systems with memory-dependent transport characteristics [8,9,10].

Challenges in the numerical solution of the time-space fractional reaction–diffusion equation (Equation (1)) stem mainly from the singular nature of its solution near boundaries. The most recent regularity findings in this regard were presented in [11,12,13]. A numerical analysis of the reaction–diffusion equation with a hyper-singular integral fractional Laplacian, with or without the Caputo time-fractional operator, has only been considered by a few authors. In [14], implicit difference schemes for solving the time-space fractional diffusion equation with integral fractional Laplacian were introduced using the graded L1 formula and special finite difference discretization. A fractional-centered difference approach was considered for the fractional diffusion equation with a high-dimensional hyper-singular integral fractional Laplacian [15]. In [16], the finite element method on a locally refined composite mesh was applied to discretize the fractional Poisson equation, the fractional-in-space Allen–Cahn equation, and the fractional Burgers equation with an integral fractional Laplacian. A numerical method with a numerical quadrature based on the singular integral representation for the fractional Laplacian was proposed in [17]. In [18], a central difference technique was utilized in conjunction with a “punched-hole” trapezoidal formula to approximate the fractional Laplacian. This method was used to study stochastic differential equations with non-Gaussian $\mu$-stable Lévy motions. The integral version of the Dirichlet homogeneous fractional Laplace equation was discussed in [19], providing weighted and fractional Sobolev estimates and proving optimal convergence for the standard linear finite element method derived for quasi-uniform and graded meshes. The study in [20] focused on constructing the components needed to build an adaptive finite element code for the integral fractional Laplacian. The goal of the study [21] was to obtain convergence results for the approximation of solutions to the fractional Laplacian equation, considering both domain truncation and spatial discretization. Mao and Shen [22] developed efficient Hermite–collocation and Hermite–Galerkin methods to solve a class of fractional PDEs in unbounded domains directly, deriving corresponding error estimates. Tang et al. [23] developed accurate spectral methods using rational basis functions for fractional partial differential equations involving fractional Laplacian in unbounded domains.

For problems involving non-smooth solutions in both time and space, hybrid finite difference/spectral approaches can face challenges, as these methods often rely on the assumption of smoothness in at least one dimension [24,25,26,27]. In cases where smoothness is present, spectral methods with smooth basis functions such as Legendre polynomials are typically employed for spatial discretization, while finite difference schemes on non-uniform graded meshes are used for time discretization [28]. However, for non-smooth solutions, these techniques may need adaptation or alternative strategies, as the lack of regularity in the spatial direction can degrade their accuracy and efficiency. To the best of our knowledge, there is essentially no work available along this line.

In this paper, we intend to fill this gap by designing and analyzing a hybrid finite difference/Pertrov–Galerkin spectral approach for a class of time-fractional reaction–diffusion equations with a spatial hyper-singular integral fractional Laplacian. We apply L1 discretization to the temporal fractional derivative, where the graded mesh can capture the model problem with weak singularity at the initial time. Moreover, the Generalized Jacobi Pertrov–Galerkin method is employed in the spatial direction to handle non-smooth spatial solutions effectively.

The remainder of this paper is organized as follows: In the next section, we present the weighted Sobolev spaces. Section 2 introduces the time discretization of the fractional reaction–diffusion equation. Section 4 discusses the stability analysis. Section 5 proves the convergence analysis of the fully discrete scheme. Section 6 describes the numerical experiments used to support our theoretical results. Finally, Section 7 summarizes the main conclusions.

2. Weighted Sobolev Spaces

Let ${\omega }^{{\mu }_1,{\mu }_2}={(1-x)}^{{\mu }_1}{(1+x)}^{{\mu }_2}$ with ${\mu }_1,{\mu }_2>-1$. Define the inner product of $L_{{\omega }^{{\mu }_1,{\mu }_2}}^2\left(I\right)$ as follows:

(4)${(u,\nu )}_{{\omega }^{{\mu }_1,{\mu }_2}}={\int }_{I}u\nu {\omega }^{{\mu }_1,{\mu }_2}dx,$

which is equipped with the norm ${\parallel u\parallel }_{{\omega }^{{\mu }_1,{\mu }_2}}$. Subscripts may be omitted when ${\omega }^{{\mu }_1,{\mu }_2}=1$, if no confusion arises.

To account for endpoint singularities, we introduce the weighted Sobolev space (cf. [29]) as follows:

(5)$B_{{\omega }^{{\mu }_1,{\mu }_2}}^q\left(I\right)=\left\{\varphi |{\partial }_x^j\varphi \in L_{{\omega }^{{\mu }_1+j,{\mu }_2+j}}^2\left(I\right),j=0,1,\dots ,q\right\},q\mathrm{a}\mathrm{nonnegative}\mathrm{integer},$

with the following norm:

(6)${\parallel \varphi \parallel }_{B_{{\omega }^{{\mu }_1,{\mu }_2}}^q}={\left({\displaystyle \sum _{i=0}^q}{\left|\varphi \right|}_{B_{{\omega }^{{\mu }_1,{\mu }_2}}^i}^2\right)}^{1/2}{,\left|\varphi \right|}_{B_{{\omega }^{{\mu }_1,{\mu }_2}}^i}={\parallel {\partial }_x^i\varphi \parallel }_{{\omega }^{{\mu }_1+i,{\mu }_2+i}}.$

These spaces are closely related to Jacobi polynomials. The orthogonality of Jacobi polynomials $P_n^{{\mu }_1,{\mu }_2}\left(x\right)$ is given by the following:

(7)${\left(P_r^{{\mu }_1,{\mu }_2}\left(x\right),P_s^{{\mu }_1,{\mu }_2}\left(x\right)\right)}_{{\omega }^{{\mu }_1,{\mu }_2}}=h_n^{{\mu }_1,{\mu }_2}{\delta }_{rs},$

where ${\delta }_{rs}$ is the Kronecker delta and:

(8)$h_s^{{\mu }_1,{\mu }_2}={\parallel P_s^{{\mu }_1,{\mu }_2}\parallel }_{{\omega }^{{\mu }_1,{\mu }_2}}={\displaystyle \frac{2^{{\mu }_1+{\mu }_2+1}}{2s+{\mu }_1+{\mu }_2+1}}{\displaystyle \frac{\mathsf{\Gamma }(s+{\mu }_1+1)\mathsf{\Gamma }(s+{\mu }_2+1)}{\mathsf{\Gamma }(s+{\mu }_1+{\mu }_2+1)\mathsf{\Gamma }(s+1)}}.$

For notational simplicity, we omit the second index when ${\mu }_1={\mu }_2$. For instance, $P_n^{\mu /2}\left(x\right)$ denotes $P_n^{\mu /2,\mu /2}\left(x\right)$, ${\omega }^{\mu /2}\left(x\right)$ represents ${\omega }^{\mu /2,\mu /2}\left(x\right)$, and $h_n^{\mu /2}$ stands for $h_n^{\mu /2,\mu /2}$. When $q=\gamma$ is a non-integer, the space is defined via interpolation (see [30]). For any $\gamma$, the norm ${\parallel \varphi \parallel }_{B_{{\omega }^{{\mu }_1,{\mu }_2}}^{\gamma }}$ is equivalent to (cf. [30]) the following:

(9)${⦀\varphi ⦀}_{B_{{\omega }^{{\mu }_1,{\mu }_2}}^{\gamma }}^2={\displaystyle \sum _{n=0}^{\infty }}{\left({\varphi }_n^{{\mu }_1,{\mu }_2}\right)}^2{\displaystyle \frac{1+n^{2\gamma }}{2n+{\mu }_1+{\mu }_2+1}},\forall \varphi \in B_{{\omega }^{{\mu }_1,{\mu }_2}}^{\gamma },$

where ${\varphi }_n^{{\mu }_1,{\mu }_2}$ are the Jacobi–Fourier coefficients of $\varphi$ with respect to $P_n^{\mu /2}\left(x\right)$. The equivalence of norms allows for their interchangeable use. These spaces coincide with Besov spaces [30].

An equivalent characterization of the weighted Sobolev space is provided in [31]. For $B_{\mu }^s\left(I\right)$ with $s=m+\sigma$, $0<\sigma <1$, and $s\ne 1+\mu$ if $-1<\mu <0$, the norm is equivalent to the following:

(10)${\parallel \varphi \parallel }_{B_{{\omega }^{\mu }}^s}={\left({\parallel \varphi \parallel }_{B_{{\omega }^{\mu }}^m}^2+{\int \int }_{{\mathsf{\Omega }}_{I,a}}{\omega }^{\mu +s}{\displaystyle \frac{|{\varphi }^{\left(m\right)}\left(x\right)-{\varphi }^{\left(m\right)}{\left(y\right)|}^2}{{|x-y|}^{1+2\sigma }}}dxdy\right)}^{1/2},$

where ${\mathsf{\Omega }}_{I,a}\subset I\times I$ for $a>1$ is defined as follows:

(11)${\mathsf{\Omega }}_{I,a}=\left\{(x,y)\in I\times I|a^{-1}(1-|x\left|\right)<1-\mathrm{sgn}\left(x\right)y<a(1-|x\left|\right)\right\}.$

Here, $a>1$ is arbitrary; we fix $a=2$.

The unweighted Sobolev space $H^{q,s}\left(I\right)$ is equipped with the norm, as follows:

(12)${\parallel \varphi \parallel }_{H^{q,s}\left(I\right)}={\left({\parallel \varphi \parallel }_{H^q}^2+\sum _{|{\mu }_2|=q}{\left|D^{{\mu }_2}\varphi \right|}_{H^s}^2\right)}^{1/2},{\left|\varphi \right|}_{H^s}={\left({\int }_I{\int }_I{\displaystyle \frac{{|\varphi \left(x\right)-\varphi \left(y\right)|}^2}{{|x-y|}^{1+2s}}}dxdy\right)}^{1/2},$

where $q\ge 0$ is an integer, $0<s<1$, and:

$H^q\left(I\right)=\left\{\varphi |\varphi ,{\partial }_x^l\varphi \in L^2\left(I\right),l=1,\dots ,q\right\}.$

Frequently, we use the following notation:

${\left|\varphi \right|}_{H_{\mu }^s}={\left({\int }_I{\int }_I{\displaystyle \frac{{|\varphi \left(x\right)-\varphi \left(y\right)|}^2{\omega }^{\mu }\left(x\right)}{{|x-y|}^{1+2s}}}dxdy\right)}^{1/2}.$

3. Numerical Scheme and Implementation

This section presents the temporal discretization for the fractional reaction–diffusion equation—Equation (1)—using Caputo derivative approximation.

3.1. Temporal Discretization on a Graded Mesh

Let K be a positive integer. Define temporal nodes $t_i=T{(i/K)}^{\rho }$ for $i=0,1,\dots ,K$ with mesh grading parameter $\rho \ge 1$. A uniform mesh corresponds to $\rho =1$. The time step sizes ${\tau }_i:=t_i-t_{i-1}$ decrease near $t=0$ to resolve solution singularities. Let ${\varphi }^i:=\varphi (·,t_i)$ denote the numerical approximation at $t_i$.

The non-uniform L1 scheme approximates the Caputo derivative ${}_0{}^{C}D_t^{\beta }\varphi (x,t_n)$ as follows:

(13)$\begin{array}{cc}\hfill {\left.{}_0{}^{C}D_t^{\beta }\varphi (x,t)\right|}_{t=t_n}& =D_K^{\beta }{\varphi }^n+R^n\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^n}{\tilde{a}}_{n-i}^n({\varphi }^i-{\varphi }^{i-1})+R^n\hfill \\ \hfill & =d_{n,1}{\varphi }^n-d_{n,n}{\varphi }^0-{\displaystyle \sum _{i=1}^{n-1}}(d_{n,n-i}-d_{n,n-i+1}){\varphi }^i+R^n,\hfill \end{array}$

where coefficients satisfy the following:

$d_{n,n-i}={\displaystyle \frac{{(t_n-t_i)}^{1-\beta }-{(t_n-t_{i+1})}^{1-\beta }}{{\tau }_{i+1}\mathsf{\Gamma }(2-\beta )}},0\le i\le n-1$

and ${\omega }_{\beta }\left(t\right)=t^{\beta -1}/\mathsf{\Gamma }\left(\beta \right)$. Equivalently:

(14)$D_K^{\beta }{\varphi }^n={\displaystyle \sum _{i=0}^n}{\tilde{b}}_{n-i}^n{\varphi }^i,n=1,\dots ,K$

with ${\tilde{b}}_0^n=d_{n,1}$, ${\tilde{b}}_n^n=-d_{n,n}$, and ${\tilde{b}}_{n-i}^n=d_{n,n-i}-d_{n,n-i+1}$ for $1\le i\le n-1$. These coefficients satisfy the following:

(15)$d_{n,n-i}>d_{n,n-i+1}>0,0\le i\le n-1.$

Lemma 1

([32]). If $\varphi \in C^2\left((0,T]\right)$ satisfies $|{\varphi }^{″}\left(t\right)|\le C_{\varphi }(1+t^{\beta -2})$ for $0<t\le T$ with non-decreasing time steps ${\tau }_{n-1}\le {\tau }_n$ ($2\le n\le K$), then:

(16)$\left|{}_0{}^{C}D_t^{\beta }\varphi \left(t_n\right)-D_K^{\beta }\varphi \left(t_n\right)\right|\le C_{\varphi }{K}^{-min\{\rho \beta ,2-\beta \}}.$

Define positive multipliers ${\mu }_{m,j}$ recursively as ${\mu }_{m,m}:=1$ and:

${\mu }_{m,j}:={\displaystyle \sum _{k=1}^{m-j}}{\displaystyle \frac{1}{d_{m-k,1}}}(d_{m,k}-d_{m,k+1}){\mu }_{m-k,j}$

where $m=1,2,\dots ,K$ and $1\le j\le m-1$.

Lemma 2

(Corollary 5.4 [33]). For $m=1,2,\dots ,K$:

(17)$\mathsf{\Gamma }(2-\beta ){\tau }_m^{\beta }{\displaystyle \sum _{j=1}^m}{\mu }_{m,j}\le {\displaystyle \frac{t_m^{\beta }}{\mathsf{\Gamma }(1+\beta )}}.$

3.2. The L1–Pertrov–Galerkin Spectral Scheme

Key results for matrix construction:

Lemma 3

(see [34]). For $i\ge 0$ and $x\in I$:

(18)${(-\Delta )}^{\mu /2}\left[{(1-x^2)}^{\mu /2}{P}_i^{\mu /2}\left(x\right)\right]=d_i^{\mu }{P}_i^{\mu /2}\left(x\right),d_i^{\mu }={\displaystyle \frac{\mathsf{\Gamma }(\mu +i+1)}{i!}}.$

Lemma 4

(see [35]). Jacobi polynomial expansion:

(19)$P_n^{{\mu }_1,{\mu }_2}\left(x\right)={\displaystyle \sum _{k=0}^n}{\widehat{c}}_k^n({\mu }_1,{\mu }_2,a,b)P_k^{a,b}\left(x\right)$

with coefficients:

(20)$\begin{array}{cc}\hfill {\widehat{c}}_k^n({\mu }_1,{\mu }_2,a,b)& ={\displaystyle \frac{(2k+a+b+1)\mathsf{\Gamma }\left(n+{\mu }_1+1\right)\mathsf{\Gamma }(k+a+b+1)}{\mathsf{\Gamma }(n+{\mu }_1+{\mu }_2+1)\mathsf{\Gamma }\left(n+a+1\right)}}\hfill \\ \hfill & \times {\displaystyle \sum _{m=0}^{n-k}}{\displaystyle \frac{{(-1)}^m\mathsf{\Gamma }(n+k+m+{\mu }_1+{\mu }_2+1)\mathsf{\Gamma }(m+k+a+1)}{m!(n-k-m)!\mathsf{\Gamma }(k+m+{\mu }_1+1)\mathsf{\Gamma }(m+2k+a+b+1)}}.\hfill \end{array}$

Define the approximation space as follows:

(21)${\omega }^{\mu /2}{\mathbb{P}}_N=\left\{g\right|g\left(x\right)={(1-x^2)}^{\mu /2}\psi \left(x\right),\psi \in {\mathbb{P}}_N\}.$

The L1–Petrov–Galerkin method seeks ${\varphi }_N^n\in {\omega }^{\mu /2}{\mathbb{P}}_N$ satisfying the following:

(22)$(D_K^{\beta }{\varphi }_N^n,v)+({(-\Delta )}^{\mu /2}{\varphi }_N^n,v)+\lambda ({\varphi }_N^n,v)=(f^n,v)\forall v\in {\mathbb{P}}_N.$

Rewriting yields the following:

(23)${\overline{b}}_0^n({\varphi }_N^n,v)+({(-\Delta )}^{\mu /2}{\varphi }_N^n,v)={\displaystyle \sum _{i=1}^{n-1}}{\tilde{b}}_{n-i}^n({\varphi }_N^i,v)+({\overline{f}}^n,v),$

where ${\overline{b}}_0^n={\tilde{b}}_0^n+\lambda$ and ${\overline{f}}^n=f^n+{\tilde{b}}_n^n{g}_0$.

Express the solution as follows:

(24)${\varphi }_N^n\left(x\right)={\displaystyle \sum _{i=0}^N}{\widehat{\varphi }}_i^n{\phi }_i\left(x\right)\in {\omega }^{\mu /2}{\mathbb{P}}_N,{\phi }_i\left(x\right)={(1-x^2)}^{\mu /2}{P}_i^{\mu /2}\left(x\right).$

Substituting into (23) with test functions $v=P_j^{\mu /2}\left(x\right)$ yields the following:

(25)$({\overline{b}}_0^n\overline{M}+\overline{S})U^n={\overline{K}}^{n-1}+{\overline{F}}^n$

Using Lemma 3 and 4, we obtain the following:

(26)$\begin{array}{cc}\hfill & {\overline{S}}_{i,j}=\left({\left(-\Delta \right)}^{\mu /2}{\phi }_i,P_j^{\mu /2}\right)=d_i^{\mu }\left(P_i^{\mu /2},P_j^{\mu /2}\right),{\overline{M}}_{i,j}=\left({\phi }_i,P_j^{\mu /2}\right)=h_i^{\mu /2}{\delta }_{i,j},\hfill \\ \hfill & {\overline{K}}^{n-1}={\displaystyle \sum _{i=1}^{n-1}}{\tilde{b}}_{n-i}^n\overline{M}{U}^i,{\overline{F}}_i^n=\left({\overline{f}}^n,{\phi }_i\right)\hfill \end{array}$

There are two methods that can be used for calculating the matrix $\overline{S}$.

(27)${\overline{S}}_{i,j}=d_i^{\mu }{\int }_{-1}^1{P}_i^{\mu /2}\left(x\right)P_j^{\mu /2}\left(x\right)dx.$

First, we precisely express the integral given above using the orthogonality of Legendre polynomials and Lemma 4 as follows:

(28)${\overline{S}}_{i,j}=d_i^{\mu }{\displaystyle \sum _{r=0}^{min\left\{i,j\right\}}}{\widehat{c}}_r^i(\mu /2,\mu /2,0,0){\widehat{c}}_r^j(\mu /2,\mu /2,0,0)h_r^0.$

The Gauss–Legendre quadrature is another tool for calculating ${\overline{S}}_{i,j}$, which can be expressed as follows:

(29)${\overline{S}}_{i,j}=d_i^{\mu }{\displaystyle \sum _{r=0}^N}{P}_j^{\mu /2}\left(x_r^0\right)P_i^{\mu /2}\left(x_r^0\right){\varpi }_r^0,$

where $\left\{{\varpi }_i^0,x_i^0\right\}$ are the Gauss–Legendre weights and nodes. As $i+j\le 2N$, the above quadrature rule is exact, where the Gauss quadrature rule is exact for all polynomials of order $(2N+1)$. To compute ${\overline{F}}_i^n=\left({\overline{f}}^n,{\phi }_i\right)$, we apply a suitable Jacobi–Gauss quadrature rule.

4. Stability Analysis

Lemma 5

(Lemma 4.2 [36]). Let ${\left\{{\xi }^m\right\}}_{m=1}^{\infty }$ and ${\left\{{\eta }^m\right\}}_{m=1}^{\infty }$ be nonnegative sequences. If the nonnegative grid function ${\left\{{\psi }^m\right\}}_{m=0}^K$ satisfies ${\psi }^0\ge 0$ and

$\left(D_K^{\beta }{\psi }^m\right){\psi }^m\le {\xi }^m{\psi }^m+{\left({\eta }^m\right)}^2\mathrm{for}m=1,2,\dots ,K,$

then for any $1\le m\le K$:

${\psi }^m\le {\psi }^0+\mathsf{\Gamma }(2-\beta ){\tau }_m^{\beta }{\displaystyle \sum _{j=1}^m}{\mu }_{m,j}({\xi }^j+{\eta }^j)+\underset{1\le j\le m}{max}\left\{{\eta }^j\right\}.$

Lemma 6.

Let u and w vanish outside $\mathsf{\Omega }\subseteq {\mathbb{R}}^d$ almost everywhere, with w such that all integrals below exist. Define $\rho \left(x\right)=c_{d,\mu }{\int }_{{\mathsf{\Omega }}^c}{\displaystyle \frac{1}{{|x-y|}^{d+\mu }}}dy$. Then:

(30) $\begin{array}{cc}\hfill 2{\int }_{\mathsf{\Omega }}wu{(-\Delta )}^{\mu /2}\varphi \left(x\right)dx=& c_{d,\mu }{\int }_{\mathsf{\Omega }}{\int }_{\mathsf{\Omega }}{\displaystyle \frac{{(\varphi \left(x\right)-\varphi \left(y\right))}^{2}w\left(y\right)}{{\left|x-y\right|}^{d+\mu }}}dydx\hfill \\ & +{\int }_{\mathsf{\Omega }}{u}^2\left(x\right)\left(\rho \left(x\right)w\left(x\right)+{(-\Delta )}^{\mu /2}w\right)dx,\hfill \end{array}$

provided all integrals are well-defined.

Theorem 1.

Let ${\varphi }_N\in {\omega }^{\mu /2}{\mathbb{P}}_N$ be the solution of (22). Then:

(31) ${∥{\varphi }_N^n∥}_{{\omega }^{-\mu /2}}\le {∥{\varphi }_N^0∥}_{{\omega }^{-\mu /2}}+{\displaystyle \frac{t_n^{\beta }}{\mathsf{\Gamma }(1+\beta )}}{∥f^j∥}_{{\omega }^{-\mu /2}}^2,\forall 1\le n\le K.$

Proof. 

Taking $\nu ={\tilde{\varphi }}_N^n={(1-x^2)}^{-\mu /2}{\varphi }_N^n$ in (22) yields the following:

(32)${(D_K^{\beta }{\tilde{\varphi }}_N^n,{\tilde{\varphi }}_N^n)}_{{\omega }^{\mu /2}}+({(-\Delta )}^{\mu /2}({\omega }^{\mu /2}{\tilde{\varphi }}_N^n,{\tilde{\varphi }}_N^n)+\lambda {({\tilde{\varphi }}_N^n,{\tilde{\varphi }}_N^n)}_{{\omega }^{\mu /2}}=(f^n,{\tilde{\varphi }}_N^n).$

Using Lemma 6 and the identity ${(-\Delta )}^{\mu /2}{\omega }^{\mu /2}=\mathsf{\Gamma }(\mu +1)$, we obtain the following:

(33)$2({(-\Delta )}^{\mu /2}\left({\omega }^{\mu /2}{\tilde{\varphi }}_N\right),{\tilde{\varphi }}_N)=c_{1,\mu }|{\tilde{\varphi }}_N{|}_{H_{\mu /2}^{\mu /2}}^2+\mathsf{\Gamma }(\mu +1){\parallel {\tilde{\varphi }}_N\parallel }^2+{\parallel {\tilde{\varphi }}_N{\rho }^{1/2}\parallel }_{{\omega }^{\mu /2}}^2.$

From the weighted norm definition in (10), we derive the following:

$|{\tilde{\varphi }}_N{|}_{B_{{\omega }^{\mu /2}}^{\mu /2}}\le {\left|{\tilde{\varphi }}_N\right|}_{H_{\mu /2}^{\mu /2}}.$

Thus:

$2({(-\Delta )}^{\mu /2}\left({\omega }^{\mu /2}{\tilde{\varphi }}_N\right),{\tilde{\varphi }}_N)\ge c_{1,\mu }|{\tilde{\varphi }}_N{|}_{B_{{\omega }^{\mu /2}}^{\mu /2}}^2+\mathsf{\Gamma }(\mu +1){\parallel {\tilde{\varphi }}_N\parallel }^2+{\parallel {\tilde{\varphi }}_N{\rho }^{1/2}\parallel }_{{\omega }^{\mu /2}}^2.$

Applying Cauchy–Schwarz inequality:

(34)$\begin{array}{cc}\hfill 2{(D_K^{\beta }{\tilde{\varphi }}_N^n,{\tilde{\varphi }}_N^n)}_{{\omega }^{\mu /2}}& + c_{1,\mu }|{\tilde{\varphi }}_N{|}_{B_{{\omega }^{\mu /2}}^{\mu /2}}^2+\mathsf{\Gamma }(\mu +1){\parallel {\tilde{\varphi }}_N\parallel }^2\hfill \\ \hfill & + {\parallel {\tilde{\varphi }}_N{\rho }^{1/2}\parallel }_{{\omega }^{\mu /2}}^2+2\lambda {∥{\tilde{\varphi }}_N^n∥}_{{\omega }^{\mu /2}}^2\le 2(f^n,{\tilde{\varphi }}_N^n).\hfill \end{array}$

This simplifies to the following:

(35)${(D_K^{\beta }{\tilde{\varphi }}_N^n,{\tilde{\varphi }}_N^n)}_{{\omega }^{\mu /2}}\le (f^n,{\tilde{\varphi }}_N^n).$

Through coefficient manipulation:

$\begin{array}{cc}\hfill d_{m,1}{\parallel {\tilde{\varphi }}_N^m\parallel }_{{\omega }^{\mu /2}}^2& \le {\displaystyle \sum _{i=0}^{m-1}}(d_{m,i}-d_{m,i+1}){({\tilde{\varphi }}_N^i,{\tilde{\varphi }}_N^m)}_{{\omega }^{\mu /2}}\hfill \\ \hfill & + d_{m,m}{({\tilde{\varphi }}_N^0,{\tilde{\varphi }}_N^m)}_{{\omega }^{\mu /2}}+(f^m,{\tilde{\varphi }}_N^m).\hfill \end{array}$

Applying Cauchy–Schwarz again:

$\begin{array}{cc}\hfill d_{m,1}{\parallel {\tilde{\varphi }}_N^m\parallel }_{{\omega }^{\mu /2}}^2& \le {\displaystyle \sum _{i=0}^{m-1}}(d_{m,m-i}-d_{m,m-i+1}){\parallel {\tilde{\varphi }}_N^i\parallel }_{{\omega }^{\mu /2}}{\parallel {\tilde{\varphi }}_N^m\parallel }_{{\omega }^{\mu /2}}\hfill \\ \hfill & + d_{m,m}{\parallel {\tilde{\varphi }}_N^0\parallel }_{{\omega }^{\mu /2}}{\parallel {\tilde{\varphi }}_N^m\parallel }_{{\omega }^{\mu /2}}+{\parallel f^m\parallel }_{{\omega }^{-\mu /2}}{\parallel {\tilde{\varphi }}_N^m\parallel }_{{\omega }^{\mu /2}}.\hfill \end{array}$

Thus:

$\left(D_K^{\mu }{\parallel {\tilde{\varphi }}_N^m\parallel }_{{\omega }^{\mu /2}}\right){\parallel {\tilde{\varphi }}_N^m\parallel }_{{\omega }^{\mu /2}}\le {\parallel f^m\parallel }_{{\omega }^{-\mu /2}}{\parallel {\tilde{\varphi }}_N^m\parallel }_{{\omega }^{\mu /2}}.$

Lemmas 2 and 5 yield the following:

$\begin{array}{cc}\hfill {\parallel {\varphi }_N^m\parallel }_{{\omega }^{-\mu /2}}& \le {\parallel {\varphi }_N^0\parallel }_{{\omega }^{-\mu /2}}+\mathsf{\Gamma }(2-\beta ){\tau }_m^{\beta }{\displaystyle \sum _{j=1}^m}{\mu }_{m,j}{\parallel f^j\parallel }_{{\omega }^{-\mu /2}}\hfill \\ \hfill & \le {\parallel {\varphi }_N^0\parallel }_{{\omega }^{-\mu /2}}+{\displaystyle \frac{t_m^{\beta }}{\mathsf{\Gamma }(1+\beta )}}\underset{1\le j\le m}{max}{\parallel f^j\parallel }_{{\omega }^{-\mu /2}}.\hfill \end{array}$

□

Remark 1.

The stability result in Theorem 1 is more accurate and β-robust, i.e., the bound does not blow up when $\beta \to 1^{-}$. In contrast, the bound in ([37], Theorem 4.1) is not β-robust, as there exists a factor $\mathsf{\Gamma }(1-\beta )$ that will blow up when $\beta \to 1^{-}$.

5. Convergence Analysis

Let the projection operator $P_N^{{\mu }_1,{\mu }_2}:L_{{\omega }^{{\mu }_1,{\mu }_2}}^2\left(I\right)\to {\mathbb{P}}_N$ be defined by ${(P_N^{{\mu }_1,{\mu }_2}\varphi -\varphi ,\nu )}_{{\omega }^{{\mu }_1,{\mu }_2}}=0$ for all $\varphi \in L_{{\omega }^{{\mu }_1,{\mu }_2}}^2\left(I\right)$ and $\nu \in {\mathbb{P}}_N$. This projection error satisfies the following property [38]:

Lemma 7

([39]). For any $\varphi \in B_{{\omega }^{{\mu }_1,{\mu }_2}}^r\left(I\right)$ and $0\le r_1\le r$,

(36)${∥P_N^{{\mu }_1,{\mu }_2}\varphi -\varphi ∥}_{B_{{\omega }^{{\mu }_1,{\mu }_2}}^{r_1}}\le C{\left(N(N+{\mu }_1+{\mu }_2)\right)}^{{\displaystyle \frac{r_1-r}{2}}}{\left|\varphi \right|}_{B_{{\omega }^{{\mu }_1,{\mu }_2}}^r},$

where $C>0$ is independent of ϕ, N, ${\mu }_1$, and ${\mu }_2$.

Theorem 2.

Let $\tilde{\varphi }={\omega }^{-\mu /2}\varphi$ and ${\tilde{\varphi }}_N={\omega }^{-\mu /2}{\varphi }_N$. Assume $f\in L^{\infty }(0,T,B_{{\omega }^{\mu /2-1}}^r\left(I\right))$ with $r\ge 0$, and that ${}_0{}^{C}D_t^{\mu }\tilde{\varphi }$ and $\tilde{\varphi }\in L^{\infty }(0,T,B_{{\omega }^{\mu /2-1}}^{r\wedge (\mu -ϵ)+\mu -ϵ}\left(I\right)\cap B_{{\omega }^{\mu /2}}^{\gamma }\left(I\right))$ for $\gamma =(r\wedge (\mu -ϵ)+\mu -ϵ)\wedge (\mu +1)\wedge r+\mu$, where $ϵ>0$ is arbitrarily small. Then:

(37) ${∥\varphi -{\varphi }_N^n∥}_{{\omega }^{-\mu /2}}\le C\left(N^{\mu -\gamma }+K^{-min\{\rho \beta ,2-\beta \}}\right),$

where $C>0$ is independent of K and N.

Proof. 

For $m=0,\dots ,K$, define $e_N^m={\varphi }_N^m-{\varphi }_N^{*m}\left(t_m\right)={\omega }^{\mu /2}{\tilde{e}}_N^n$ and ${\eta }_N^m=\varphi \left(t_m\right)-{\varphi }_N^{*m}\left(t_m\right)={\omega }^{\mu /2}{\tilde{\eta }}_N^m$, where ${\varphi }_N^{*}$ will be specified later.

From (1) and (22), we derive the error equation for $n=1,\dots ,K$:

$\begin{array}{c}\hfill \left(D_K^{\beta }{e}_N^n+{(-\Delta )}^{\mu /2}{e}_N^n+\lambda e_N^n,v_N\right)=\left(R^n+D_K^{\beta }{\eta }_N^n+{(-\Delta )}^{\mu /2}{\eta }_N^n+\lambda {\eta }_N^n,v_N\right),\end{array}$

where $R^n={}_0{}^{C}D_t^{\beta }{\varphi }^n-D_K^{\beta }{\varphi }^n$. Choosing $v_N={\tilde{e}}_N^n$ and using (35):

(38)$\begin{array}{cc}\hfill 2{(D_K^{\beta }{\tilde{e}}_N^n,{\tilde{e}}_N^n)}_{{\omega }^{\mu /2}}& + c_{1,\mu }|{\tilde{e}}_N^n{|}_{B_{{\omega }^{\mu /2}}^{\mu /2}}^2+\mathsf{\Gamma }(\mu +1){\parallel {\tilde{e}}_N^n\parallel }^2+{\parallel {\tilde{e}}_N^n{\rho }^{1/2}\parallel }_{{\omega }^{\mu /2}}^2\hfill \\ \hfill & + 2\lambda {∥{\tilde{e}}_N^n∥}_{{\omega }^{\mu /2}}^2\le 2\left(D_K^{\beta }{\eta }_N^n+{(-\Delta )}^{\mu /2}{\eta }_N^n+\lambda {\eta }_N^n+R^n,{\tilde{e}}_N^n\right).\hfill \end{array}$

Using ${\omega }^{-\mu }\le \rho \left(x\right)$ and $2{(D_K^{\beta }{\tilde{e}}_N^n,{\tilde{e}}_N^n)}_{{\omega }^{\mu /2}}\ge D_K^{\beta }{∥{\tilde{e}}_N^n∥}_{{\omega }^{\mu /2}}^2$:

(39)$\begin{array}{cc}\hfill D_K^{\beta }{∥{\tilde{e}}_N^n∥}_{{\omega }^{\mu /2}}^2 +& c_{1,\mu }{\left|{\tilde{e}}_N^n\right|}_{B_{{\omega }^{\mu /2}}^{\mu /2}}^2+\mathsf{\Gamma }(\mu +1){∥{\tilde{e}}_N^n∥}^2+{∥{\tilde{e}}_N^n∥}_{{\omega }^{-\mu /2}}^2+2\lambda {∥{\tilde{e}}_N^n∥}_{{\omega }^{\mu /2}}^2\hfill \\ \hfill & \le {∥{(-\Delta )}^{\mu /2}{\eta }_N^n)∥}_{{\omega }^{\mu /2}}^2+{∥{\tilde{e}}_N^n∥}_{{\omega }^{-\mu /2}}^2+2\left(D_K^{\beta }{\eta }_N^n+\lambda {\eta }_N^n+R^n,{\tilde{e}}_N^n\right).\hfill \end{array}$

(40)$\begin{array}{cc}\hfill D_K^{\beta }{∥{\tilde{e}}_N^n∥}_{{\omega }^{\mu /2}}^2 +& c_{1,\mu }{\left|{\tilde{e}}_N^n\right|}_{B_{{\omega }^{\mu /2}}^{\mu /2}}^2+\mathsf{\Gamma }(\mu +1){∥{\tilde{e}}_N^n∥}^2+2\lambda {∥{\tilde{e}}_N^n∥}_{{\omega }^{\mu /2}}^2\hfill \\ \hfill & \le {∥{(-\Delta )}^{\mu /2}{\eta }_N^n∥}_{{\omega }^{\mu /2}}^2+2\left(D_K^{\beta }{\eta }_N^n+\lambda {\eta }_N^n+R^n,{\tilde{e}}_N^n\right)\hfill \\ \hfill & \le {∥{(-\Delta )}^{\mu /2}{\eta }_N^n∥}_{{\omega }^{\mu /2}}^2+2\left(∥D_K^{\beta }{\eta }_N^n∥+\lambda ∥{\eta }_N^n∥+∥R^n∥\right)∥{\tilde{e}}_N^n∥,\hfill \end{array}$

(41)$\begin{array}{cc}\hfill D_K^{\beta }{∥{\tilde{e}}_N^n∥}_{{\omega }^{\mu /2}}^2 +& c_{1,\mu }{\left|{\tilde{e}}_N^n\right|}_{B_{{\omega }^{\mu /2}}^{\mu /2}}^2+\left(\mathsf{\Gamma }(\mu +1)-1\right){∥{\tilde{e}}_N^n∥}^2+2\lambda {∥{\tilde{e}}_N^n∥}_{{\omega }^{\mu /2}}^2\hfill \\ \hfill & \le {∥{(-\Delta )}^{\mu /2}{\eta }_N^n∥}_{{\omega }^{\mu /2}}^2+{\left(∥D_K^{\beta }{\eta }_N^n∥+\lambda ∥{\eta }_N^n∥+∥R^n∥\right)}^2,\hfill \end{array}$

(42)$\begin{array}{cc}\hfill D_K^{\beta }{∥{\tilde{e}}_N^n∥}_{{\omega }^{\mu /2}}^2& \le C{\left({∥{(-\Delta )}^{\mu /2}{\eta }_N^n∥}_{{\omega }^{\mu /2}}+∥D_K^{\beta }{\eta }_N^n∥+\lambda ∥{\eta }_N^n∥+∥R^n∥\right)}^2,\hfill \end{array}$

Taking ${\varphi }_N^{*n}={\omega }^{\mu /2}{\mathcal{P}}_N^{\mu /2}{\tilde{\varphi }}^n$ and by Lemma 7, we obtain the following:

(43)$\begin{array}{cc}\hfill D_K^{\beta }{∥{\tilde{e}}_N^n∥}_{{\omega }^{\mu /2}}^2& \le C{\left({∥(I-{\mathcal{P}}_N^{\mu /2}){(-\Delta )}^{\mu /2}\varphi ∥}_{{\omega }^{\mu /2}}+∥D_K^{\beta }{\eta }_N^n∥+\lambda ∥{\eta }_N^n∥+∥R^n∥\right)}^2\hfill \\ & \le C{\left({∥(I-{\mathcal{P}}_N^{\mu /2}){(-\Delta )}^{\mu /2}\varphi ∥}_{{\omega }^{\mu /2}}+∥D_K^{\beta }{\eta }_N^n∥+\lambda ∥{\eta }_N^n∥+∥R^n∥\right)}^2\hfill \\ \hfill & \le C{\left(C{N}^{\mu -\gamma }{\left|{\tilde{\varphi }}^n\right|}_{B_{{\omega }^{\mu /2}}^{\gamma }}+C{N}^{\mu -\gamma }{\left|D_K^{\beta }{\tilde{\varphi }}^n\right|}_{B_{{\omega }^{\mu /2}}^{\gamma }}+∥R^n∥\right)}^2.\hfill \end{array}$

Using the non-uniform discrete fractional Grönwall inequality and the maximum principle, we derive the following:

(44)${∥e_N^n∥}_{{\omega }^{-\mu /2}}\le C\left(N^{\mu -\gamma }+R^n\right),\forall 1\le n\le K.$

Further by the triangle inequality, we obtain the following equation:

${∥\varphi -{\varphi }_N^n∥}_{{\omega }^{-\mu /2}}={∥e_N^n-{\eta }_N^n∥}_{{\omega }^{-\mu /2}}\le \mathcal{C}\left({∥e_N^n∥}_{{\omega }^{-\mu /2}}+{∥\eta ∥}_{{\omega }^{-\mu /2}}\right)\le C\left(N^{\mu -\gamma }+K^{-min\{\rho \beta ,2-\beta \}}\right),$

where C is a positive constant independent of K and N. □

6. Numerical Results

We consider the following time-space fractional diffusion equation with a non-smooth solution in both time and space:

(45)${}_0{}^{C}D_t^{\beta }\varphi (x,t)+{(-\Delta )}^{\mu /2}\varphi (x,t)+\lambda \varphi (x,t)=f(x,t),x\in I,t\in (0,1],$

where $f(x,t)={(1-x^2)}^{\mu /2}\left({\displaystyle \frac{t^{\sigma -\beta }}{\mathsf{\Gamma }(\sigma -\beta +1)}}+{\displaystyle \frac{t^{\sigma }}{\mathsf{\Gamma }(\sigma +1)}}\right)+\mathsf{\Gamma }(1+\mu ){\displaystyle \frac{t^{\sigma }}{\mathsf{\Gamma }(\sigma +1)}}$. The exact solution is $u(x,t)={(1-x^2)}^{\mu /2}\left({\displaystyle \frac{t^{\sigma }}{\mathsf{\Gamma }(\sigma +1)}}\right).$

We demonstrate the convergence order and temporal errors by selecting N to be large in order to eliminate the spatial error. For $N=20$, the errors for various $\rho$ and $\beta$ are illustrated in Table 1 and Table 2. The theoretical analysis is consistent with the finding that $min\left\{\rho \sigma ,2-\beta \right\}$-order temporal accuracy has been achieved. Additionally, we confirm the spatial accuracy by selecting a sufficiently large K to prevent contamination of the temporal error. Table 3 illustrates the errors with regard to the polynomial degree N for various $\mu$ when $K=1000$ and $\rho =3$.

7. Conclusions

We have established a finite difference/Petrov–Galerkin spectral approach to solve the space-time fractional reaction–diffusion equation, where the solution is inherently non-smooth in both the temporal and spatial directions. The stability and convergence of the fully discrete scheme are analyzed. The results show that the temporal convergence order is $min\{\rho \sigma ,2-\beta \}$. Future work will extend this framework to nonlinear fractional diffusion equations, where coupling between nonlinearities and fractional operators introduces additional analytical and computational complexity. This extension aims to broaden applicability to problems such as reaction–diffusion systems with nonlocal interactions or phase-field models with memory effects, further bridging the gap between theoretical fractional calculus and practical computational science.

Data Availability Statement

All data supporting the findings of this study are available within the article.

Conflicts of Interest

The author declares 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.

Tables

Table 1

Numerical temporal accuracy for $N=20$, $\sigma =\beta =0.4$, and $\mu =1.3$.

K	$\mathit{\rho }=1$		$\mathit{\rho }=2$		$\mathit{\rho }=4$
	Error	Order	Error	Order	Error	Order
50	$4.084\times {10}^{-2}$	-	$1.134\times {10}^{-2}$	-	$1.217\times {10}^{-3}$	-
100	$3.334\times {10}^{-2}$	0.287	$6.771\times {10}^{-3}$	0.745	$4.313\times {10}^{-4}$	1.496
200	$2.704\times {10}^{-2}$	0.308	$3.978\times {10}^{-3}$	0.767	$1.496\times {10}^{-4}$	1.528
400	$2.157\times {10}^{-2}$	0.326	$2.316\times {10}^{-3}$	0.781	$5.130\times {10}^{-5}$	1.544
800	$1.703\times {10}^{-2}$	0.341	$1.340\times {10}^{-3}$	0.789	$1.746\times {10}^{-5}$	1.555
1000	$1.575\times {10}^{-2}$	0.350	$1.123\times {10}^{-4}$	0.792	$1.232\times {10}^{-5}$	1.560
$min\left\{r\sigma ,2-\beta \right\}$		$0.400$		$0.800$		$1.600$

Table 2

Numerical temporal accuracy for $N=20,$$\sigma =\beta =0.8$, and $\mu =1.3$.

K	$\mathit{\rho }=1$		$\mathit{\rho }=\frac{3}{2}$		$\mathit{\rho }=2$
	Error	Order	Error	Order	Error	Order
50	$7.641\times {10}^{-3}$	-	$2.755\times {10}^{-3}$	-	$1.681\times {10}^{-3}$	-
100	$4.701\times {10}^{-3}$	0.700	$1.357\times {10}^{-3}$	1.021	$7.598\times {10}^{-4}$	1.146
200	$2.814\times {10}^{-3}$	0.740	$6.527\times {10}^{-4}$	1.057	$3.379\times {10}^{-4}$	1.169
400	$1.665\times {10}^{-3}$	0.756	$3.077\times {10}^{-4}$	1.085	$1.489\times {10}^{-4}$	1.182
800	$9.762\times {10}^{-4}$	0.770	$1.428\times {10}^{-4}$	1.107	$6.529\times {10}^{-5}$	1.190
1000	$8.2034\times {10}^{-4}$	0.780	$1.112\times {10}^{-4}$	1.119	$5.003\times {10}^{-5}$	1.193
$min\left\{r\sigma ,2-\beta \right\}$		$0.8000$		$1.200$		$1.200$

Table 3

Numerical spacial accuracy for $\rho =3$, and $M=1000$.

N	$\mathit{\mu }=1.5$	$\mathit{\mu }=1.9$
10	$2.153\times {10}^{-5}$	$2.079\times {10}^{-5}$
20	$2.153\times {10}^{-5}$	$2.079\times {10}^{-5}$
30	$2.153\times {10}^{-5}$	$2.079\times {10}^{-5}$
40	$2.153\times {10}^{-5}$	$2.079\times {10}^{-5}$
50	$2.153\times {10}^{-5}$	$2.079\times {10}^{-5}$

References

1. Murray, J.D. Mathematical Biology: I. An Introduction; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2007; Volume 17.

2. del Castillo-Negrete, D.; Carreras, B.; Lynch, V. Front Dynamics in Reaction-Diffusion Systems with Levy Flights: A Fractional Diffusion Approach. Phys. Rev. Lett.; 2003; 91, 018302. [DOI: https://dx.doi.org/10.1103/PhysRevLett.91.018302]

3. Hanert, E.; Schumacher, E.; Deleersnijder, E. Front dynamics in fractional-order epidemic models. J. Theor. Biol.; 2011; 279, pp. 9-16. [DOI: https://dx.doi.org/10.1016/j.jtbi.2011.03.012]

4. Valdinoci, E. From the long jump random walk to the fractional Laplacian. SeMA J. Boletín Soc. Espa Nola Matemática Apl.; 2009; 355, pp. 33-44.

5. Kullberg, A.; del Castillo-Negrete, D.; Morales, G.; Maggs, J. Isotropic model of fractional transport in two-dimensional bounded domains. Phys. Rev. E-Stat. Nonlinear Soft Matter Phys.; 2013; 87, 052115. [DOI: https://dx.doi.org/10.1103/PhysRevE.87.052115] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23767495]

6. Izsák, F.; Szekeres, B.J. Models of space-fractional diffusion: A critical review. Appl. Math. Lett.; 2017; 71, pp. 38-43. [DOI: https://dx.doi.org/10.1016/j.aml.2017.03.006]

7. Duo, S.; Zhang, Y. Accurate numerical methods for two and three dimensional integral fractional Laplacian with applications. Comput. Methods Appl. Mech. Eng.; 2019; 355, pp. 639-662. [DOI: https://dx.doi.org/10.1016/j.cma.2019.06.016]

8. Wei, L.; Yang, Y. Optimal order finite difference/local discontinuous Galerkin method for variable-order time-fractional diffusion equation. J. Comput. Appl. Math.; 2021; 383, 113129. [DOI: https://dx.doi.org/10.1016/j.cam.2020.113129]

9. Kheybari, S.; Darvishi, M.T.; Alizadeh, F.; Hincal, E.; Hosseini, K.; Chen, X. Solving Multi-Point Boundary Value Differential Equations with a Novel Semi-Analytical Algorithm. Int. J. Comput. Methods; 2025; 22, 2450086. [DOI: https://dx.doi.org/10.1142/S0219876224500865]

10. Zhang, X.; Feng, Y.; Luo, Z.; Liu, J. A spatial sixth-order numerical scheme for solving fractional partial differential equation. Appl. Math. Lett.; 2025; 159, 109265. [DOI: https://dx.doi.org/10.1016/j.aml.2024.109265]

11. Alsaedi, A.; Kirane, M.; Torebek, B.T. Global existence and blow-up for a space and time nonlocal reaction-diffusion equation. Quaest. Math.; 2021; 44, pp. 747-753. [DOI: https://dx.doi.org/10.2989/16073606.2020.1745923]

12. Kolokoltsov, V.N.; Veretennikova, M. Well-posedness and regularity of the Cauchy problem for nonlinear fractional in time and space equations. Fract. Differ. Calc.; 2014; 4, pp. 1-30. [DOI: https://dx.doi.org/10.7153/fdc-04-01]

13. Jia, J.; Li, K. Maximum principles for a time–space fractional diffusion equation. Appl. Math. Lett.; 2016; 62, pp. 23-28. [DOI: https://dx.doi.org/10.1016/j.aml.2016.06.010]

14. Gu, X.M.; Sun, H.W.; Zhang, Y.; Zhao, Y.L. Fast implicit difference schemes for time-space fractional diffusion equations with the integral fractional Laplacian. Math. Methods Appl. Sci.; 2021; 44, pp. 441-463. [DOI: https://dx.doi.org/10.1002/mma.6746]

15. Hao, Z.; Zhang, Z.; Du, R. Fractional centered difference scheme for high-dimensional integral fractional Laplacian. J. Comput. Phys.; 2021; 424, 109851. [DOI: https://dx.doi.org/10.1016/j.jcp.2020.109851]

16. Zhou, J.; Chen, H. Finite Element Method on locally refined composite meshes for Dirichlet fractional Laplacian. J. Comput. Sci.; 2024; 82, 102433. [DOI: https://dx.doi.org/10.1016/j.jocs.2024.102433]

17. Huang, Y.; Oberman, A. Numerical methods for the fractional Laplacian: A finite difference-quadrature approach. SIAM J. Numer. Anal.; 2014; 52, pp. 3056-3084. [DOI: https://dx.doi.org/10.1137/140954040]

18. Gao, T.; Duan, J.; Li, X.; Song, R. Mean exit time and escape probability for dynamical systems driven by Lévy noises. SIAM J. Sci. Comput.; 2014; 36, pp. A887-A906. [DOI: https://dx.doi.org/10.1137/120897262]

19. Acosta, G.; Borthagaray, J.P. A fractional Laplace equation: Regularity of solutions and finite element approximations. SIAM J. Numer. Anal.; 2017; 55, pp. 472-495. [DOI: https://dx.doi.org/10.1137/15M1033952]

20. Ainsworth, M.; Glusa, C. Aspects of an adaptive finite element method for the fractional Laplacian: A priori and a posteriori error estimates, efficient implementation and multigrid solver. Comput. Methods Appl. Mech. Eng.; 2017; 327, pp. 4-35. [DOI: https://dx.doi.org/10.1016/j.cma.2017.08.019]

21. Tian, X.; Du, Q.; Gunzburger, M. Asymptotically compatible schemes for the approximation of fractional Laplacian and related nonlocal diffusion problems on bounded domains. Adv. Comput. Math.; 2016; 42, pp. 1363-1380. [DOI: https://dx.doi.org/10.1007/s10444-016-9466-z]

22. Mao, Z.; Shen, J. Hermite spectral methods for fractional PDEs in unbounded domains. SIAM J. Sci. Comput.; 2017; 39, pp. A1928-A1950. [DOI: https://dx.doi.org/10.1137/16M1097109]

23. Tang, T.; Wang, L.L.; Yuan, H.; Zhou, T. Rational spectral methods for PDEs involving fractional Laplacian in unbounded domains. SIAM J. Sci. Comput.; 2020; 42, pp. A585-A611. [DOI: https://dx.doi.org/10.1137/19M1244299]

24. Abbaszadeh, M.; Zaky, M.A.; Dehghan, M. Virtual element approximation and BDF2 time-discrete scheme for a partial integro-differential equation with a singular Abel’s kernel. Appl. Math. Comput.; 2025; 501, 129451. [DOI: https://dx.doi.org/10.1016/j.amc.2025.129451]

25. Moussa, H.; Saker, M.; Zaky, M.; Babatin, M.; Ezz-Eldien, S. Mapped Legendre-spectral method for high-dimensional multi-term time-fractional diffusion-wave equation with non-smooth solution. Comput. Appl. Math.; 2025; 44, 167. [DOI: https://dx.doi.org/10.1007/s40314-025-03123-z]

26. Ghoreyshi, A.; Abbaszadeh, M.; Zaky, M.A.; Dehghan, M. Finite block method for nonlinear time-fractional partial integro-differential equations: Stability, convergence, and numerical analysis. Appl. Numer. Math.; 2025; 214, pp. 82-103. [DOI: https://dx.doi.org/10.1016/j.apnum.2025.03.002]

27. 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]

28. Zaky, M.A.; Al Kenany, A.; Alhazmi, S.; Ezz-Eldien, S.S. Numerical treatment of tempered space–fractional Zeldovich–Frank–Kamenetskii equation. Rom. Rep. Phys.; 2025; Available online: https://rrp.nipne.ro/IP/AP821.pdf (accessed on 10 May 2025).

29. Guo, B.y.; Wang, L.l. Jacobi approximations in non-uniformly Jacobi-weighted Sobolev spaces. J. Approx. Theory; 2004; 128, pp. 1-41. [DOI: https://dx.doi.org/10.1016/j.jat.2004.03.008]

30. Babuska, I.; Guo, B. Direct and inverse approximation theorems for the p-version of the finite element method in the framework of weighted Besov spaces. Part I: Approximability of functions in the weighted Besov spaces. SIAM J. Numer. Anal.; 2002; 39, pp. 1512-1538. [DOI: https://dx.doi.org/10.1137/S0036142901356551]

31. Fdez-Manin, G.; Munoz-Sola, R. Polynomial approximation of some singular solutions in weighted Sobolev space. Proceedings of the Third International Conference on Spectral and Higher Order Methods; Houston, TX, USA, 5–9 June 1995.

32. Stynes, M.; O’Riordan, E.; Gracia, J.L. Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation. SIAM J. Numer. Anal.; 2017; 55, pp. 1057-1079. [DOI: https://dx.doi.org/10.1137/16M1082329]

33. Chen, H.; Stynes, M. Blow-up of error estimates in time-fractional initial-boundary value problems. IMA J. Numer. Anal.; 2021; 41, pp. 974-997. [DOI: https://dx.doi.org/10.1093/imanum/draa015]

34. Acosta, G.; Borthagaray, J.P.; Bruno, O.; Maas, M. Regularity theory and high order numerical methods for the (1D)-fractional Laplacian. Math. Comput.; 2018; 87, pp. 1821-1857. [DOI: https://dx.doi.org/10.1090/mcom/3276]

35. Zeng, F.; Liu, F.; Li, C.; Burrage, K.; Turner, I.; Anh, V. A Crank–Nicolson ADI spectral method for a two-dimensional Riesz space fractional nonlinear reaction-diffusion equation. SIAM J. Numer. Anal.; 2014; 52, pp. 2599-2622. [DOI: https://dx.doi.org/10.1137/130934192]

36. Huang, C.; Stynes, M. α-robust error analysis of a mixed finite element method for a time-fractional biharmonic equation. Numer. Algorithms; 2021; 87, pp. 1749-1766. [DOI: https://dx.doi.org/10.1007/s11075-020-01036-y]

37. Chen, H.; Hu, X.; Ren, J.; Sun, T.; Tang, Y. L1 scheme on graded mesh for the linearized time fractional KdV equation with initial singularity. Int. J. Mod. Sim. Sci. Comp.; 2019; 10, 18. [DOI: https://dx.doi.org/10.1142/S179396231941006X]

38. Ma, H.; Sun, W. Optimal error estimates of the Legendre–Petrov–Galerkin method for the Korteweg–de Vries equation. SIAM J. Numer. Anal.; 2001; 39, pp. 1380-1394. [DOI: https://dx.doi.org/10.1137/S0036142900378327]

39. Zhang, Z. Error estimates of spectral Galerkin methods for a linear fractional reaction–diffusion equation. J. Sci. Comput.; 2019; 78, pp. 1087-1110. [DOI: https://dx.doi.org/10.1007/s10915-018-0800-0]