# A Robust Numerical Scheme for Solving Riesz-Tempered Fractional Reaction-
Diffusion Equations

Mohammad Partohaghighi111Department of Applied Mathematics, University of
California Merced, Merced, CA, USA; e-mail<EMAIL_ADDRESS>,
Emmanuel Asante-Asamani222Department of Mathematics, Clarkson University,
Potsdam, NY, USA<EMAIL_ADDRESS>, Olaniyi S. Iyiola333Department of
Mathematics, Morgan State University, Baltimore, MD, USA;
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The Fractional Diffusion Equation (FDE) is a mathematical model that describes
anomalous transport phenomena characterized by non-local and long-range
dependencies which deviate from the traditional behavior of diffusion. Solving
this equation numerically is challenging due to the need to discretize
complicated integral operators which increase the computational costs. These
complexities are exacerbated by nonlinear source terms, nonsmooth data and
irregular domains. In this study, we propose a second order Exponential Time
Differencing Finite Element Method (ETD-RDP-FEM) to efficiently solve
nonlinear FDE, posed in irregular domains. This approach discretizes matrix
exponentials using a rational function with real and distinct poles, resulting
in an L-stable scheme that damps spurious oscillations caused by non-smooth
initial data. The method is shown to outperform existing second-order methods
for FDEs with a higher accuracy and faster computational time.

Keywords: Exponential time differencing; Finite element method; tempered
fractional operator; Reaction-diffusion equation.

2010 MSC classification:

##  1 Introduction

Fractional calculus is a developing field in mathematics that explores non-
integer order derivatives and integrals. Recently, this area has attracted
significant attention due to its broad range of applications in various
scientific areas including chemistry, biology, engineering, and finance, among
others [1, 2, 3, 4, 5, 6, 7, 8]. One of the main drivers for the study of
fractional calculus is its ability to provide models of physical phenomena
that are difficult to explain using classical calculus. For example, it has
been used to model viscoelastic materials [9] which show time-dependent
deformation and stress relaxation as well as to describe particle diffusion in
complex media [10]. Caputo, Riemann-Liouville, Reisz, Grünwald-Letnikov,
Tempered fractional derivatives are among the several types of non-integer
order derivatives with specific associated applications and/or properties, see
[11, 12].

The fractional diffusion equation describes the behavior of particles
undergoing anomalous diffusion and/or chemical reactions. This phenomenon is
characterized by a non-linear correlation between mean squared displacement
and time. Anomalous diffusion is prevalent in various physical, biological,
and chemical systems including porous media, polymers, and biological tissues
[13, 14, 16, 15]. The tempered fractional diffusion equation (TFDE) adjusts
the standard second order spatial derivative into a tempered fractional
derivative. This adjustment refines the parameters of random walk models ruled
by an exponentially tempered power-law jump distribution. In finance, one
typically observes the creation of semi-heavy tails due to the restricted
tempered stable probability densities. Tempered fractional time derivatives
emerge from waiting times governed by a tempered power law, which have proven
significantly valuable in geophysics. When it comes to Brownian motion, the
associated tempered fractional derivative or integral, otherwise known as
tempered fractional Brownian motion, can exhibit semi-long range dependence
[17].

Numerical solution of TFDE comes with several notable challenges. The first is
the so-called curse of non-locality. The characteristic of non-locality in
fractional diffusion equations stems from the employment of fractional
derivatives. These are an extension of conventional derivatives but to non-
integer orders. Consequently, unlike the typical derivatives that depend only
on the function’s behavior in the immediate vicinity of a specific point,
fractional derivatives rely on the function’s characteristics across the
entire domain. This implies that the value of a fractional derivative at a
given point is influenced by the function’s behavior everywhere, not just near
that particular point. Another layer of complexity comes from the memory
effect of the tempered fractional operator, where the future state of the
system is dependent on all of its past states. The memory effect implies long-
range dependence, whereby states from the distant past can still significantly
influence the current state. This feature complicates the truncation method
commonly used in numerical solutions to limit the influence of previous
states, making the system far more intricate to solve. The third challenge
lies in nonlinearity in the reaction term, which typically requires costly
iterative methods resulting from an implicit discretization of temporal
derivatives. The fourth challenge, non smoothness in the initial data, results
in spurious oscillations which can blow up the numerical solution if not
sufficiently damped. Finally, TFDEs posed in irregular spatial domains are
usually difficult to discretize with finite difference schemes.

Several strategies have been implemented to solve the tempered fractional
diffusion equation. In [20] a Forward Euler (FE) method is combined with a
weighted and shifted Lubich difference (WSLD) scheme to solve a linear FTDE.
The use of an explicit FE method introduces a stability restriction on the
time step with can lead to costly simulations.  In [19] the Crank-Nicolson
(CN) method is combined with a Finite Element Method (FEM) to resolve the
challenge with solving TFDE in irregular domains and remove the stability
restriction on the time step. However, the study did not investigate the
nonlinear problem, which is likely to be costly since it will require an
iterative scheme to evolve the implicit CN method. The finite element method
as been applied to solve other fractional boundary value problems involving
Riemann-Liouville or Caputo derivative [28], fractional advection-diffusion
reaction equations[29], and fractional diffusion-wave equation [30].

The family of Exponential Time Differencing (ETD) schemes avoid costly
iterative procedures such as Newton’s method by utilizing an implicit
discretization of the stiff, typically linear, terms in a differential
equation (DE) and an explicit treatment of the nonstiff, typically nonlinear,
terms to achieve excellent stability properties at a reduced computational
cost[31, 32]. In [21, 22] a second order L-stable exponential time
differencing (ETD-RDP) method is introduced for non-linear equations of
reaction-diffusion type (RDE) and is shown to efficiently damp out spurious
oscillations resulting from nonsmooth initial data. The method was initially
applied to integer RDE and later extended to fractional RDE in [23, 25].

In this paper, we consider the following tempered fractional diffusion
equation with the order $\alpha$

$\begin{split}u_{t}(x,t)-\partial_{\mid
x\mid}^{\alpha,\lambda}u(x,t)&=f(x,t,u(x,t)),\\\ \end{split}$ (1)

where

$\partial_{\mid
x\mid}^{\alpha,\lambda}u(x,t)=C_{\alpha}\bigg{(}\partial_{-x}^{\alpha,\lambda}u(x,t)+\partial_{+x}^{\alpha,\lambda}u(x,t)\bigg{)},\
C_{\alpha}=\frac{1}{2cos(0.5\pi)},$

with the initial and boundary conditions

$\begin{split}u(0,t)&=u(1,t)=0,\\\ u(x,0)&=g(x),\end{split}$ (2)

where $\lambda$ is the tempering parameter, the source function here is
$f(x,t,u(x,t))$, and the derivatives $\partial_{-x}^{\alpha,\lambda}u(x,t)$
and $\partial_{+x}^{\alpha,\lambda}u(x,t)$ are defined as follows

$\displaystyle\partial_{+x}^{\alpha,\lambda}u(x,t)$
$\displaystyle=D_{+x}^{\alpha,\lambda}u(x,t)-\lambda^{\alpha}u(x,t)-\alpha\lambda^{\alpha-1}D_{x}u(x,t)$
$\displaystyle\partial_{-x}^{\alpha,\lambda}u(x,t)$
$\displaystyle=D_{-x}^{\alpha,\lambda}u(x,t)-\lambda^{\alpha}u(x,t)+\alpha\lambda^{\alpha-1}D_{x}u(x,t),$

also, $D_{+x}^{\alpha,\lambda}u(x,t)$ and $D_{-x}^{\alpha,\lambda}u(x,t)$
denote the left and right tempered fractional derivatives of the function
$u(x,t)$, respectively and $D_{x}$ is the general derivatie operator. We
propose to utilize a novel combination of Exponential Time Differencing (ETD)
together with Finite Element Method (FEM) to solve the tempered fractional
diffusion equation. We adopt a method of lines approach and begin with an FEM
discretization of the model in space to yield a system of ordinary
differential equations which are then solved in time with the ETDRDP scheme.

Our contribution. This paper introduces an innovative strategy, the
Exponential Time Differencing Finite Element Method (ETDRDP-FEM), which unites
the Finite Element Method (FEM) and Exponential Time Differencing (ETDRDP)
scheme to solve the nonlinear TFDE. This combined approach is designed to
solve various tempered  fractional models efficiently, encompassing linear,
nonlinear, and nonsmooth problems in regular or irregular domains. By
amalgamating the strengths of both FEM and ETDRDP, the method we propose
aspires to offer superior accuracy and computational efficacy in tackling
intricate fractional problems, whilst preserving its versatility across a
broad spectrum of situations.

Outline. The organization of this paper is as follows: Section 2 introduces
the essential concepts and primary outcomes of fractional calculus necessary
for deriving the FEM method. The numerical discretizations of the fractional
diffusion equation, which includes the variational formulation of the primary
equation and its corresponding discretizations in space and time, are detailed
in Section 3. Section 5 delves into the numerical experiments performed.
Finally, Section 6 wraps up the paper with concluding observations.

##  2 Preliminaries

This section is devoted to recall some definitions and basic results necessary
to derive the variational formulation of the TFDE.

######  Definition 1.

For any $\alpha>0$, gamma function, denoted as $\Gamma$, is defined as

$\Gamma(\alpha)=\int_{0}^{\infty}\eta^{\alpha-1}e^{-\eta}d\eta,$

with the property of $\Gamma(\alpha+1)=\alpha\Gamma(\alpha).$

######  Definition 2.

For $\operatorname{Re}(\alpha)>0$ and $\operatorname{Re}(\beta)>0$, Beta
function, denoted as $B$, is defined as [26]

$B(\alpha,\beta)=\int_{0}^{1}\eta^{\alpha-1}(1-\eta)^{\beta-1}d\eta.$

The beta function is closely related to the gamma function, and it can be
expressed in terms of gamma function using the formula:

$B(\alpha,\beta)=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}.$
(3)

######  Definition 3 (Riemann-Liouville Integral).

[27] Let $f\in L^{1}[a,b]$ be a real-valued locally integrable function. Then
for $-\infty\leq a<w<b\leq\infty$ and $\alpha>0$, the left-sided Riemann-
Liouville (R-L) fractional integrals of a function $f$ is defined as

$I_{+w}^{\alpha}f(w)=\frac{1}{\Gamma(\alpha)}\int_{-\infty}^{w}(w-\eta)^{\alpha-1}f(\eta)\mathrm{d}\eta{\color[rgb]{1,0,0}.},$

and right-sided case is defined as

$I_{-w}^{\alpha}f(w)=\frac{1}{\Gamma(\alpha)}\int_{w}^{\infty}(w-\eta)^{\alpha-1}f(\eta)\mathrm{d}\eta.$

The order of the derivative here is $\alpha$.

######  Definition 4 (Riemann-Liouville Derivative, [27]).

Let $f\in L^{1}[a,b]$ and $n=\lfloor\alpha\rfloor+1$, $\alpha>0$. The Riemann-
Liouville (RL) fractional derivative (left- and right-sided) of order $\alpha$
of the function $f$ are given as

$D_{+w}^{\alpha}f(w)=\frac{1}{\Gamma(n-\alpha)}\frac{\mathrm{d}^{n}}{\mathrm{d}w^{n}}\int_{-\infty}^{w}(w-\eta)^{n-\alpha-1}f(\eta)\mathrm{d}\eta,$

and

$D_{-w}^{\alpha}f(w)=\frac{(-1)^{n}}{\Gamma(n-\alpha)}\frac{\mathrm{d}^{n}}{\mathrm{d}w^{n}}\int_{w}^{\infty}(\eta-w)^{n-\alpha-1}f(\eta)\mathrm{d}\eta,$

respectively.

######  Definition 5 (Tempered fractional integrals, [19]).

For $1\leq p<\infty$, and for any $f\in L^{p}(\mathbb{R})$, the left tempered
integral is defined as

$\begin{split}I_{+w}^{\alpha,\lambda}f(w)=\frac{1}{\Gamma(\alpha)}\int_{-\infty}^{w}(w-\eta)^{\alpha-1}e^{-\lambda(w-\eta)}f(\eta)\mathrm{d}\eta,\end{split}$
(4)

and the right tempered integral is defined as

$I_{-w}^{\alpha,\lambda}f(w)=\frac{1}{\Gamma(\alpha)}\int_{w}^{\infty}(\eta-w)^{\alpha-1}e^{-\lambda(\eta-w)}f(\eta)\mathrm{d}\eta,$
(5)

where $\alpha>0$ , $\lambda>0$ (the tempering parameter).

######  Definition 6 (Tempered Fractional derivative, [19]).

For any $\alpha>0$ and $\lambda>0$, the left fractional tempered derivatives
of a function $f:\mathbb{R}\to\mathbb{R}$ with order $\alpha$ is defined as

$\begin{split}D_{+w}^{\alpha,\lambda}f(w):&=e^{-\lambda
w}D_{+w}^{\alpha}e^{\lambda w}f(w)\\\ &=\frac{e^{-\lambda
w}}{\Gamma(n-\alpha)}\frac{\mathrm{d}^{n}}{\mathrm{d}w^{n}}\int_{-\infty}^{w}(w-\eta)^{n-\alpha-1}e^{\lambda\eta}f(\eta)\mathrm{d}\eta,\end{split}$
(6)

and the right is given as

$\begin{split}D_{-w}^{\alpha,\lambda}f(w):&=e^{\lambda
w}D_{-w}^{\alpha}e^{-\lambda w}f(w)\\\ &=\frac{(-1)^{n}e^{\lambda
x}}{\Gamma(n-\alpha)}\frac{\mathrm{d}^{n}}{\mathrm{d}w^{n}}\int_{w}^{\infty}(\eta-w)^{n-\alpha-1}e^{-\lambda\eta}f(\eta)\mathrm{d}\eta.\end{split}$
(7)

######  Lemma 1.

[18] Let $\lambda,\alpha>0$ and $f,g\in L^{2}(\mathbb{R})$. The following
identity is true for the left and right-tempered fractional integrals

$\bigg{\langle}I_{+}^{\alpha,\lambda}f,g\bigg{\rangle}=\bigg{\langle}f,I_{-}^{\alpha,\lambda}g\bigg{\rangle},$
(8)

which indicates that the two integrals are adjoint to each other.

######  Definition 7 (Tempered fractional derivative spaces, [19]).

Let $\alpha,\lambda>0$. $J_{+}^{\alpha,\lambda}$ and $J_{-}^{\alpha,\lambda}$
denote the left and right fractional tempered derivative spaces respectively
which are defined as

$J_{\mp}^{\alpha,\lambda}(\mathbb{R})=\\{z\in L^{2}(\mathbb{R});D{{}_{\mp
w}^{\alpha,\lambda}}z\in L^{2}(\mathbb{R})\\},$ (9)

with the norm

$\lVert z\lVert_{J_{\mp}^{\alpha,\lambda}(\mathbb{R})}=\lVert D{{}_{\mp
w}^{\alpha,\lambda}}z\rVert_{L^{2}(\mathbb{R})}.$ (10)

where $J_{\mp}^{\alpha,\lambda}(\mathbb{R})$ denotes the closure of
$\mathbb{C}_{0}^{\infty}(\mathbb{R})$ with respect to $\lVert
z\lVert_{J_{\mp}^{\alpha,\lambda}(\mathbb{R})}$.

######  Definition 8.

[19] The space $H^{\alpha,\lambda}$ for $\alpha,\lambda>0$ is define as
follows

$H^{\alpha,\lambda}(\mathbb{R})=\\{z\in
L^{2}(\mathbb{R});(\lambda^{2}+\omega^{2})^{\alpha/2}\hat{z}\in
L^{2}(\mathbb{R})\\}.$ (11)

and the associated norm is given by

$\lVert
z\lVert_{H^{\alpha,\lambda}(\mathbb{R})}=\lVert(\lambda^{2}+\omega^{2})^{\alpha/2}\hat{z}\lVert_{L^{2}(\mathbb{R})},$
(12)

where $\hat{z}$ represents the Fourier transformation $z$, and the Fourier
transformation variable $\omega$. In this case, the closure of
$\mathbb{C}_{0}^{\infty}(\mathbb{R})$ with respect to $\lVert
z\lVert_{H^{\alpha,\lambda}(\mathbb{R})}$ is denoted by
$H^{\alpha,\lambda}(\mathbb{R})$.

######  Definition 9.

[19] Let $\alpha>0$. Define the spaces $H^{\alpha}$ as follows

$H^{\alpha}(\mathbb{R})=\\{z\in
L^{2}(\mathbb{R});(1+\omega^{2})^{\alpha/2}\hat{z}\in L^{2}(\mathbb{R})\\}.$
(13)

with the norm

$\lVert
z\lVert_{H^{\alpha}(\mathbb{R})}=\lVert(1+\omega^{2})^{\alpha/2}\hat{z}\lVert_{L^{2}(\mathbb{R})}{\color[rgb]{1,0,0}.}$
(14)

######  Lemma 2.

[19] The following hold

$D_{\mp w}^{\alpha,\lambda}I_{\mp
w}^{\alpha,\lambda}f(w)=f(w){\color[rgb]{1,0,0},}$ (15)

for $\alpha,\,\lambda>0$.

######  Lemma 3.

[19] The following mappings

$I^{\alpha,\lambda}_{\mp x}:L_{2}(\Omega)\rightarrow
L_{2}(\Omega),\,\,I^{\alpha,\lambda}_{\mp x}:L_{2}(\Omega)\rightarrow
J^{\alpha,\lambda}_{\mp}(\Omega),\,\,\text{and}\,\,D^{\alpha,\lambda}_{\mp
x}:J^{\alpha,\lambda}_{\mp}(\Omega)\rightarrow L_{2}(\Omega)$

are bounded linear operators for $\lambda,\,\alpha>0$.

######  Lemma 4.

[19] The following hold

$I^{\alpha,\lambda}_{\mp w}D^{\alpha,\lambda}_{\mp w}z=z$

for every $z\in J^{\alpha,\lambda}_{\mp,0}(\Omega)$.

##  3 Spatial Discretization: Finite Element Method

The Galerkin finite element method (GFEM) is a widely utilized numerical
approach in solving partial differential equations (PDEs) in complicated
domains. The method begins with a variational formulation of the PDE followed
by an approximation of the weak solution in a finite dimensional function
space.

For completeness, we give the detailed proofs of some of the following results
about tempered fractional calculus. Readers should see [18] for further
details.

######  Theorem 5.

Let $\lambda,\alpha,\beta>0$ and $f\in L^{p}(\mathbb{R})$. The left-tempered
fractional integrations have the following property

$I_{+x}^{\alpha,\lambda}I_{+x}^{\beta,\lambda}f(x)=I_{+x}^{\alpha+\beta,\lambda}f(x).$
(16)

###### Proof.

Applying the left fractional operator $I_{+x}^{\alpha,\lambda}$ on
$I_{+x}^{\beta,\lambda}f(x)$, we have

$\displaystyle\begin{aligned}
&I_{+x}^{\alpha,\lambda}I_{+x}^{\beta,\lambda}f(x)=\frac{1}{\Gamma(\alpha)}\int_{-\infty}^{x}(x-t)^{\alpha-1}e^{-\lambda(x-t)}I_{+x}^{\beta,\lambda}f(t)\mathrm{d}t\\\
&=\frac{1}{\Gamma(\alpha)\Gamma(\beta)}\int_{-\infty}^{x}\int_{-\infty}^{t}(x-t)^{\alpha-1}e^{-\lambda(x-t)}(t-s)^{\beta-1}e^{-\lambda(t-s)}f(s)\mathrm{d}s\mathrm{d}t\\\
&=\frac{1}{\Gamma(\alpha)\Gamma(\beta)}\int_{-\infty}^{t}f(s)\int_{s}^{x}(x-t)^{\alpha-1}e^{-\lambda(x-t)}(t-s)^{\beta-1}e^{-\lambda(t-s)}\mathrm{d}t\mathrm{d}s,\end{aligned}$
(17)

Next, by changing variables $t=s+(x-s)r$, we have the following

$\begin{split}I_{+x}^{\alpha,\lambda}I_{+x}^{\beta,\lambda}f(x)&=\frac{1}{\Gamma(\alpha)\Gamma(\beta)}\int_{-\infty}^{t}(x-s)^{\alpha-1}(x-s)^{\beta-1}(x-s)f(s)\\\
&\times\int_{0}^{1}e^{-\lambda(x-s)(1-r)}e^{-\lambda(x-s)r}r^{\beta-1}(1-r)^{\alpha-1}\mathrm{d}r\mathrm{d}s\\\
&=\frac{1}{\Gamma(\alpha)\Gamma(\beta)}\int_{-\infty}^{t}(x-s)^{\alpha+\beta-1}e^{-\lambda(x-s)}f(s)\int_{0}^{1}r^{\beta-1}(1-r)^{\alpha-1}\mathrm{d}r\mathrm{d}s.\end{split}$
(18)

The result follows by using Equation 3.
Similar procedure can be used for the right-tempered fractional integral. ∎

######  Theorem 6.

[17] Let $\lambda,\alpha>0$ and $f\in L^{2}(\mathbb{R})$. The Fourier
transformations of the left and right-tempered fractional integrals are

$\mathcal{F}\bigg{[}I_{\pm x}^{\alpha,\lambda}f(x)\bigg{]}(\omega)=(\lambda\pm
i\omega)^{-\alpha}\hat{f}(x).$ (19)

###### Proof.

Recall that, for any $f\in L^{2}(\mathbb{R})$, its Fourier transform and
inverse Fourier transform are defined by

$\mathcal{F}\bigg{[}f(x)\bigg{]}(\omega)=\int_{-\infty}^{\infty}f(x)e^{-i\omega
x}\mathrm{d}x,$ (20)

and

$f(x)=\mathcal{F}^{-1}\bigg{[}\mathcal{F}[f](\omega)\bigg{]}(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathcal{F}[f](\omega)e^{-i\omega
x}\mathrm{d}\omega,$ (21)

Using the Heaviside function, $H(x)$, the the left-tempered fractional
integral can be reformulated as

$\begin{split}I_{+x}^{\alpha,\lambda}f(x)&=\frac{1}{\Gamma(\alpha)}\int_{-\infty}^{x}(x-\xi)^{\alpha-1}H(\xi)e^{-\alpha(x-\xi)}f(x)\mathrm{d}\xi\\\
&=(K*f)(x),\end{split}$ (22)

where

$K(x):=x^{\alpha-1}e^{-\lambda x}\frac{H(x)}{\Gamma(\alpha)}.$

Claim:

$\mathcal{F}[K](\omega)=(\lambda+i\omega)^{-\alpha}.$ (23)

###### Proof.

Using the following identity

$\mathcal{L}{(x^{\alpha-1};s)}=\frac{\Gamma(\alpha+1)}{s^{\alpha+1}},\ \
\alpha>0,$

we have

$\displaystyle\mathcal{F}[K](\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\alpha)}\int_{0}^{\infty}x^{\alpha-1}e^{-\lambda
x}e^{-i\omega x}\mathrm{d}x$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\alpha)}\int_{0}^{\infty}x^{\alpha-1}e^{-x(\lambda+i\omega)}\mathrm{d}x$
$\displaystyle=$ $\displaystyle(\lambda+i\omega)^{-\alpha}.$

∎

Then, using 23, we get

$\displaystyle\mathcal{F}\bigg{[}I_{\pm
x}^{\alpha,\lambda}f(x)\bigg{]}(\omega)$ $\displaystyle=$
$\displaystyle\mathcal{F}[K*f](\omega)$ $\displaystyle=$
$\displaystyle\mathcal{F}[K](\omega)\mathcal{F}[f](\omega)$ $\displaystyle=$
$\displaystyle(\lambda\pm i\omega)^{-\alpha}\mathcal{F}[f](\omega).$

∎

Hence, we obtain

$\mathcal{F}\bigg{[}I_{\pm x}^{\alpha,\lambda}f(x)\bigg{]}(\omega)=(\lambda\pm
i\omega)^{-\alpha}\hat{f}(x).$

Similarly, we have

$\mathcal{F}\bigg{[}D_{\pm}^{\alpha,\lambda}u(x)\bigg{]}(\omega)=\widehat{D_{\pm}^{{\alpha},\lambda}u(x)}=(\lambda\pm
i\omega)^{\alpha}\hat{u}(\omega).$ (24)

######  Theorem 7.

[19] The spaces $J_{\pm}^{\alpha,\lambda}(\mathbb{R})$ and
$H^{\alpha,\lambda}(\mathbb{R})$ are equal with equivalent norms.

###### Proof.

Using Equation 24, we have

$\begin{split}\lVert u\lVert_{J_{\pm}^{\alpha,\lambda}(\mathbb{R})}&=\lVert
D{{}_{\pm x}^{\alpha,\lambda}}u\lVert_{L_{2}(\mathbb{R})}\\\
&=\bigg{(}\int_{\mathbb{R}}D{{}_{+x}^{\alpha,\lambda}}u\overline{D{{}_{\mp
x}^{\alpha,\lambda}}u}dx\bigg{)}^{1/2}\\\
&=\bigg{(}\int_{\mathbb{R}}\widehat{D{{}_{\pm
x}^{\alpha,\lambda}}u}\overline{\widehat{D{{}_{\mp
x}^{\alpha,\lambda}}u}}d\omega\bigg{)}^{1/2}\\\
&=\left(\int_{\mathbb{R}}(\lambda\pm
i\omega)^{\alpha}\widehat{u(\omega)}(\lambda\mp
i\omega)^{\alpha}\widehat{u(\omega)}d\omega\right)^{1/2}\\\
&=\left(\int_{\mathbb{R}}(\lambda^{2}+\omega^{2})^{\alpha}\lvert\widehat{u(\omega)}\rvert^{2}d\omega\right)^{1/2}\\\
&=\lVert(\lambda^{2}+\omega^{2})^{\alpha/2}\widehat{u(\omega)}\lVert_{L^{2}(\mathbb{R})}\\\
&=\lVert u\lVert_{H^{\alpha,\lambda}(\mathbb{R})}.\end{split}$

∎

###  3.1 Variational Formulation Involving Fractional Operator

We seek a solution  $u\in L^{2}(0,T:H_{0}^{\frac{\alpha}{2}}(\Omega))$ such
that

$\langle u_{t},v\rangle-\langle\partial_{\mid
x\mid}^{\alpha,\lambda}u,v\rangle=\langle f,v\rangle,$ (25)

for all $v\in L^{2}(0,T:H_{0}^{\frac{\alpha}{2}}(\Omega))$.

Applying the definition of the tempered derivative to (25) we obtain

$\langle
u_{t},v\rangle-\mathcal{C}_{\alpha}\bigg{\langle}D_{+x}^{\alpha,\lambda}+D_{-x}^{\alpha,\lambda}u-2\lambda^{\alpha}u,v\bigg{\rangle}=\langle
f,v\rangle.$ (26)

This can be further simplified to

$\langle u_{t},v\rangle-\mathcal{C}_{\alpha}\langle
D_{+x}^{\alpha,\lambda}u,v\rangle-\mathcal{C}_{\alpha}\langle
D_{-x}^{\alpha,\lambda}u,v\rangle+2\mathcal{C}_{\alpha}\lambda^{\alpha}\langle
u,v\rangle=\langle f,v\rangle.$ (27)

By the linearity of the inner product. It remains to apply integration by
parts to simplify the inner products involving the tempered derivative. While
a discussion of this can be found in [19], we have added critical details
which we present here for completeness. Using that fact that
$D_{+x}^{\alpha,\lambda}u:=e^{-\lambda x}D_{+x}^{\alpha}e^{\lambda x}u$, we
have

$\displaystyle\langle
D_{+x}^{\alpha,\lambda}u,v\rangle=\int_{\Omega}(D_{+x}^{\alpha,\lambda}u)v\mathrm{d}x$
$\displaystyle=\int_{\Omega}(D_{+x}^{\alpha}e^{\lambda x}u)e^{-\lambda
x}v\mathrm{d}x.$

Using integrating by parts and the fact that

$D_{+x}^{\alpha}e^{\lambda x}u=D^{2}I_{+x}^{2-\alpha}(e^{\lambda
x}u)\quad\text{and}\quad v\in L^{2}(0,T:H_{0}^{\frac{\alpha}{2}}(\Omega)),$

we obtain

$\displaystyle\langle D_{+x}^{\alpha,\lambda}u,v\rangle$
$\displaystyle=\bigg{[}(e^{-\lambda x}v)I(D_{+x}^{\alpha}e^{\lambda
x}u)\bigg{]}_{\partial\Omega}-\int_{\Omega}ID_{+x}^{\alpha}(e^{\lambda
x}u)D(e^{-\lambda x}v)\mathrm{d}x$
$\displaystyle=-\bigg{\langle}ID_{+x}^{\alpha}(e^{\lambda x}u),D(e^{-\lambda
x}v)\bigg{\rangle}$
$\displaystyle=-\bigg{\langle}ID^{2}I_{+x}^{2-\alpha}(e^{\lambda
x}u),D(e^{-\lambda x}v)\bigg{\rangle}$
$\displaystyle=-\bigg{\langle}DI_{+x}^{2-\alpha}(e^{\lambda x}u),D(e^{-\lambda
x}v)\bigg{\rangle}$
$\displaystyle=-\bigg{\langle}I_{+x}^{2-\alpha}D(e^{\lambda x}u),D(e^{-\lambda
x}v)\bigg{\rangle}.$ (28)

Using Theorem 5 and Lemma 1 and the fact that
$I_{+x}^{\frac{2-\alpha}{2}}D=D_{+x}^{\frac{\alpha}{2}}$, we have

$\begin{split}\bigg{\langle}D_{+x}^{\alpha,\lambda}u,v\bigg{\rangle}&=-\bigg{\langle}I_{+x}^{2-\alpha}D(e^{\lambda
x}u),D(e^{-\lambda x}v)\bigg{\rangle}\\\
&=-\bigg{\langle}I_{+x}^{\frac{2-\alpha}{2}}I_{+x}^{\frac{2-\alpha}{2}}D(e^{\lambda
x}u),D(e^{-\lambda x}v)\bigg{\rangle}\\\
&=-\bigg{\langle}I_{+x}^{\frac{2-\alpha}{2}}D(e^{\lambda
x}u),I_{-x}^{\frac{2-\alpha}{2}}D(e^{-\lambda x}v)\bigg{\rangle}\\\
&=-\bigg{\langle}D_{+x}^{\frac{\alpha}{2}}(e^{\lambda
x}u),D_{-x}^{\frac{\alpha}{2}}(e^{-\lambda x}v)\bigg{\rangle}\\\
&=-\bigg{\langle}\mathcal{D}_{+x}^{\frac{\alpha}{2},\lambda}u,\mathcal{D}_{-x}^{\frac{\alpha}{2},\lambda}v\bigg{\rangle}.\end{split}$
(29)

Similarly,

$\displaystyle\begin{aligned}
\bigg{\langle}D_{-x}^{\alpha,\lambda}u,v\bigg{\rangle}=-\bigg{\langle}D_{-x}^{\frac{\alpha}{2},\lambda}u,D_{+x}^{\frac{\alpha}{2},\lambda}v\bigg{\rangle}.\end{aligned}$
(30)

Hence, we have from (27) for any $u,v\in H_{0}^{\frac{\alpha}{2}}(\Omega)$,
that

$\langle
u_{t},v\rangle-\mathcal{C}_{\alpha}\bigg{\langle}D_{+x}^{\frac{\alpha}{2},\lambda}u,D_{-x}^{\frac{\alpha}{2},\lambda}v\bigg{\rangle}-\mathcal{C}_{\alpha}\bigg{\langle}\mathcal{D}_{-x}^{\frac{\alpha}{2},\lambda}u,\mathcal{D}_{+x}^{\frac{\alpha}{2},\lambda}v\bigg{\rangle}+2\mathcal{C}_{\alpha}\lambda^{\alpha}\langle
u,v\rangle=\langle f,v\rangle.$ (31)

Then for each $t\in[0,T]$, we can write the variational form of the problem
presented in (1) as

$\langle u_{t},v\rangle+R\langle u,v\rangle=\langle f,\nu\rangle,\qquad\nu\in
H_{0}^{\frac{\alpha}{2}}(\Omega),$ (32)

where $R(u,v)$ (the bilinear form) is defined as

$R(u,v)=-\mathcal{C}_{\alpha}\bigg{\langle}D_{+x}^{\frac{\alpha}{2},\lambda}u,D_{-x}^{\frac{\alpha}{2},\lambda}v\bigg{\rangle}-\mathcal{C}_{\alpha}\bigg{\langle}D_{-x}^{\frac{\alpha}{2},\lambda}u,D_{+x}^{\frac{\alpha}{2},\lambda}v\bigg{\rangle}+2\mathcal{C}_{\alpha}\lambda^{\alpha}\langle
u,v\rangle,$ (33)

and the initial condition (2) is satisfied by $u\in
H_{0}^{\frac{\alpha}{2}}(\Omega)$.

Using Equation 24 for any $\alpha>0$ and $\lambda>0$, we have

$\begin{split}\bigg{\langle}D_{+x}^{\frac{\alpha}{2},\lambda}u(\omega),D_{-x}^{\frac{\alpha}{2},\lambda}v(\omega)\bigg{\rangle}&=\int_{\omega}D_{+x}^{\frac{\alpha}{2},\lambda}u(\omega)D_{-x}^{\frac{\alpha}{2},\lambda}v(\omega)\mathrm{d}x\\\
&=\int_{-\infty}^{\infty}\widehat{D_{+x}^{\frac{\alpha}{2},\lambda}u(\omega)}\widehat{D_{-x}^{\frac{\alpha}{2},\lambda}v(\omega)}\mathrm{d}x\\\
&=\int_{-\infty}^{\infty}\widehat{D_{+x}^{\frac{\alpha}{2},\lambda}u(\omega)}\overline{\widehat{D_{-x}^{\frac{\alpha}{2},\lambda}v(\omega)}}\mathrm{d}x\\\
&=\int_{-\infty}^{\infty}(\lambda+i\omega)^{\frac{\alpha}{2}}\hat{u}(\omega)\overline{(\lambda-i\omega)^{\frac{\alpha}{2}}\hat{v}(\omega)}\mathrm{d}\omega.\end{split}$
(34)

Note that

$\overline{(\lambda-i\omega)^{\frac{\alpha}{2}}}=\left\\{\begin{array}[]{ll}\overline{(\lambda+i\omega)^{\frac{\alpha}{2}}}e^{i\theta\alpha},&\omega\geq
0,\\\
\overline{(\lambda+i\omega)^{\frac{\alpha}{2}}}e^{-i\theta\alpha},&\omega<0,\end{array}\right.$
(35)

and

$\overline{(\lambda+i\omega)^{\frac{\alpha}{2}}}=\left\\{\begin{array}[]{ll}\overline{(\lambda-i\omega)^{\frac{\alpha}{2}}}e^{-i\theta\alpha},&\omega\geq
0,\\\
\overline{(\lambda-i\omega)^{\frac{\alpha}{2}}}e^{i\theta\alpha},&\omega<0.\end{array}\right.$
(36)

Then

$\begin{split}\bigg{\langle}D_{+x}^{\frac{\alpha}{2},\lambda}u(\omega),D_{-x}^{\frac{\alpha}{2},\lambda}v(\omega)\bigg{\rangle}&=\int_{-\infty}^{0}(\lambda+i\omega)^{\frac{\alpha}{2}}\hat{u}(\omega)\
\overline{(\lambda+i\omega)^{\frac{\alpha}{2}}\hat{v}(\omega)}e^{-i\theta\alpha}\mathrm{d}\omega\\\
&+\int_{0}^{\infty}(\lambda+i\omega)^{\frac{\alpha}{2}}\hat{u}(\omega)\
\overline{(\lambda+i\omega)^{\frac{\alpha}{2}}\hat{v}(\omega)}e^{i\theta\alpha}\mathrm{d}\omega,\\\
&=\int_{-\infty}^{0}(\lambda+i\omega)^{\frac{\alpha}{2}}\hat{u}(\omega)\
\overline{(\lambda+i\omega)^{\frac{\alpha}{2}}}\ \
\overline{\hat{v}(\omega)}e^{-i\theta\alpha}\mathrm{d}\omega\\\
&+\int_{0}^{\infty}(\lambda+i\omega)^{\frac{\alpha}{2}}\hat{u}(\omega)\
\overline{(\lambda+i\omega)^{\frac{\alpha}{2}}}\overline{\hat{v}(\omega)}e^{i\theta\alpha}\mathrm{d}\omega.\\\
\end{split}$ (37)

Equivalently, we have

$\displaystyle\begin{aligned}
\bigg{\langle}D_{+x}^{\frac{\alpha}{2},\lambda}u(\omega),D_{-x}^{\frac{\alpha}{2},\lambda}v(\omega)\bigg{\rangle}&=\int_{-\infty}^{0}(\lambda+i\omega)^{\frac{\alpha}{2}}\hat{u}(\omega)(\lambda-i\omega)^{\frac{\alpha}{2}}\
\ \overline{\hat{v}(\omega)}e^{-i\theta\alpha}\mathrm{d}\omega\\\
&+\int_{0}^{\infty}(\lambda+i\omega)^{\frac{\alpha}{2}}\hat{u}(\omega)(\lambda-i\omega)^{\frac{\alpha}{2}}\overline{\hat{v}(\omega)}e^{i\theta\alpha}\mathrm{d}\omega.\\\
\end{aligned}$ (38)

Hence,

$\begin{split}\bigg{\langle}D_{+x}^{\frac{\alpha}{2},\lambda}u(\omega),D_{-x}^{\frac{\alpha}{2},\lambda}v(\omega)\bigg{\rangle}&=\int_{-\infty}^{0}(\lambda^{2}+\omega^{2})^{\frac{\alpha}{2}}\hat{u}(\omega)\overline{\hat{v}(\omega)}e^{-i\theta\alpha}\mathrm{d}\omega\\\
&+\int_{0}^{\infty}(\lambda^{2}+\omega^{2})^{\frac{\alpha}{2}}\hat{u}(\omega)\overline{\hat{u}(\omega)}e^{i\theta\alpha}\mathrm{d}\omega.\\\
\end{split}$ (39)

Similarly, we have

$\begin{split}\bigg{\langle}D_{-x}^{\frac{\alpha}{2},\lambda}u(\omega),D_{+x}^{\frac{\alpha}{2},\lambda}v(\omega)\bigg{\rangle}&=\int_{-\infty}^{0}(\lambda^{2}+\omega^{2})^{\frac{\alpha}{2}}\hat{u}(\omega)\overline{\hat{v}(\omega)}e^{i\theta\alpha}\mathrm{d}\omega\\\
&+\int_{0}^{\infty}(\lambda^{2}+\omega^{2})^{\frac{\alpha}{2}}\hat{u}(\omega)\overline{\hat{u}(\omega)}e^{-i\theta\alpha}\mathrm{d}\omega.\\\
\end{split}$ (40)

From (39) and (40), we get

$\begin{split}\bigg{\langle}D_{+x}^{\frac{\alpha}{2},\lambda}u(\omega),D_{-x}^{\frac{\alpha}{2},\lambda}v(\omega)\bigg{\rangle}&+\bigg{\langle}D_{-x}^{\frac{\alpha}{2},\lambda}u(\omega),D_{+x}^{\frac{\alpha}{2},\lambda}v(\omega)\bigg{\rangle}\\\
&=\int_{-\infty}^{\infty}\bigg{(}\lambda^{2}+\omega^{2}\bigg{)}^{\frac{\alpha}{2}}\hat{u}(\omega)\overline{\hat{v}(\omega)}\bigg{(}e^{i\theta\alpha}+e^{-i\theta\alpha}\bigg{)}\mathrm{d}\omega.\end{split}$
(41)

Using

$(e^{i\theta\alpha}+e^{-i\theta\alpha})=2cos(\theta\alpha),$

and Equation 42, the relation (33) can be written as follows

$R(u,v)=-2\mathcal{C}_{\alpha}\int_{-\infty}^{\infty}\bigg{(}\lambda^{2}+\omega^{2}\bigg{)}^{\frac{\alpha}{2}}\hat{u}(\omega)\
\overline{\hat{v}}(\omega)\cos(\theta\alpha)\mathrm{d}\omega+2\mathcal{C}_{\alpha}\lambda^{\alpha}\bigg{\langle}u(\omega),v(\omega)\bigg{\rangle}.$
(42)

###  3.2 Discretization Of Variational Formulation

Here, we discretize the variational form in (32) to obtain a system of
ordinary differential equations. Suppose we have the following
representations.

$u(x,t)=\sum_{i=1}^{N-1}u_{i}(t)\varphi_{i}(x),\ \
v(x,t)=\sum_{j=1}^{N-1}v_{j}(t)\varphi_{j}(x),\ \
f(x,t)=\sum_{i=1}^{N-1}f_{i}\varphi_{i}(x).$ (43)

Substituting (43) into (32), we obtain

$\begin{split}\bigg{\langle}\sum_{i=1}^{N-1}u^{\prime}_{i}(t)\varphi_{i}(x),\sum_{j=1}^{N-1}v_{j}(t)\varphi_{j}(x)\bigg{\rangle}&+R\bigg{(}\sum_{i=1}^{N-1}u_{i}(t)\varphi_{i}(x),\sum_{j=1}^{N-1}v_{j}(t)\varphi_{j}(x)\bigg{)}\\\
&=\bigg{\langle}\sum_{i=1}^{N-1}f_{i}\varphi_{i}(x),\sum_{j=1}^{N-1}v_{j}(t)\varphi_{j}(x)\bigg{\rangle}.\end{split}$
(44)

We can write (44) as follows

$\begin{split}\sum_{j=1}^{N-1}v_{j}(t)\bigg{\langle}\sum_{i=1}^{N-1}u^{\prime}_{i}(t)\varphi_{i}(x),\varphi_{j}(x)\bigg{\rangle}&+\sum_{j=1}^{N-1}v_{j}(t)R\bigg{(}\sum_{i=1}^{N-1}u_{i}(t)\varphi_{i}(x),\varphi_{j}(x)\bigg{)}\\\
&=\sum_{j=1}^{N-1}v_{j}(t)\bigg{\langle}\sum_{i=1}^{N-1}f_{i}\varphi_{i}(x),\varphi_{j}(x)\bigg{\rangle}.\end{split}$
(45)

So, we have

$\begin{split}\bigg{\langle}\sum_{i=1}^{N-1}u^{\prime}_{i}(t)\varphi_{i}(x),\varphi_{j}(x)\bigg{\rangle}&+R\bigg{(}\sum_{i=1}^{N-1}u_{i}(t)\varphi_{i}(x),\varphi_{j}(x)\bigg{)}\\\
&=\bigg{\langle}\sum_{i=1}^{N-1}f_{i}\varphi_{i}(x),\varphi_{j}(x)\bigg{\rangle},\
\ j=1,...,N-1.\end{split}$ (46)

Then

$\begin{split}\sum_{i=1}^{N-1}u^{\prime}_{i}(t)\bigg{\langle}\varphi_{i}(x),\varphi_{j}(x)\bigg{\rangle}&+\sum_{i=1}^{N-1}u_{i}(t)R\bigg{(}\varphi_{i}(x),\varphi_{j}(x)\bigg{)}\\\
&=\sum_{i=1}^{N-1}f_{i}\bigg{\langle}\varphi_{i}(x),\varphi_{j}(x)\bigg{\rangle}.\end{split}$
(47)

Consider $(N-1)\times(N-1)$ matrices $P$ and $G$ having entries

$p_{i,j}=\bigg{\langle}\varphi_{j}(x),\varphi_{i}(x)\bigg{\rangle},\quad\text{and}\quad
g_{i,j}=R\bigg{(}\varphi_{j}(x),\varphi_{i}(x)\bigg{)}.$

respectively.

Also,

$U=\bigg{[}u_{1}(t),\cdots,u_{N-1}(t)\bigg{]}^{T},\ \
U^{\prime}=\bigg{[}u^{\prime}_{1}(t),\cdots,u^{\prime}_{N-1}(t)\bigg{]}^{T},$

and $F=\bigg{[}f_{1}(t),\cdots,f_{N-1}(t)\bigg{]}^{T}$, we get, from Equation
(47)

$PU^{\prime}+GU=PF.$ (48)

Taking $B=P^{-1}G$, we obtained the system of ordinary differential equations
as follows

$\begin{split}U^{\prime}+BU=F.\end{split}$

##  4 Time Discretization: ETD-RDP and ETD-CN

In this section, our goal is to present overview of the time discretization
methods adopted to solve the system of ordinary differential equations derived
in the previous section

$U^{\prime}+BU=F(t,U).$ (49)

The focus is on exponential time differencing with real distinct poles
rational approximation to the exponential. For completeness, we also give the
Crank-Nicolson method. We should note here that both ETD-RDP and CN schemes
are second order accurate schemes.

###  4.1 ETD-RDP Scheme

We give brief steps in constructing the ETD-RDP scheme solving system of
ordinary differential equation in (49). Using a variation of the constants
formula on the time interval $[t_{m},t_{m+1}]$, we have

$U(t_{m+1})=e^{-\tau
B}U(t_{m})+\int_{t_{m}}^{t_{m+1}}e^{-B(t_{m+1}-s)}F(s,U(s))ds.$ (50)

Taking $s=t_{m}+z\tau$ with $t_{m}=m\tau$, $z\in[0,1]$ and $\tau$ is the time-
step, the following recurrence formula is obtained:

$U(t_{m+1})=e^{-\tau B}U(t_{m})+\tau\int_{0}^{1}e^{-\tau
B(1-z)}F(t_{m}+z\tau,U(t_{m}+z\tau))dz.$ (51)

When dealing with the discretizations of the integral part in (51), many
approximations result into different schemes with different properties such as
stability and convergence orders. In the case of ETD-RDP scheme, we seek a
second-order scheme in which a linear approximation of the nonlinear function
is used as follows:

$F(t_{m}+z\tau,U(t_{m}+z\tau))\approx
F(t_{m},U(t_{m}))+z\tau\frac{F(t_{m+1},U(t_{m+1}))-F(t_{m},U(t_{m}))}{\tau}.$
(52)

Using the above approximation in Equation (51), we obtain a semi-discretized
scheme with $U(t_{m})\approx U_{m}$

$\displaystyle U_{m+1}$ $\displaystyle=$ $\displaystyle
e^{-B\tau}U_{m}+B^{-1}\left[I-e^{-B\tau}\right]F(t_{m},U_{m})$
$\displaystyle+\frac{B^{-2}}{\tau}\left[\tau
B-I+e^{-B\tau}\right]\bigg{[}F(t_{m+1},U_{m+1})-F(t_{m},U_{m})\bigg{]}.$

By employing the constant approximation in (51) and integrating, first order
accuracy approximation is obtained as

$\displaystyle U_{m+1}^{*}=e^{-\tau B}U_{m}+B^{-1}\left[I-e^{-\tau
B}\right]F(t_{m},U_{m}),$ (53)

The final semi-discretized scheme is then obtained as:

$\displaystyle U_{m+1}^{*}$ $\displaystyle=$ $\displaystyle
e^{-B\tau}U_{n}+B^{-1}\left(I-e^{-B\tau}\right)F(t_{m},U_{m})$ (54)
$\displaystyle U_{m+1}$ $\displaystyle=$ $\displaystyle
e^{-B\tau}U_{m}+B^{-1}\left(I-e^{-B\tau}\right)F(t_{m},U_{m})$ (55)
$\displaystyle+\frac{B^{-2}}{\tau}\left[\tau
B-I+e^{-B\tau}\right]\bigg{[}F(t_{n+1},U^{*}_{m+1})-F(t_{m},U_{m})\bigg{]}.$

The primary computational incentive for approximating the exponential function
lies in its efficiency and practicality. This leads us to introduce the RDP
rational function as an approximation method [21, 25]. This approximation
technique is L-acceptable, ensuring a reliable and satisfactory outcome [21,
25]. Consider the following rational function approximation of the exponential
function.

$r(z)=\frac{1+\frac{5}{12}z}{(1-\frac{1}{4}z)(1-\frac{1}{3}z)}=\frac{9}{1-\frac{1}{3}z}-\frac{8}{1-\frac{1}{4}z}\approx
e^{z}.$

The function $r(z)$ is second order accurate. Assume $r(z)$ is defined and
bounded on the spectrum of $\tau B$, then

$\displaystyle e^{-\tau B}\approx r(-\tau B)$ $\displaystyle=$
$\displaystyle\left[I-\frac{5}{12}\tau B\right]\left[I+\frac{1}{4}\tau
B\right]^{-1}\left[I+\frac{1}{3}\tau B\right]^{-1}$ (56) $\displaystyle=$
$\displaystyle 9\left[I+\frac{1}{3}\tau B\right]^{-1}-8\left[I+\frac{1}{4}\tau
B\right]^{-1}.$

One major advantage of the above approximation is the fact that it can be
decomposed into partial fraction, allowing basic algebra to simplify the
scheme and also permit parallel implementation of the final scheme. Consider
the following partial fraction decomposition

$\displaystyle\left[I+\frac{B\tau}{4}\right]^{-1}\left[I+\frac{B\tau}{3}\right]^{-1}$
$\displaystyle\approx$ $\displaystyle 4\left[I+\frac{1}{3}\tau
B\right]^{-1}-3\left[I+\frac{1}{4}\tau B\right]^{-1}$ (57)
$\displaystyle\left[I+\frac{1}{6}\tau B\right]\left[I+\frac{1}{4}\tau
B\right]^{-1}\left[I+\frac{1}{3}\tau B\right]^{-1}$ $\displaystyle\approx$
$\displaystyle 2\left[I+\frac{1}{3}\tau B\right]^{-1}-\left[I+\frac{1}{4}\tau
B\right]^{-1}.$ (58)

Combining Equations (56), (57) and (58), we have the final ETD-RDP scheme as

$\displaystyle W_{m+1}^{*}$ $\displaystyle=$
$\displaystyle\left[I+B\tau\right]^{-1}\left[W_{n}+\tau
F(t_{m},W_{m})\right],$ (59) $\displaystyle W_{m+1}$ $\displaystyle=$
$\displaystyle\left[I+\frac{1}{3}B\tau\right]^{-1}\bigg{(}9W_{m}+2\tau
F(t_{m},W_{m})+\tau F(t_{m+1},W^{*}_{m+1})\bigg{]}$ (60)
$\displaystyle+\left[I+\frac{1}{4}B\tau\right]^{-1}\bigg{[}-8W_{m}-\frac{3\tau}{2}F(t_{m},W_{m})-\frac{\tau}{2}F(t_{m+1},W^{*}_{m+1})\bigg{]}.$

###  4.2 Crank-Nicolson Scheme

We present overview of Crank-Nicolson scheme here which is used for comparing
purposes in the implementation section. Consider the same system of ordinary
differential equations in (43) and the following approximations:

$\bar{\partial}U_{m}=\frac{U_{m}-U_{m-1}}{\tau},\ \
\tilde{U}_{m}=\frac{U_{m}+U_{m-1}}{2},\ \
F^{m-\frac{1}{2}}=\frac{1}{2}\bigg{[}F(U_{m},t_{m})+F(U_{m+1},t_{m+1})\bigg{]}.$

Hence, the Crank-Nicolson final scheme is given as

$\bigg{[}I+\frac{1}{2}B\tau\bigg{]}W_{m}=\bigg{[}I-\frac{1}{2}B\tau\bigg{]}W_{m-1}+F^{m-\frac{1}{2}},$
(61)

where Newton’s method is used to handle the nonlinear function.

##  5 Numerical Examples and Discussion

This section is focused on the implementation of the proposed scheme ETD-RDP
together with Finite Element Method (ETD-RDP-FEM) for solving fractional
reaction diffusion equations. The impact of the fractional order in the
problems considered on the efficacy of the scheme is further investigated. The
order of convergence is calculated using the formula [35]

$p\approx\frac{\log\left(\frac{E(\tau)}{E(\tau/2)}\right)}{\log(2)},\
\tau=\frac{T}{M},$

where $M$ denotes the number of discretization points in the temporal
direction and error $E(\tau)$ is given by:

$E(\tau)=\frac{\lVert u-\tilde{u}\rVert_{\infty}}{\lVert u\rVert_{\infty}},$

where $u$ represents the exact solution and $\tilde{u}$ represents the
approximate solution.

###  5.1 Numerical Examples

######  Example 5.1 (Non-homogeneous Linear Model).

Consider the following non-homogeneous linear model of Riesz-tempered
fraactional order derivative:

$\begin{split}u_{t}(x,t)-\partial_{\mid
x\mid}^{\alpha,\lambda}u(x,t)&=f(x,t),\ (x,t)\in[0,1]\times(0,1],\\\
u(0,t)&=u(1,t)=0,\\\ u(x,0)&=g(x),\end{split}$ (62)

with the exact solution

$u(x,t)=-2\lambda^{6-\alpha}\cos\left(\frac{\pi\alpha}{2}\right)(-\alpha)e^{-t}x^{3}(1-x)^{3},$

and function $f(x,t)$ is given by:

$\displaystyle f(x,t)$ $\displaystyle=$ $\displaystyle
2\lambda^{6-\alpha}\cos\left(\frac{\pi\alpha}{2}\right)\Gamma(-\alpha)e^{-t}x^{3}(1-x)^{3}$
$\displaystyle+e^{-t}\left(g_{1}(x)+\lambda\left[g_{2}(x)+g_{3}(x)+g_{4}(x)+g_{5}(x)+g_{6}(x)+g_{7}(x)\right]+g_{8}(x)+g_{9}(x)\right).$

where

$\displaystyle g_{1}(x)$ $\displaystyle=$ $\displaystyle
3\lambda^{5}(x-1)^{2}x^{2}(2x-1)\gamma(1-\alpha,x\lambda)$ $\displaystyle
g_{2}(x)$ $\displaystyle=$
$\displaystyle-\lambda^{5}x^{6}\gamma(-\alpha,x\lambda)+3\lambda^{5}x^{5}\gamma(-\alpha,x\lambda)-3\lambda^{5}x^{4}\gamma(-\alpha,x\lambda)-15\lambda^{3}x^{4}\gamma(2-\alpha,\lambda-x\lambda)$
$\displaystyle+\lambda^{5}x^{3}\gamma(-\alpha,x\lambda)-\lambda^{5}(x-1)^{3}x^{3}\gamma(-\alpha,\lambda-x\lambda)+30\lambda^{3}x^{3}\gamma(2-\alpha,\lambda-x\lambda)$
$\displaystyle g_{3}(x)$ $\displaystyle=$
$\displaystyle+20\lambda^{2}x^{3}\gamma(3-\alpha,x\lambda)-20\lambda^{2}x^{3}\gamma(3-\alpha,\lambda-x\lambda)-30\lambda^{2}x^{2}\gamma(3-\alpha,x\lambda)$
$\displaystyle-3\lambda^{4}(x-1)^{2}(2x-1)x^{2}\gamma(1-\alpha,\lambda-x\lambda)-18\lambda^{3}x^{2}\gamma(2-\alpha,\lambda-x\lambda)$
$\displaystyle+30\lambda^{2}x^{2}\gamma(3-\alpha,\lambda-x\lambda)$
$\displaystyle g_{5}(x)$ $\displaystyle=$ $\displaystyle-15\lambda
x^{2}\gamma(4-\alpha,x\lambda)-3\lambda^{3}(x-1)(5(x-1)x+1)x\gamma(2-\alpha,x\lambda)$
$\displaystyle-15\lambda x^{2}\gamma(4-\alpha,\lambda-x\lambda)$
$\displaystyle g_{6}(x)$ $\displaystyle=$
$\displaystyle+3\lambda^{3}x\gamma(2-\alpha,\lambda-x\lambda)+12\lambda^{2}x\gamma(3-\alpha,x\lambda)-12\lambda^{2}x\gamma(3-\alpha,\lambda-x\lambda)$
$\displaystyle-\lambda^{2}\gamma(3-\alpha,x\lambda)+\lambda^{2}\gamma(3-\alpha,\lambda-x\lambda)+15\lambda
x\gamma(4-\alpha,x\lambda)$ $\displaystyle+15\lambda
x\gamma(4-\alpha,\lambda-x\lambda)$ $\displaystyle g_{7}(x)$ $\displaystyle=$
$\displaystyle+6x\gamma(5-\alpha,x\lambda)-6x\gamma(5-\alpha,\lambda-x\lambda)-3\lambda\gamma(4-\alpha,x\lambda)-3\lambda\gamma(4-\alpha,\lambda-x\lambda)$
$\displaystyle-3\gamma(5-\alpha,x\lambda)+3\gamma(5-\alpha,\lambda-x\lambda)$
$\displaystyle g_{8}(x)$ $\displaystyle=$
$\displaystyle+2\Gamma(2-\alpha)\bigg{(}\alpha^{4}-14\alpha^{3}+71\alpha^{2}-154\alpha+3(\alpha-3)(\alpha-2)\lambda^{2}(5(x-1)x+1)$
$\displaystyle+3\lambda^{4}(x-1)x(5(x-1)x+1)+120\bigg{)}$ $\displaystyle
g_{9}(x)$ $\displaystyle=$
$\displaystyle-\gamma(6-\alpha,x\lambda)-\gamma(6-\alpha,\lambda-x\lambda),$

and the incomplete Gamma function, denoted as $\gamma(z,x)$, is defined as:

$\gamma(z,x)=\int_{0}^{x}t^{z-1}e^{-t}dt.$

As a proof of concept, we first introduce a test case example where the exact
solution is known. Although this is a linear model, the fractional order in
space makes it a difficult problem to solve due to nonlocal property effect.

Table 1: Example (5.1): Time rate of convergence and $L_{\infty}$-norm of ETD-RDP-FEM and CN-FEM with different values of $\alpha$ at fixed $t=1$. |  |  | Norm | Order
---|---|---|---|---
$\alpha$ | $h$ | $\tau$ | CN-FEM | ETD-RDP-FEM | CN-FEM | ETD-RDP-FEM
1.2 | 1/4 | 1/4 | $8.7000\times 10^{-3}$ | $8.4801\times 10^{-3}$ | 1.95623 | 1.91356
| 1/8 | 1/8 | $2.2000\times 10^{-3}$ | $2.3000\times 10^{-3}$ | 1.99015 | 2.07254
| 1/16 | 1/16 | $5.5377\times 10^{-4}$ | $5.4680\times 10^{-4}$ | 2.08542 | 2.05254
| 1/32 | 1/32 | $1.3048\times 10^{-4}$ | $1.3181\times 10^{-4}$ | 2.04652 | 2.03512
1.4 | 1/4 | 1/4 | $7.4000\times 10^{-4}$ | $6.1000\times 10^{-3}$ | 2.12561 | 2.08985
| 1/8 | 1/8 | $1.7000\times 10^{-3}$ | $1.4000\times 10^{-3}$ | 2.11021 | 2.06416
| 1/16 | 1/16 | $3.9374\times 10^{-4}$ | $8.1890\times 10^{-5}$ | 2.08141 | 2.01553
| 1/32 | 1/32 | $9.3034\times 10^{-5}$ | $2.0253\times 10^{-5}$ | 2.05750 | 2.01035
1.6 | 1/4 | 1/4 | $6.3352\times 10^{-4}$ | $1.2544\times 10^{-4}$ | 2.09511 | 2.08252
| 1/8 | 1/8 | $1.4828\times 10^{-4}$ | $2.9617\times 10^{-5}$ | 2.08856 | 2.05658
| 1/16 | 1/16 | $3.4863\times 10^{-5}$ | $7.1195\times 10^{-6}$ | 2.06758 | 2.03521
| 1/32 | 1/32 | $8.3169\times 10^{-6}$ | $1.7370\times 10^{-6}$ | 2.04651 | 2.01042
1.8 | 1/4 | 1/4 | $1.1100\times 10^{-4}$ | $1.6543\times 10^{-5}$ | 2.06253 | 2.04523
| 1/8 | 1/8 | $2.6573\times 10^{-5}$ | $4.0081\times 10^{-6}$ | 2.05568 | 2.03435
| 1/16 | 1/16 | $6.3917\times 10^{-6}$ | $9.7845\times 10^{-7}$ | 2.03658 | 2.00525
| 1/32 | 1/32 | $1.5579\times 10^{-6}$ | $2.4372\times 10^{-7}$ | 2.02546 | 2.00358

Figure 1: Example (5.1): Convergence and Efficiency Plots of ETD-RDP-FEM vs
CN-FEM Schemes for Different Values of $\alpha$.

Figure 2: Example 5.1: Exact Solution vs ETD-RDP-FEM vs CN-FEM Schemes for
$\alpha=1.6$.

It is not surprising that both ETD-RDP-FEM and CN-FEM methods perform well in
this case as this is a linear model. However, the efficiency of CN-FEM method
declines as the fractional order $\alpha$ increases as can be seen in Table 1
and Figures 1 and 2. Overall, a second order convergence is attained by both
methods as well.

######  Example 5.2 (Nonlinear Model).

Here, we investigate the performance of the proposed ETD-RDP-FEM scheme for
the following nonlinear Riesz-tempered fractional reaction diffusion equation:

$\begin{split}u_{t}(x,t)-\partial_{\mid
x\mid}^{\alpha,\lambda}u(x,t)&=f(x,t,u)\\\ u(0,t)&=u(1,t)=0,\\\
u(x,0)&=g(x),\end{split}$ (64)

where

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

and

$\displaystyle f(x,t,u)$ $\displaystyle=$ $\displaystyle
u^{2}-x^{2}(1-x)^{2}e^{-t}-x^{4}(1-x)^{4}e^{-2t}$
$\displaystyle+\frac{e^{-t}}{2cos(\alpha\pi/2)}\bigg{[}e^{-\lambda
x}H_{k,m,\lambda}(x)+e^{-\lambda(x+1)}\bar{H}_{k,m,\lambda}(x)-2\lambda^{\alpha}x^{2}(1-x)^{2}\bigg{]},$

with

$\displaystyle H_{k,m,\lambda}(x)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\sum_{m=2}^{4}\frac{\lambda^{k}\Gamma(k+m+1)A_{m}}{\Gamma(k+1)\Gamma(m+1-\alpha))}x^{k+m},$
$\displaystyle\bar{H}_{k,m,\lambda}(x)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\sum_{m=2}^{4}\frac{\lambda^{k}\Gamma(k+m+1)A_{m}}{\Gamma(k+1)\Gamma(m+1-\alpha))}(1-x)^{k+m-\alpha},$

and $A_{2}=1$, $A_{3}=-2$, and $A_{4}=1$.
The nonlinearity stems from the polynomial $u^{2}$ term in the function $f$.
The numerical solution from both schemes agree very well with the exact
solution (see Figure 4). The second order convergence of the ETD-RDP-FEM
scheme is evident in Table 2 for all values of $\alpha$. The scheme is also
more accurate and significantly faster than the existing CN-FEM scheme (see
Figure 3).

Table 2: Example (5.2): Time rate of convergence and $L_{\infty}$-norm of ETD-
RDP-FEM and CN-FEM with different values of $\alpha$ at fixed $t=1$.

.  Norm Order $\alpha$ $h$ $\tau$ CN-FEM ETD-RDP-FEM CN-FEM ETD-RDP-FEM 1.2
1/4 1/4 $7.5000\times 10^{-3}$ $4.5000\times 10^{-3}$ 2.18640 2.09653 1/8 1/8
$1.6000\times 10^{-3}$ $1.1000\times 10^{-3}$ 2.09254 2.07329 1/16 1/16
$3.7515\times 10^{-4}$ $2.6138\times 10^{-4}$ 2.08117 2.04752 1/32 1/32
$8.8656\times 10^{-5}$ $6.3228\times 10^{-5}$ 2.05712 2.02367 1.4 1/4 1/4
$5.5000\times 10^{-3}$ $6.1000\times 10^{-3}$ 2.12561 2.08561 1/8 1/8
$1.7000\times 10^{-3}$ $1.30000\times 10^{-3}$ 2.11021 2.07142 1/16 1/16
$3.9374\times 10^{-4}$ $3.0930\times 10^{-4}$ 2.08141 2.04653 1/32 1/32
$9.3034\times 10^{-5}$ $7.4871\times 10^{-5}$ 2.0575 2.02142 1.6 1/4 1/4
$3.1000\times 10^{-3}$ $9.0000\times 10^{-4}$ 2.08589 2.07125 1/8 1/8
$7.3021\times 10^{-4}$ $2.1416\times 10^{-4}$ 2.07661 2.05415 1/16 1/16
$1.7311\times 10^{-4}$ $5.1568\times 10^{-5}$ 2.06910 2.03536 1/32 1/32
$4.1254\times 10^{-5}$ $1.2580\times 10^{-5}$ 2.04113 2.01250 1.8 1/4 1/4
$5.5631\times 10^{-4}$ $1.6328\times 10^{-4}$ 2.06862 2.05004 1/8 1/8
$1.3262\times 10^{-4}$ $3.9428\times 10^{-5}$ 2.05108 2.04129 1/16 1/16
$3.2002\times 10^{-5}$ $9.5789\times 10^{-6}$ 2.05056 2.02258 1/32 1/32
$7.7250\times 10^{-6}$ $2.3575\times 10^{-6}$ 2.03401 2.01961

Figure 3: Example (5.2): Convergence and Efficiency Plots of ETD-RDP-FEM vs
CN-FEM Schemes

Figure 4: Example 5.2: Exact Solution vs ETD-RDP-FEM vs CN-FEM Schemes

######  Example 5.3 (Model with Nonsmooth Data).

Consider the non-smooth initial-boundary value problem in the following
tempered fractional diffusion equation:

$\begin{split}u_{t}(x,t)-\partial_{\mid
x\mid}^{\alpha,\lambda}u(x,t)&=f(x,t),\\\ u(0,t)&=u(1,t)=0,\\\ \end{split}$
(65)

with the initial condition

$u(x,0)=\left\\{\begin{array}[]{ll}0\ \ \quad 0\leq x<0.25,\\\ 1\ \ \quad
0.25\leq x<0.75,\\\ 0\ \ \quad 0.75\leq x\leq 1.\end{array}\right.$ (66)

Figure 5: Plot of the piecewise function $u(x,0)$.

The source term, $f(x,t)$ here is the same as defined in Equation 5.1. Note
that Example 5.3 has a nonsmooth initial data which could lead to spurious
oscillation. It is worthwhile to understand how the both methods considered
here handle a problem with such initial condition. We use a finer grid
($h=1/512$) for our reference solution in this case since there is no known
exact solution.

Table 3: Comparing order of convergence and $L_{\infty}$ norm for ETD-RDP-FEM
and CN-FEM for $h=1/512$ for Experiment (5.3)

.  Norm Order $\alpha$ $\tau$ CN-FEM ETD-RDP-FEM CN-FEM ETD-RDP-FEM 1.2 1/4
6.20970$\times 10^{-5}$ 2.85963$\times 10^{-5}$ 2.01558 1.83061 1/8
1.53575$\times 10^{-5}$ 8.0398$\times 10^{-6}$ 2.02613 2.09858 1/16
4.7408$\times 10^{-6}$ 1.8772$\times 10^{-6}$ 2.08905 2.01992 1/32
1.1143$\times 10^{-6}$ 5.713$\times 10^{-7}$ 1.9765 2.05117 1.4 1/4
1.21690$\times 10^{-4}$ 1.15635$\times 10^{-4}$ 1.73515 2.01515 1/8
3.65528$\times 10^{-5}$ 2.8607$\times 10^{-5}$ 1.99108 2.20616 1/16
9.1949$\times 10^{-6}$ 7.1194$\times 10^{-6}$ 1.89055 2.00655 1/32
2.8074$\times 10^{-6}$ 1.7534$\times 10^{-6}$ 2.10271 2.02159 1.6 1/4
8.59655$\times 10^{-4}$ 2.15255$\times 10^{-4}$ 1.98175 2.02479 1/8
2.1765$\times 10^{-4}$ 5.2897$\times 10^{-5}$ 2.08775 2.04210 1/16
5.1202$\times 10^{-5}$ 1.2903$\times 10^{-5}$ 2.09112 2.03548 1/32
1.2017$\times 10^{-5}$ 3.2220$\times 10^{-6}$ 2.21559 2.00170 1.8 1/4
3.45616$\times 10^{-4}$ 1.01124$\times 10^{-4}$ 1.99570 2.11327 1/8
8.6662$\times 10^{-5}$ 2.3372$\times 10^{-5}$ 1.93860 2.02320 1/16
2.2608$\times 10^{-5}$ 5.7985$\times 10^{-6}$ 2.13117 2.00015 1/32
5.1608$\times 10^{-6}$ 1.4495$\times 10^{-6}$ 2.10124 2.00100

Figure 6: Example 5.3: ETD-RDP-FEM vs CN-FEM Schemes

Figure 7: Example (5.3): Convergence and Efficiency Plots of ETD-RDP-FEM vs
CN-FEM Schemes.

We observe from the above results in Example 5.3 also that second-order
convergence is achieved by both ETD-RDP-FEM and CN-FEM, see 3. However, ETD-
RDP-FEM outperformed CN-FEM across different values of $\alpha$ as well here.
This observation agrees with what has been established in the literature
regarding CN method.

###  5.2 Discussion

The focus here is to discuss the results obtained from the above Examples
(5.1)-(5.3) in the previous subsection especially the impact of the fractional
order $\alpha$ on the performance of numerical schemes. In particular, we
discuss the effect of $\alpha$ on CPU time and accuracy. This is very
important to explore in order to gain insight into the impact of the non-
integer order on the considered problems in the previous Section.

####  5.2.1 Effect of Fractional Order $\alpha$ on CPU Time

It is widely known that the efficiency of a method could be greatly affected
depending on the values of fractional order $alpha$. The impact is usually
noticed in the CPU time. We have performed experiments across different values
of $alpha$ here to understand the effect on the efficiency of the methods of
interest here.

Table 4: Comparing Efficiency of ETD-RDP-FEM and CN-FEM. The time step $\tau=1/32$ and mesh size $h=1/512$ |  CPU time (in secs)
---|---
$\alpha$ | Methods | Linear | Nonlinear | Nonsmooth
1.2 | ETD-RDP-FEM | 1.687328 | 8.415730 | 0.285885
| CN-FEM | 2.783019 | 17.379510 | 0.506686
1.4 | ETD-RDP-FEM | 1.791548 | 8.204377 | 0.319528
| CN-FEM | 2.862972 | 20.911838 | 0.532703
1.6 | ETD-RDP-FEM | 1.703123 | 8.018364 | 0.270826
| CN-FEM | 2.881956 | 24.915706 | 0.520299
1.8 | ETD-RDP-FEM | 1.833103 | 8.377576 | 0.266322
| CN-FEM | 3.047225 | 33.042088 | 0.491660

$(a)$ Linear Model

$(b)$ Nonlinear Model

$(c)$ Nonsmooth Model

Figure 8: Computational Efficiency Plots for Examples (5.1) - (5.3) Table 5: Percentage Improvement of ETD-RDP-FEM over CN-FEM $\alpha$ | Nonsmooth | Nonlinear | Linear
---|---|---|---
1.2 | 43.57% | 51.57% | 39.42%
1.4 | 40.02% | 60.78% | 37.41%
1.6 | 47.97% | 67.81% | 40.86%
1.8 | 45.84% | 72.10% | 39.86%
2 | 44.47% | 75.61% | -
Figure 9: Plot of the percentage improvements in CPU time of ETD-RDP-FEM over
CN-FEM across different Models and fractional order $\alpha$.

The comparison of the ETD-RDP-FEM with CN-FEM is computationally more
efficient for different values of fractional orders in all the Examples
considered as can be seen in Table 4. It is notable to see that the efficiency
of both the ETD-RDP-FEM and the CN-FEM in the linear and nonsmooth models is
approximately linear, see Figures 8 $(a)$ and $(c)$. However, in the case of
nonlinear model, CN-FEM is much slower than ETD-RDP-FEM as shown in Figure
8$(b)$. This is expected as ETD-RDP-FEM requires no nonlinear solver unlike
CN-FEM. Newton method is used to handle the nonlinear term here with
predefined tolerance $10^{-6}$ which is fixed across different values of
$alpha$. We further demonstrate the superiority of efficiency of ETD-RDP-FEM
over CN-FEM in Table 5 and Figure 9. The difference is efficiency is calculate
in percentage across different values of $a\alpha$. On average, ETD-RDP-FEM
shows notable improvement over CN-FEM of approximately $47.97\%$, $67.81\%$,
and $40.86\%$ in the linear, nonlinear and nonsmooth cases, respectively.

####  5.2.2 Effect of Fractional Order $\alpha$ on Accuracy

Our numerical experiments for all test problems show that ETD-RDP-FEM produces
more accurate solutions than the existing CN-FEM scheme. Comparing the errors
for each scheme across different values of the fractional order $\alpha$, we
observe from Table 6 that numerical errors tends to increase as $\alpha$
approaches 1 and decreases as $\alpha$ approaches 2 for the linear and
nonlinear model (Figure 10 a,b). This observation is consistent with the
theoretical results and performance of many numerical schemes for solving
fractional differential equations [36, 37, 38]. However, for the nonsmooth
problem the error increases with increasing alpha up to $\alpha=1.6$ and then
decreases. This is an interesting observation, worthy of further exploration
(Figure 10 c).

Table 6: Comparing the accuracy of ETD-RDP-FEM and CN-FEM across different fractional order $\alpha$. The time step $\tau=1/32$ and mesh size $h=1/512$ |  $L_{\infty}$ Error
---|---
$\alpha$ | Methods | Linear | Nonlinear | Nonsmooth
1.2 | ETD-RDP-FEM | $1.0640\times 10^{-5}$ | $6.7603\times 10^{-6}$ | $5.7130\times 10^{-7}$
| CN-FEM | $3.1163\times 10^{-5}$ | $4.7221\times 10^{-5}$ | $1.1143\times 10^{-6}$
1.4 | ETD-RDP-FEM | $6.5507\times 10^{-6}$ | $3.2633\times 10^{-6}$ | $1.7534\times 10^{-6}$
| CN-FEM | $1.4364\times 10^{-5}$ | $4.5148\times 10^{-5}$ | $2.8074\times 10^{-6}$
1.6 | ETD-RDP-FEM | $1.5818\times 10^{-6}$ | $1.4703\times 10^{-6}$ | $3.2200\times 10^{-6}$
| CN-FEM | $2.4090\times 10^{-6}$ | $4.4593\times 10^{-6}$ | $1.2017\times 10^{-5}$
1.8 | ETD-RDP-FEM | $1.0600\times 10^{-6}$ | $1.0243\times 10^{-6}$ | $1.4495\times 10^{-6}$
| CN-FEM | $1.7028\times 10^{-6}$ | $1.3288\times 10^{-6}$ | $5.1608\times 10^{-6}$

$(a)$ Linear Model

$(b)$ Nonlinear Model

$(c)$ Nonsmooth Model

Figure 10: Computational accuracy Plots for Examples (5.1) - (5.3)

##  6 Conclusion and Recommendation

In conclusion, the Fractional Diffusion Equation (FDE) serves as a
mathematical model that characterizes anomalous transport processes with non-
local and long-range dependencies, diverging from traditional diffusion
patterns. This makes its numerical solution challenging, as it requires
complex integral operators and incurs high computational expenses for precise
approximations. In this work, we introduced an Exponential Time Differencing
Finite Element Method (ETD-RDP-FEM) to efficiently address both linear and
nonlinear problems. Our method, which is grounded on a rational function with
distinct real poles to discretize matrix exponentials, leads to an L-stable
scheme. The proposed technique exhibits second-order convergence and
resiliency when applied to problems with non-smooth initial conditions,
underscoring its efficacy and adaptability in handling intricate scenarios.
Our findings indicate that the developed approach surpasses the Crank-Nicolson
technique in terms of computational efficiency, as demonstrated by reduced CPU
time and increased accuracy.

In the future, the ETD-FEM scheme may be a promising tool for solving more
complex problems, as it offers improved CPU efficiency compared to existing
numerical methods. Specifically, future work could focus on extending the ETD-
RDP-FEM scheme to solve fractional reaction diffusion equation with an
advection term, as well as using it to solve systems of reaction-diffusion-
advection equations. Additionally, the scheme’s potential applicability to
multidimensional problems in irregular domains should also be explored
further.

## References

  * [1] H. Sun, Y. Zhang, D. Baleanu, W. Chen W, Y. Chen, (2018), A new collection of real-world applications of fractional calculus in science and engineering, Commun Nonlinear Sci Numer Simulat, 64, 213-231.
  * [2] P. Kumar, D. Baleanu, V.S. Erturk, M. Inc, and V. Govindaraj, (2022), A delayed plant disease model with Caputo fractional derivatives. Advances in Continuous and Discrete Models, 1, 1-22.
  * [3] M. Partohaghighi, A. Akgül, L. Guran, and M.F. Bota, (2022), Novel Mathematical Modelling of Platelet-Poor Plasma Arising in a Blood Coagulation System with the Fractional Caputo–Fabrizio Derivative. Symmetry, 14(6), 1128.
  * [4] M.S. Hashemi, M. Partohaghighi, H. Ahmad, (2023), mathematical modelings of the human liver and hearing loss systems with fractional derivatives, International Journal of Biomathematics, 16 (01), 2250068.
  * [5] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, (2006), Theory and Applications of Fractional Differential Equations, 204, 1-521.
  * [6] I. Podlubny, (1998), Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications. Academic Press.
  * [7] D. Baleanu, K. Diethelm, E. Scalas, J.J. Trujillo, (2012), Fractional calculus: models and numerical methods, 7, 3, 1-383.
  * [8] O.S. Iyiola, F.D. Zaman, (2014), A fractional diffusion equation model for cancer tumor, AIP Advances, 4(10).
  * [9] F. Mainardi, (2010), Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models, World Scientific, 1-368.
  * [10] R. Metzler and J. Klafter, (2000), The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Physics Reports, 339(1), 1-77.
  * [11] K.M. Furati, O.S. Iyiola, M. Kirane, (2014), An inverse problem for a generalized fractional diffusion, Applied Mathematics and Computation, 249, 24-31.
  * [12] M. Etefa, G. Guerekata, P. Ngnepieba, O.S. Iyiola, (2023), On a generalized fractional differential Cauchy problem, Malaya Journal of Matematik, 11(1), 80-93.
  * [13] H. Sun, W. Chen, and Y. Chen, (2021), Modeling anomalous diffusion in biological tissues using a fractional-order reaction-diffusion model. Journal of Mathematical Biology, 82(4), 42.
  * [14] C. Ionescu, A. Lopes, D. Copot, J. Machado and J. Bates, (2017), The role of fractional calculus in modeling biological phenomena: A review, Communications in Nonlinear Science and Numerical Simulation, 51, 141-159.
  * [15] O.S. Iyiola, (2017), Exponential integrator methods for nonlinear fractional reaction-diffusion models, University of WIsconsin, WI.
  * [16] J. Liang, Y. Li, and D. Baleanu, (2021), Time-fractional mathematical models for anomalous transport phenomena in biological tissues. Communications in Nonlinear Science and Numerical Simulation, 93, 105498.
  * [17] F. Sabzikar, M.M. Meerschaert, J. Chen, (2015), Tempered fractional calculus, Journal of Computational Physics 293, 14-28.
  * [18] M.M. Meerschaert, F. Sabzikar (2014), Stochastic integration for tempered fractional Brownian motion, Stochastic Processes and their Applications 124, 2363-2387.
  * [19] C. Celik and M. Duman, (2017), Finite element method for a symmetric tempered fractional diffusion equation, Applied Numerical Mathematics, 120, 270–286.
  * [20] X. Guo, Y. Li, H. Wang, (2018), A high order finite difference method for tempered fractional diffusion equations with applications to the CGMY model, SIAM Journal on Scientific Computing, 40,5, A3322-A3343.
  * [21] E.O. Asante-Asamani, (2016), An exponential time differencing scheme with a real distinct poles rational function for advection-diffusion-reaction systems.
  * [22] E.O. Asante-Asamani, A.Q.M. Khaliq and B.A. Wade, (2016), A real distinct poles exponential time differencing scheme for reaction–diffusion systems, Journal of Computational and Applied Mathematics 299, 19, 24-34.
  * [23] O.S. Iyiola, E.O. Asante-Asamani, K.M Furati, A.Q,M. Khaliq and B.A. Wade, (2018), Efficient time discretization scheme for nonlinear space fractional reaction–diffusion equations, International Journal of Computer Mathematics 95, 9, 1274-1291.
  * [24] E.O. Asante-Asamani, A. Kleefeld and B.A. Wade, (2020), A second-order exponential time differencing scheme for non-linear reaction-diffusion systems with dimensional splitting, Journal of computational physics, 415, 109490, 6.
  * [25] O.S. Iyiola and B.A. Wade, (2018), Exponential integrator methods for systems of non-linear space-fractional models with super-diffusion processes in pattern formation, Computers & Mathematics with Applications 75, 12, 3719-3736.
  * [26] Abramowitz, M., Stegun, I. A. (1972). Handbook of mathematical functions: with formulas, graphs, and mathematical tables. US Government Printing Office.
  * [27] Abdon Atangana, Fractional Operators and Their Applications, R. Garra, R. Gorenflo, F. Polito, Z. Tomovski, Hilfer-Prabhakar derivatives and some applications. Applied Mathematics and Computation 242 (2014), 576–589.
  * [28] B. Jin, R. Lazarov, Z. Zhou, (2016), A Petrov–Galerkin finite element method for fractional convection-diffusion equations, SIAM Journal on Numerical Analysis, 54,1, 10.1137/140992278.
  * [29] Y. Lian, Y. Ying, S. Tang, S. Lin, G. Wagner, (2016), A Petrov–Galerkin finite element method for the fractional advection–diffusion equation, Computer Methods in Applied Mechanics and Engineering, 309, 388-410.
  * [30] L. Li, D. Xu, M. Luo, (2013), Alternating direction implicit Galerkin finite element method for the two-dimensional fractional diffusion-wave equation, Journal of Computational Physics, 255, 471-485.
  * [31] S.M.Cox and P.C.Matthews, (2002), Exponential Time Differencing for Stiff Systems, Journal of Computational Physics, 176, 430-455.
  * [32] Kassam, Aly-Khan and Trefethen, Lloyd N, (2005), Fourth-order time-stepping for stiff PDEs, SIAM Journal on Scientific Computing, 26, 1214-1233.
  * [33] U. Ali, F. Abdullah, A. Ismail, (2017), Crank-Nicolson finite difference method for two-dimensional fractional sub-diffusion equation, Journal of Interpolation and Approximation in Scientific Computing, 2, 18-29.
  * [34] N. Sweilam, H. Moharram, N. Moniem, (2014), A parallel Crank–Nicolson finite difference method for time-fractional parabolic equation, Journal of Numerical Mathematics, 22, 4, 363–382.
  * [35] R.J. Leveque, (2007), Finite difference methods for ordinary and partial differential equations, Society for Industrial and Applied Mathematics, 1-356.
  * [36] Y.O. Afolabi, T.A. Biala, O.S. Iyiola, A.QM Khaliq, B.A. Wade, (2022), A Second-Order Crank-Nicolson-Type scheme for nonlinear space-time reaction–diffusion equations on time-graded meshes, Fractal and Fractional, 7(1), 40.
  * [37] K. Wang, Z. Zhou, (2020), High-order time-stepping schemes for semilinear subdiffusion equations, SIAM J. Numer. Anal., 58, 3226-3250.
  * [38] K. Mustapha, (2020), An L1 approximation for a fractional reaction-diffusion equation, a second-order error analysis over time-graded meshes. SIAM J. Numer. Anal., 58, 1319-1338.