Analytical and numerical investigation on the tempered time-fractional operator with application to the Bloch equation and the two-layered problem

Abstract

In the continuous-time random walk model, the time-fractional operator usually expresses an infinite waiting time probability density. Different from that usual setting, this work considers the tempered time-fractional operator, which reflects a finite waiting time probability density. Firstly, we analyse the solution of a tempered benchmark problem, which shows a weak singularity near the initial time. The L1 scheme on graded mesh and the weighted shifted Grünwald–Letnikov formula with correction terms are adapted to deal with the non-smooth solution, in which we compare these two methods systematically in terms of the convergence and consumed CPU time. Furthermore, a fast calculation for the time tempered Caputo fractional derivative is developed based on a sum-of-exponentials approximation, which significantly reduces the running time. Moreover, the tempered operator is applied to the Bloch equation in nuclear magnetic resonance and a two-layered problem with composite material exhibiting distinct memory effects, for which both the analytical (or semi-analytical) and numerical solutions are derived using transform techniques and finite difference methods. Data-fitting results verify that the tempered time-fractional model is much effective to describe the MRI data. An important finding is that, compared with the fractional index, the tempered operator parameter could further accelerate the diffusion. The tempered model with two parameters $\alpha$ and $\rho$ are more flexible, which can avoid choosing a too small fractional index leading to low regularity and strong heterogeneity.

1 Introduction

We start with a continuous-time random walk (CTRW) model, which uses a waiting time probability density function (pdf) $\psi (t)$ and a jump length pdf $\eta (x)$ to describe the particles having different jumps ($x_{1},x_{2},\ldots ,x_{i},\ldots$) at different times ($t_{1},t_{2},\ldots ,t_{i},\ldots$). Define a jump pdf $\Phi (x,t)$, then the jump length pdf and the waiting time pdf can be expressed as [1, 2]:

$$\begin{aligned} \eta (x)=\int _{0}^{\infty }\Phi (x,t)dt,\quad \psi (t)=\int _{-\infty }^{\infty }\Phi (x,t)dx. \end{aligned}$$

We are interested in the case that the jump length is independent of the waiting time, which leads to $\Phi (x,t)=\psi (t)\eta (x)$. For a given $\psi (t)$ and $\eta (x)$, the probability of finding a particle at position x and time t satisfies the Montroll–Weiss (MW) equation [1, 2], which has the following form in the Fourier–Laplace space:

$$\begin{aligned} \widehat{{\overline{P}}}(\omega ,s)=\frac{1-{\overline{\psi }}(s)}{s}\frac{1}{1-{\overline{\psi }}(s){\hat{\eta }}(\omega )}, \end{aligned}$$

where ${\overline{\psi }}(s)=\int _{0}^{\infty }e^{st}\psi (t)dt$ and ${\hat{\eta }}(\omega )=\int _{-\infty }^{\infty }e^{i\omega x}\eta (x)dx$. When a long-tailed waiting time pdf with the asymptotic behaviour [2] $\psi (t)\sim A_{\alpha }(C_{\tau }/t)^{1+\alpha }$, $0<\alpha <1$, is considered, the Caputo time-fractional operator ${_{0}^{C}D_{t}^{\alpha }}$ can be introduced in the setting. In this case, the corresponding characteristic waiting time is:

$$\begin{aligned} {T=\int _{0}^{\infty }t\psi (t)dt \sim A_{\alpha }C_{\tau }^{1+\alpha }\int _{0}^{\infty }t^{-\alpha }dt=\infty }. \end{aligned}$$

In the real world, since the waiting times may not be arbitrarily long with a pure power-law distribution, the concept of tempered waiting time pdf was proposed [3, 4], which is expected to have broader applicability [5]. Now consider a tempered long-tailed waiting time pdf [6] $\psi (t)\sim A_{\alpha }(C_{\tau }/t)^{1+\alpha }e^{-\rho t}$ with Laplace transform ${\overline{\psi }}(s)\sim 1-(C_{\tau }(s+\rho ))^{\alpha }$, then the Caputo-tempered time-fractional operator ${_{0}^{C}D_{t}^{(\alpha ,\rho )}}$ can be introduced as

$$\begin{aligned}&{_{0}^{C}D_{t}^{(\alpha ,\rho )}}u(t)=e^{-\rho t}{_{0}^{C}D_{t}^{\alpha }} (e^{\rho t}u(t))\\&\quad =\frac{e^{-\rho t}}{\Gamma (1-\alpha )}\int _{0}^{t}\frac{ d(e^{\rho s}u(s))}{ds}\frac{1}{(t-s)^{\alpha }}ds. \end{aligned}$$

The corresponding characteristic waiting time turns out to be

$$\begin{aligned}&T=\int _{0}^{\infty }t\psi (t)dt \sim A_{\alpha }C_{\tau }^{1+\alpha }\int _{0}^{\infty }t^{-\alpha }e^{-\rho t}dt\\&\quad =A_{\alpha }C_{\tau }^{1+\alpha }\frac{\Gamma (1-\alpha )}{\rho ^{1-\alpha }}, \end{aligned}$$

which is finite. Motivated by this, we investigate the Caputo-tempered fractional operator in this paper.

In the past two decades, some attention has been paid to the investigation of the tempered operator. In [3, 4], the truncated Lévy process was introduced to eliminate the arbitrarily large flights produced by Lévy stable distributions. In [7], Cartea and del-Castillo-Negrete considered the CTRWs with exponentially truncated Lévy jump pdf, which led to a transport equation with space-tempered fractional derivatives to describe the interaction between long jumps and memory in the intermediate asymptotic regime. In [5], Meerschaert et al. proposed a novel tempered anomalous diffusion model with exponentially tempered waiting times to capture the natural cut-off of retention times in heterogeneous systems. In [8], Zhang et al. presented a tempered fractional mobile–immobile model with a tempered stable pdf for both particle jump length and resting time able to explain the sediment transport involving super-diffusive, sub-diffusive, and regular diffusive dynamics. In [9], Wu et al. derived the tempered fractional Feynman–Kac equation to describe the functional distribution of the paths of tempered anomalous dynamics. In [10], Boniece et al. constructed two classes of second-order, non-Gaussian transient anomalous diffusion models to depict the tempered fractional Lévy process.

In addition, different numerical methods were proposed to deal with the tempered fractional models [11,12,13,14], most of which are based on the spatial tempered operator. Here, we concentrate on a numerical treatment of the tempered time-fractional model. In [15], Hanert and Piret applied a Chebyshev pseudospectral method to discretise a time-space-tempered fractional diffusion equation. In [16], Chen and Deng proposed some high-order finite difference schemes for the time-tempered fractional Feynman–Kac equation. In [17], Ding and Li developed a tempered fractional-compact difference formula to deal with a time Caputo-tempered partial differential equation. In [18], Cao et al. considered the finite element method for a tempered time fractional advection–dispersion equation.

Most existing numerical methods for the time-tempered models assume a smooth solution, while the works on numerical methods dealing with non-smooth solutions are sparse [19, 20]. The time-tempered benchmark problem shows that the solution has a weak singularity near the initial time, which motivates us to develop numerical methods to handle non-smooth solutions. To observe the behaviour of the tempered operator, we adapt the tempered operator to the Bloch equation and the two-layered problem. The main contributions of the present work are as follows:

• Beginning with a tempered benchmark problem, the exact expression of the solution is derived, which incorporates the Mittag–Leffler function and an exponentially tempered factor. When $t\rightarrow 0$, the solution shows weak singularity near the initial time, for which most existing high-order numerical schemes would fail. Based on the regularity of the solution to the benchmark problem, a modified L1 scheme on graded mesh and the weighted shifted Grünwald–Letnikov (WSGL) formula with correction terms are developed, for which a systematical comparison in terms of the convergence and consumed CPU time is presented. It finds that for the L1 scheme on graded mesh, changing the regularity index $\alpha$ does not affect the total consumed CPU time too much, which facilitates the implementation of the fast algorithm. The WSGL scheme with correction terms shows high accuracy and needs more correction terms when $\alpha$ is small, which increases the consumed CPU time and leads to an ill-conditioned matrix. A fast calculation for the time-tempered Caputo fractional derivative is developed based on a sum-of-exponentials approximation, in which the kernel function $t^{-\beta }$ is approximated under desired precision $\varepsilon$ using three different quadratures. Numerical results demonstrate that this fast method can reduce the running time significantly.

• For the tempered Bloch equation, the analytical solution is presented using Laplace transform. The numerical solution is also deduced, by which the effects of the fractional index $\alpha$ and tempered parameter $\rho$ are analysed in detail. The fractional index delays the ‘longitudinal’ relaxation and promotes the ‘transversal’ relaxation. In comparison, the tempered parameter accelerates the ‘longitudinal’ relaxation and alters the asymptotic value and further advances the ‘transversal’ relaxation. In addition, the classical monoexponential model, the time-fractional model and the tempered time-fractional model are compared to fit the MRI signal data in the human brain, of which the tempered time-fractional model performs the best.

• For the tempered diffusion problem, the associated mean squared displacement has one more factor $e^{-\rho t}$ compared to that of the pure time-fractional diffusion problem. The analytical solution using finite Fourier transform and numerical solution with stability and convergence analysis are applied to solve the problem, which is also extended to the two-layered problem composing two different materials. An important finding is that, compared with the fractional index, the tempered parameter could further accelerate the diffusion. The tempered model with two parameters $\alpha$ and $\rho$ are more flexible, which can avoid choosing a too small fractional index leading to low regularity and strong heterogeneity.

The structure of this paper is as follows: In Sect. 2, we consider a tempered benchmark problem, in which two classical numerical schemes are applied and the fast calculation for the time tempered Caputo fractional derivative is developed. In Sect. 3, the tempered Bloch equation is proposed, which is solved analytically and numerically. In Sect. 4, the tempered diffusion problem is presented, for which the analytical and numerical solutions are established. In Sect. 5, a tempered two-layered problem is introduced. Semi-analytical solution and numerical solution are derived. Finally, some conclusions are drawn in Sect. 6.

2 Benchmark problems

In this work, we denote C as a general constant independent of the grid size and may take distinct values in different contexts. Firstly, we start with the following benchmark problem:

$$\begin{aligned}&{_0^CD_t^{(\alpha ,\rho )}}u(t)=-k_0 u(t),~0<t\le T,\nonumber \\&u(0)=u_0,\quad 0<\alpha <1,\quad \rho \ge 0, \end{aligned}$$

(1)

where $k_{0}$ and $u_{0}$ are some constants. If $\rho =0$, problem (1) reduces to the fractional benchmark problem [21].

2.1 Analytical solution

Define the Laplace transform ${\mathcal {L}}\left\{ f(t)\right\} ={\bar{f}}(s)=\int _{0}^{\infty }e^{-st}f(t)dt$ and the inverse Laplace transform $f(t)={\mathcal {L}}^{-1}\left\{ {\bar{f}}(s)\right\} =\frac{1}{2\pi i}\int _{c-i\infty }^{c+i\infty }e^{st}{\bar{f}}(s)ds$. Applying the Laplace transform to (1) using the formula ${\mathcal {L}}\left\{ {_{a}^{C}D_{t}^{(\alpha ,\rho )}}u(t)\right\} =(s+\rho )^{\alpha }{\bar{u}}(s)-(s+\rho )^{\alpha -1}u(0)$, we have

$$\begin{aligned} (s+\rho )^{\alpha }{\bar{u}}(s)-(s+\rho )^{\alpha -1}u(0)=-k_{0}{\bar{u}}(s), \end{aligned}$$

which leads to

$$\begin{aligned} {\bar{u}}(s)=\frac{(s+\rho )^{\alpha -1}}{(s+\rho )^{\alpha }+k_{0}}u_{0}. \end{aligned}$$

Utilising the property ${\mathcal {L}}\left\{ e^{\rho t}f(t)\right\} ={\bar{f}}(s-\rho )$ and ${\mathcal {L}}\left\{ E_{\alpha }(-k_{0}t^{\alpha })\right\} =\frac{s^{\alpha -1}}{s^{\alpha }+k_{0}}$, we can derive

$$\begin{aligned} u(t)=u_{0}e^{-\rho t}E_{\alpha }(-k_{0}t^{\alpha }), \end{aligned}$$

where $E_{\alpha }(z)=\sum _{n=0}^{\infty }\frac{z^{n}}{\Gamma (n\alpha +1)}$ is the Mittag–Leffler function. When $t\rightarrow 0$, using the Taylor series $e^{-\rho t}=\sum _{j=0}^{\infty }\frac{(-\rho t)^{j}}{j!}$, we obtain

$$\begin{aligned} u(t)&=u_0e^{-\rho t}E_{\alpha }(-k_0 t^\alpha )\nonumber \\&=u_0\sum _{j=0}^{\infty }\frac{(-\rho t)^j}{j!}\sum _{n=0}^{\infty }\frac{(-k_0t^\alpha )^n}{\Gamma (n\alpha +1)}\nonumber \\&=u_0\sum _{j=0}^{\infty }\sum _{n=0}^{\infty }C_{nj}t^{j+n\alpha }. \end{aligned}$$

(2)

We can see that, similar to the fractional benchmark problem, the solution shows a weak singularity near the initial time $t=0$ since $\frac{du(t)}{dt}$ blows up as $t\rightarrow 0$. We will extend two classical methods to deal with the weak singularity problem: the L1 scheme on graded mesh and the WSGL scheme with correction terms.

2.2 L1 scheme on graded mesh

Let N be a positive integer. Denote $t_{n}=T\left( \frac{n}{N}\right) ^{r}$, $r\ge 1$, $n=0,1,2,\ldots ,N$, $\tau _{n}=t_{n}-t_{n-1}$, $n=1,2,\ldots ,N$. When $r=1$, the mesh recovers the uniform mesh. Recall the L1 formula on graded mesh [22] to approximate the Caputo derivative at $t=t_{n}$:

$$\begin{aligned}&{_{0}^{C}D_{t}^{\alpha }}u(t_{n})=\frac{1}{\Gamma (2-\alpha )}\sum _{k=0}^{n-1}\frac{u(t_{k+1})-u(t_{k})}{\tau _{k+1}}\left[ (t_{n}-t_{k})^{1-\alpha }\right. \\&\quad \left. -(t_{n}-t_{k+1})^{1-\alpha }\right] +R_{1}^{n}. \end{aligned}$$

According to the definition of the Caputo-tempered fractional derivative, it is straightforward to derive

$$\begin{aligned}&{_{0}^{C}D_{t}^{(\alpha ,\rho )}}u(t_{n})=e^{-\rho t_{n}}{_{0}^{C}D_{t}^{\alpha }}(e^{\rho t_{n}}u(t_{n}))\\&\quad =\frac{\tau _{n}^{-\alpha }}{\Gamma (2-\alpha )}u(t_{n})-\frac{\tau _{n}^{-\alpha }}{\Gamma (2-\alpha )} e^{-\rho \tau _{n}}u(t_{n-1}) \\&\qquad + \sum _{k=0}^{n-2}\frac{e^{-\rho (t_{n}-t_{k+1})}u(t_{k+1})-e^{-\rho (t_{n}-t_{k})}u(t_{k})}{\Gamma (2-\alpha )\tau _{k+1}}\left[ (t_{n}-t_{k})^{1-\alpha }\right. \\&\qquad \left. -(t_{n}-t_{k+1})^{1-\alpha }\right] +R_{1}^{n} \\&:= {_{0}^{C}{\mathbb {D}}_{t}^{(\alpha ,\rho )}}u(t_{n})+R_{1}^{n}. \end{aligned}$$

For the bound of the error term $R_{1}^{n}$, we have the following lemma.

Lemma 1

Suppose $|\frac{du}{dt}|\le C(1+t^{\alpha -l})$ for $l=0,1,2$ and $0<\alpha <1$, then the truncation error of the L1 scheme satisfies

$$\begin{aligned} |R_1^n|&=\left| {_0^CD_t^{(\alpha ,\rho )}}u(t_n)- {_0^C{\mathbb {D}}_t^{(\alpha ,\rho )}}u(t_n)\right| \nonumber \\&\le C n^{-\min \{2-\alpha ,r\alpha \}}. \end{aligned}$$

(3)

Proof

Using the definition of the derivative and integration by parts, we have

$$\begin{aligned}&{_0^CD_t^{(\alpha ,\rho )}}u(t_n)- {_0^C{\mathbb {D}}_t^{(\alpha ,\rho )}}u(t_n)\\&\quad =\frac{e^{-\rho t_n}}{\Gamma (1-\alpha )}\sum _{k=0}^{n-1}\int _{t_k}^{t_{k+1}}\frac{ \left( e^{\rho s}[u(s)-u_h(s)]\right) '}{(t_n-s)^{\alpha }}ds\\&\quad =\frac{-\alpha e^{-\rho t_n}}{\Gamma (1-\alpha )}\sum _{k=0}^{n-1}\int _{t_k}^{t_{k+1}}(t_n-s)^{-(\alpha +1)}\\&\qquad e^{\rho s}[u(s)-u_h(s)]ds\\&\quad =\sum _{k=0}^{n-1}T_{n,k}, \end{aligned}$$

where $u_h(t)$ is a linear interpolation function to approximate u(t) in $[t_k,t_{k+1}]$ and

$$\begin{aligned} T_{n,k}=\frac{-\alpha e^{-\rho t_n}}{\Gamma (1-\alpha )}\int _{t_k}^{t_{k+1}}(t_n-s)^{-(\alpha +1)} e^{\rho s}[u(s)-u_h(s)]ds. \end{aligned}$$

As $u(s)-u_h(s)=\frac{1}{2}u''(\xi _k)(s-t_{k-1})(t_k-s)$, $\xi _k\in [t_{k-1},t_k]$, we have

$$\begin{aligned} |T_{n,k}|&\le C\tau _{k+1}^2\max \limits _{t\in [t_k,t_{k+1}]}|u''(t)| \int _{t_k}^{t_{k+1}}\\&\quad \times (t_n-s)^{-(\alpha +1)}ds. \end{aligned}$$

Similar to the proof in Lemma 5.2 in [22], when $1\le k\le n-1$, we have

$$\begin{aligned}&\sum _{k=1}^{\lceil n/2\rceil -1}|T_{n,k}|\le \left\{ \begin{array}{ll} Cn^{-r(\alpha +1)},&{} \mathrm{{if}}~ r(\alpha +1)<2,\\ Cn^{-2}\ln n,&{} \mathrm{{if}}~ r(\alpha +1)=2,\\ Cn^{-2},&{} \mathrm{{if}}~ r(\alpha +1)>2, \end{array}\right. \\&\quad \text {and}\quad \sum _{k=\lceil n/2\rceil }^{n-2}|T_{n,k}|\le Cn^{-(2-\alpha )}. \end{aligned}$$

Table 1 The error and convergence order of the L1 scheme on graded mesh for the tempered ODE (1) at $t=1$, where the parameters are $k_0=2$, $\rho =0.5$, $u_0=1$

Now we consider the bound of $T_{n,0}$. If $n=1$, then

$$\begin{aligned} T_{1,0}&=\frac{\tau _1^{-\alpha }}{\Gamma (2-\alpha )}[F(t_1)-F(t_0)]\\&\quad -\frac{1}{\Gamma (1-\alpha )}\int _0^{t_1}(t_1-s)^{-\alpha }F'(s)ds, \end{aligned}$$

where $F(t)=e^{\rho t}u(t)$. As $F'(t)=e^{\rho t}(\rho u(t)+u'(t))$, then $|F'(t)|\le Ct^{\alpha -1}$, $t\in (0,t_1)$. Similar to the proof in [22], we can obtain $|T_{1,0}|\le C$. For $n>1$, we have

$$\begin{aligned} |T_{n,0}|&=\left| \frac{\alpha e^{-\rho t_n}}{\Gamma (1-\alpha )}\int _{t_0}^{t_{1}}(t_n-s)^{-(\alpha +1)} e^{\rho s}[u(s)-u_h(s)]ds\right| \\&\le C\int _{t_0}^{t_{1}}(t_n-s)^{-(\alpha +1)}|u(s)-u_h(s)| ds\\&=C\int _{t_0}^{t_{1}}(t_n-s)^{-(\alpha +1)}\left| \int _0^s(u(\theta )-u_h(\theta ))'d\theta \right| ds\\&\le C\int _{t_0}^{t_{1}}(t_n-s)^{-(\alpha +1)}\left( \int _0^s|u'(\theta )|d\theta \right. \\&\left. +\int _0^st_1^{-1}\int _0^{t_1}|u'(\eta )|d\eta d\theta \right) ds\\&\le C\int _{t_0}^{t_{1}}(t_n-s)^{-(\alpha +1)}(s^\alpha +st_1^{\alpha -1})ds\\&\le Ct_1^{\alpha }\int _{t_0}^{t_{1}}(t_n-s)^{-(\alpha +1)}ds\\&\le C \left( \frac{t_n-t_1}{t_1} \right) ^{-\alpha }\le Cn^{-r\alpha }. \end{aligned}$$

Finally, we consider the term $T_{n,n-1}$

$$\begin{aligned} |T_{n,n-1}|&\le C\tau _{n}^2\max \limits _{t\in [t_{n-1},t_{n}]}|u''(t)| \int _{t_{n-1}}^{t_{n}}(t_n-s)^{-(\alpha +1)}ds\\&\le C\tau _{n}^2 t_{n-1}^{\alpha -2}\tau _{n}^{-\alpha }\le C\tau _{n}^{2-\alpha } t_{n-1}^{\alpha -2} \le Cn^{-(2-\alpha )}. \end{aligned}$$

Combining these error bounds, we can derive (3). $\square$

Denote $u^{n}$ as the numerical approximation to $u(t_{n})$. We can derive the numerical scheme of (1) as

$$\begin{aligned} {_{0}^{C}{\mathbb {D}}_{t}^{(\alpha ,\rho )}}u^{n}=-k_{0}u^{n}. \end{aligned}$$

We give an example to illustrate the effectiveness of the numerical scheme by choosing $k_{0}=2$, $\rho =0.5$, $u_{0}=1$ and $T=1$. The error and convergence order of the scheme with varying N and $\alpha$ are presented in Table 1. It can be seen that the L1 scheme fails for a uniform mesh ($r=1$), while the highest possible convergence order $2-\alpha$ is obtained for the optimal graded mesh $r=\frac{2-\alpha }{\alpha }$, which means the L1 scheme is effective as well to approximate the time-tempered operator. As suggested by [22], it is better to choose $r=\frac{2(2-\alpha )}{\alpha }$ when the fractional index $\alpha$ is unknown.

2.3 The WSGL formula with correction terms

Generally, the analytical solution of the time-fractional differential equation contains the term $t^{\sigma _{m}}$, which exhibits low regularity when ${\sigma _{m}}$ is small. The WSGL formula with correction terms is a useful method to deal with the solution with low regularity term $t^{\sigma _{m}}$, in which the correction term can be exact or highly accurate to approximate the term $t^{\sigma _{m}}$ [23,24,25]. Even when the regularity indices in the correction terms do not match the singularity of the analytical solution exactly, the accuracy of the WSGL formula can still be improved significantly [24]. As shown in (2), the solution incorporates a similar term $t^{\sigma _{m}}$, therefore the WSGL formula could be used to handle the problem. Here we use a uniform grid in the temporal direction. Denote $t_{n}=n\tau$, $n=0,1,2,\ldots ,N$, where $\tau =\frac{T}{N}$ is the uniform temporal step and $N\in {\mathbb {Z}}^+$. Introduce the Riemann–Liouville fractional derivative

$$\begin{aligned} {_{0}^{R}D_{t}^{\alpha }}v(t)=\frac{1}{\Gamma (1-\alpha )}\frac{d}{dt} \int _{0}^{t}\frac{v(\xi )}{(t-\xi )^{\alpha }}d\xi ,\quad 0<\alpha <1. \end{aligned}$$

Then, the WSGL formula to discretise the Riemann–Liouville time-fractional derivative is

$$\begin{aligned} {_0^RD_t^{\alpha }}u(t_n)&=\tau ^{-\alpha }\sum _{k=0}^n\omega _{n-k}^{(\alpha )}u(t_k) \nonumber \\&\quad + \tau ^{-\alpha }\sum _{k=1}^m W_k^{(n,\alpha )}u(t_k), \end{aligned}$$

(4)

where $\omega _{0}^{(\alpha )}=\frac{2+\alpha }{2}g_{0}^{(\alpha )}$ $\omega _{k}^{(\alpha )}=\frac{2+\alpha }{2}g_{k}^{(\alpha )}-\frac{\alpha }{2}g_{k-1}^{(\alpha )}$, $g_{k}^{(\alpha )}=(-1)^{k}{\left( {\begin{array}{c}\alpha \\ k\end{array}}\right) }$, and the starting weights $W_{k}^{(n,\alpha )}$ are chosen such that (4) is exact for $v(t)=t^{\sigma _{m}}$ ($m=0,1,2,\ldots$), which leads to the system:

$$\begin{aligned} \sum _{k=1}^{m}W_k^{(n,\alpha )}k^{\sigma _m}=\frac{\Gamma (\sigma _m+1)}{\Gamma (\sigma _m+1-\alpha )}n^{\sigma _m-\alpha } -\sum _{k=0}^{n}\omega _{n-k}^{(\alpha )}k^{\sigma _m}. \end{aligned}$$

(5)

We apply the formula to the Caputo-tempered operator at $t=t_{n}$ to obtain

$$\begin{aligned} {_{0}^{C}D_{t}^{(\alpha ,\rho )}}u(t_{n})&=e^{-\rho t_{n}}{_{0}^{R}D_{t}^{\alpha }}\left( F(t_{n})-F(t_{0})\right) \\&=e^{-\rho t_{n}}\tau ^{-\alpha }\sum _{k=0}^{n}\omega _{n-k}^{(\alpha )}\left( F(t_{k})-F(t_{0})\right) \\&\quad +e^{-\rho t_{n}}\tau ^{-\alpha }\sum _{k=1}^{m}W_{k}^{(n,\alpha )}\\&\quad \times \left( F(t_{k})-F(t_{0})\right) +R_{2}^{n}\\&=\tau ^{-\alpha }\sum _{k=0}^{n}\omega _{n-k}^{(\alpha )}e^{-\rho t_{n-k}}u(t_{k})\\&\quad -e^{-\rho t_{n}}\tau ^{-\alpha }\sum _{k=0}^{n}\omega _{n-k}^{(\alpha )}u(t_{0}) \\&\quad +\tau ^{-\alpha }\sum _{k=1}^{m}W_{k}^{(n,\alpha )}e^{-\rho t_{n-k}}u(t_{k})\\&-e^{-\rho t_{n}}\tau ^{-\alpha }\sum _{k=1}^{m}W_{k}^{(n,\alpha )}u(t_{0})+R_{2}^{n} \\&:={_{0}^{C}{\mathfrak {D}}_{t}^{(\alpha ,\rho )}}u(t_{n})+R_{2}^{n}, \end{aligned}$$

where $F(t)=e^{\rho t}u(t)$ and $|R_{2}^{n}|\le C\tau ^{2}t_{n}^{\sigma _{m+1}-2-\alpha }$. To guarantee the WSGL formula has a global second-order convergence, we need $\sigma _{m+1}\ge 2+\alpha$. For more details about choosing the correction terms, one can refer to [24]. For our problem, at $t=t_{n}$, we have

$$\begin{aligned} {_{0}^{C}{\mathfrak {D}}_{t}^{(\alpha ,\rho )}}u^{n}=-k_{0}u^{n}. \end{aligned}$$

To verify the WSGL formula, we also give an example using the same parameters $k_{0}=2$, $\rho =0.5$, $u_{0}=1$ and $T=1$, to which the numerical results are given in Table 2. We can observe that without the correction terms ($m=0$), the WSGL formula only exhibits convergence rate $O(\tau ^{\alpha })$. While adding the correction terms into the WSGL formula, the second-order convergence is obtained for the case $\alpha =0.8$ with $m=2$. Although the optimal second-order convergence is not reached for the case $\alpha =0.4$ with $m=4$, the accuracy has been improved significantly. It reveals that the WSGL formula is also feasible to tackle the time-fractional tempered problem.

Table 2 The error and convergence order of the WSGL scheme for the tempered ODE (1) at $t=1$, where the parameters are $k_0=2$, $\rho =0.5$, $u_0=1$

Table 3 The error and convergence order of the L1 scheme for the example (6) at $t=1$, where the parameters are $\rho =0.5$, $u_0=1$, $r=\frac{2(2-\alpha )}{\alpha }$

2.4 Comparison of the two methods

In this section, we compare the L1 scheme on graded mesh with the WSGL scheme. To avoid the error caused by discretising the Mittag–Leffler function, we consider the following example:

$$\begin{aligned} {_0^CD_t^{(\alpha ,\rho )}}u(t)=f(t),~ 0<t\le T, \quad u(0)=u_0, \end{aligned}$$

(6)

with an exact solution $u(t)=u_{0}e^{-\rho t}\sum \limits _{k=0}^{8}t^{k\alpha }$. In the calculation, we choose the parameters $\rho =0.5$, $u_{0}=1$ and $T=1$. The related numerical results are shown in Tables 3 and 4, respectively. From Table 3, we can see that the optimal convergence order $2-\alpha$ is obtained using mesh grading $r=\frac{2(2-\alpha )}{\alpha }$. An important property is that the CPU time for the cases $\alpha =0.4$ and $\alpha =0.8$ is almost the same, which means the fractional index $\alpha$ does not affect the CPU time too much for a fixed N. This will facilitate the implementation of the fast calculation of the tempered operator in the subsequent section. Compared to the L1 scheme, the WSGL scheme is of high order, which shows better accuracy and convergence. However, when $\alpha$ is small, we need to add more correction terms, which increases the CPU time. In addition, the system (5) leads to an exponential Vandermonde type matrix. When the number of the correction terms is large, the matrix will be ill-conditioned, which may cause big round-off errors [24]. We can conclude that both methods have merits and drawbacks.

Table 4 The error and convergence order of the WSGL scheme for the example (6) at $t=1$, where the parameters are $\rho =0.5$, $u_0=1$

2.5 Fast calculation for the time-tempered Caputo fractional derivative

In this part, based on the L1 scheme, we consider a fast evaluation for the time-tempered Caputo fractional derivative, in which the derivative is to be split into two different parts: the local part and the history part. A fast approximation will be applied to the history part. We now give a detailed implementation. At $t=t_{n}$, we have

$$\begin{aligned} {^C_0D_t^{(\alpha ,\rho )}}u(t)\Big |_{t=t_n}&=e^{-\rho t}\,{^C_0D^{\alpha }_t}(e^{\rho t}u(t))\Big |_{t=t_n}\\&=\frac{e^{-\rho t_n}}{\Gamma (1-\alpha )} \int _0^{t_n}\frac{\left( e^{\rho s}u(s)\right) ' }{(t_n-s)^\alpha } ds\\&=\frac{e^{-\rho t_n}}{\Gamma (1-\alpha )} \int _0^{t_{n-1}}\frac{\left( e^{\rho s}u(s)\right) ' }{(t_n-s)^\alpha } ds\\&\quad +\frac{e^{-\rho t_n}}{\Gamma (1-\alpha )} \int _{t_{n-1}}^{t_n}\frac{\left( e^{\rho s}u(s)\right) ' }{(t_n-s)^\alpha } ds\\&:=C_h(t_n)+C_l(t_n), \end{aligned}$$

where $C_{h}(t_{n})$ is the history part and $C_{l}(t_{n})$ is the local part. For the local part $C_{l}(t_{n})$, we use the backward difference scheme to approximate the first-order derivative, which gives

$$\begin{aligned} C_l(t_n)&\approx \frac{e^{-\rho t_n}}{\Gamma (1-\alpha )}\frac{e^{\rho t_n}u(t_n)-e^{\rho t_{n-1}}u(t_{n-1})}{\tau _n}\\&\quad \times \int _{t_{n-1}}^{t_n}\frac{1 }{(t_n-s)^\alpha } ds\\&=\frac{u(t_n)-e^{-\rho \tau _{n}}u(t_{n-1})}{\tau _n^\alpha \Gamma (2-\alpha )}. \end{aligned}$$

With regard to the history part $C_{h}(t_{n})$, using integration by parts, we obtain

$$\begin{aligned} C_h(t_n)&=\frac{e^{-\rho t_n}}{\Gamma (1-\alpha )} \int _0^{t_{n-1}}\frac{\left( e^{\rho s}u(s)\right) ' }{(t_n-s)^\alpha } ds \nonumber \\&=\frac{e^{-\rho t_n}}{\Gamma (1-\alpha )} \left[ \frac{e^{\rho t_{n-1}}u(t_{n-1})}{\tau _n^\alpha }-\frac{u(t_{0})}{t_n^\alpha }\right. \nonumber \\&\quad \left. -\alpha \int _0^{t_{n-1}}\frac{e^{\rho s}u(s) }{(t_n-s)^{1+\alpha }} ds\right] . \end{aligned}$$

(7)

Next, we approximate the kernel $t^{-\beta }$ $(0<\beta <2)$ via a sum-of-exponentials (SOE) approximation [26]. For $t\in [\sigma ,T]$, $\sigma =\min \limits _{1\le n\le N}\{\tau _n \}$, there exist some positive real numbers $s_i$ and $\omega _i$ $(i=1,2,\ldots ,N_\mathrm{{exp}})$ such that

$$\begin{aligned} \left| t^{-\beta }- \sum _{i=1}^{N_\mathrm{{exp}}}\omega _ie^{-s_{i}t}\right| \le \varepsilon ,\quad t\in [\sigma ,T],\quad \beta \in (0,2), \end{aligned}$$

where

$$\begin{aligned}&N_\mathrm{{exp}}={\mathcal {O}}\left( \log \frac{1}{\varepsilon }\left( \log \log \frac{1}{\varepsilon }+\log \frac{T}{\sigma }\right) \right. \\&\left. +\log \frac{1}{\sigma }\left( \log \log \frac{1}{\varepsilon }+\log \frac{1}{\sigma }\right) \right) . \end{aligned}$$

Replacing the kernel $t^{-1-\alpha }$ in the history part (7) with the SOE approximation, we obtain

$$\begin{aligned}&\alpha \int _0^{t_{n-1}}\frac{e^{\rho s}u(s) }{(t_n-s)^{1+\alpha }} ds=\alpha \\&\qquad \times \int _0^{t_{n-1}} \sum _{i=1}^{N_\mathrm{{exp}}}\omega _ie^{-s_{i}(t_n-s)}e^{\rho s}u(s) ds\\&\quad =\alpha \sum _{i=1}^{N_\mathrm{{exp}}}\omega _i e^{\rho t_n}{\mathcal {U}}_{his,i}(t_n), \end{aligned}$$

where ${\mathcal {U}}_{his,i}(t_n)=\int _0^{t_{n-1}}e^{-(\rho +s_i)(t_n-s)}u(s) ds$. For the term ${\mathcal {U}}_{his,i}(t_n)$, we have the following recurrence relation

$$\begin{aligned} {\mathcal {U}}_{his,i}(t_n)&=e^{-(\rho +s_i)\tau _n}{\mathcal {U}}_{his,i}(t_{n-1})\\&\quad +\int _{t_{n-2}}^{t_{n-1}}e^{-(\rho +s_i)(t_n-s)}u(s) ds. \end{aligned}$$

Using a linear interpolation function to approximate u(t) on the interval $[t_{n-2},t_{n-1}]$, we can derive

$$\begin{aligned}&\int _{t_{n-2}}^{t_{n-1}}e^{-(\rho +s_i)(t_n-s)}u(s) ds\\&\quad \approx \frac{e^{-(\rho +s_i)\tau _n}}{(\rho +s_i)^2\tau _{n-1}} \bigg [ (e^{-(\rho +s_i)\tau _{n-1}}-1\\&\qquad +(\rho +s_i)\tau _{n-1})u(t_{n-1})\\&\qquad +(1-e^{-(\rho +s_i)\tau _{n-1}}\\&\qquad -e^{-(\rho +s_i)\tau _{n-1}}(\rho +s_i)\tau _{n-1})u(t_{n-2}) \bigg ]. \end{aligned}$$

Finally, the fast approximation of ${^C_0D_t^{(\alpha ,\rho )}}u(t_n)$, $n\ge 2$ can be written as

$$\begin{aligned}&{^{FC}_0{\mathbb {D}}_t^{(\alpha ,\rho )}}u(t_n)\nonumber \\&\quad :=\frac{u(t_n)-e^{-\rho \tau _{n}}u(t_{n-1})}{\tau _n^\alpha \Gamma (2-\alpha )}\nonumber \\&\quad \quad +\frac{e^{-\rho t_n}}{\Gamma (1-\alpha )} \left[ \frac{e^{\rho t_{n-1}}u(t_{n-1})}{\tau _n^\alpha }\right. \nonumber \\&\quad \quad \left. -\frac{u(t_{0})}{t_n^\alpha }-\alpha \sum _{i=1}^{N_\mathrm{{exp}}}\omega _i e^{\rho t_n}{\mathcal {U}}_{his,i}(t_n)\right] \nonumber \\&\quad =\frac{u(t_n)}{\tau _n^\alpha \Gamma (2-\alpha )}-\frac{\alpha e^{-\rho \tau _{n}}u(t_{n-1})}{\tau _n^\alpha \Gamma (2-\alpha )} \nonumber \\&\qquad -\frac{e^{-\rho t_{n}}u(t_{0})}{\Gamma (1-\alpha )t_n^\alpha }-\frac{\alpha }{\Gamma (1-\alpha )}\sum _{i=1}^{N_\mathrm{{exp}}}\omega _i {\mathcal {U}}_{his,i}(t_n). \end{aligned}$$

(8)

When $n=1$, we have

$$\begin{aligned} {^{FC}_0{\mathbb {D}}_t^{(\alpha ,\rho )}}u(t_{1}):=\frac{u(t_{1})-e^{-\rho \tau _{1}}u(t_{0})}{\tau _1^\alpha \Gamma (2-\alpha )}. \end{aligned}$$

(9)

For the truncation error of the fast algorithm, similar to the proof in [18, 26], it is straightforward to conclude the following lemma.

Lemma 2

Suppose $|\frac{du}{dt}|\le C(1+t^{\alpha -l})$ for $l=0,1,2$ and $0<\alpha <1$ and let $\varepsilon$ be the desired precision, then we have the truncation error of the fast schemes (8)–(9)

$$\begin{aligned}&|R_3^n|=\left| {_0^CD_t^{(\alpha ,\rho )}}u(t_n)- {^{FC}_0{\mathbb {D}}_t^{(\alpha ,\rho )}}u(t_n)\right| \\&\le C \left( n^{-\min \{2-\alpha ,r\alpha \}}+\varepsilon \right) . \end{aligned}$$

In the subsequent sections, we adapt the tempered operator into different models to observe the behaviour of the tempered models.

3 Application I: tempered Bloch equations

Ordinary matter all has nuclei with induced magnetic moments, leading to a nuclear paramagnetic polarisation in a constant magnetic field when an equilibrium is established. Imposing an external radiofrequency field to the constant field at a right angle can cause the Larmor precession of the moments around the constant field, which can induce an electromotive force in a coil [27, 28]. Furthermore, the induced electromotive force can be transferred into visible signals. This is the principle of nuclear magnetic resonance (NMR) and magnetic resonance imaging (MRI). NMR has been widely used to analyse complex biological materials in chemistry, medicine, and engineering. The phenomenological Bloch equation is utilised to capture the magnetisation dynamics. The empirical vector form of the Bloch equation is:

$$\begin{aligned} \frac{d\mathbf{M}(t)}{dt}&=\gamma \mathbf{M}(t)\times \mathbf{B}(t) +\frac{M_0-M_z(t)}{T_1}{} \mathbf{k}\\&\quad -\frac{M_x(t)\mathbf{i}+M_y(t)\mathbf{j}}{T_2}, \end{aligned}$$

where ${\mathbf {M}}=[M_{x},M_{y},M_{z}]^{T}$ is the magnetisation, ${\mathbf {B}}$ is the static magnet field, ${\mathbf {B}}(t)=B_{0}{\mathbf {k}}$, $M_{0}$ is the equilibrium magnetisation, $\gamma$ is the gyromagnetic ratio with a relation to the Larmor frequency $\omega _{0}=\gamma B_{0}$. The relaxation terms refer to the time return to equilibrium, where $T_{1}$ is the ‘longitudinal’ relaxation time and $T_{2}$ is the ‘transversal’ relaxation time. To investigate heterogeneous, porous and complex materials exhibiting memory, Magin et al. [29, 30] generalise the Bloch equation to the time-fractional Bloch equation. From the view of CTRW, the time-fractional operator relates to an infinite waiting time. Since the waiting time for a water proton cannot be arbitrarily long, we modify the equation as the tempered Bloch equation with finite waiting time, which is written in the form

$$\begin{aligned}&\mu _1^{\alpha -1}{^{C}_0D_t^{(\alpha ,\rho )}}M_z(t)=\frac{M_0-M_z(t)}{T_1},\\&\mu _2^{\alpha -1}{^{C}_0D_t^{(\alpha ,\rho )}}M_x(t)=\omega _0M_y(t)-\frac{M_x(t)}{T_2},\\&\mu _2^{\alpha -1}{^{C}_0D_t^{(\alpha ,\rho )}}M_y(t)=-\omega _0M_x(t)-\frac{M_y(t)}{T_2}. \end{aligned}$$

These equations can be recast into

$$\begin{aligned} {^{C}_0D_t^{(\alpha ,\rho )}}M_z(t)&=\frac{M_0-M_z(t)}{T'_1}, \end{aligned}$$

(10)

$$\begin{aligned} {^{C}_0D_t^{(\alpha ,\rho )}}M_x(t)&=\varpi _0M_y(t)-\frac{M_x(t)}{T'_2}, \end{aligned}$$

(11)

$$\begin{aligned} {^{C}_0D_t^{(\alpha ,\rho )}}M_y(t)&=-\varpi _0M_x(t)-\frac{M_y(t)}{T'_2}, \end{aligned}$$

(12)

where $\mu _{1}$ and $\mu _{2}$ are constants to preserve units, $\alpha$ is the fractional index, $\rho$ is the tempered parameter, $\varpi _{0}=\omega _{0}\mu _{2}^{1-\alpha }$, $T_{1}^{\prime }=\mu _{1}^{\alpha -1}T_{1}$ and $T_{2}^{\prime }=\mu _{2}^{\alpha -1}T_{2}$. When $\rho =0$, the system recovers the time-fractional Bloch equation proposed in [29]. We consider the analytical and numerical solutions for the system.

3.1 Analytical solution

Applying the Laplace transform to (10) leads to

$$\begin{aligned} \overline{M_z}(s)=\frac{(s+\rho )^{\alpha -1}}{(s+\rho )^{\alpha }+k_1}M_z(0) +\frac{k_2}{s[(s+\rho )^{\alpha }+k_1]}, \end{aligned}$$

where $k_1=\frac{1}{T'_1}$ and $k_2=\frac{M_0}{T'_1}$. Using the property ${\mathcal {L}}\{e^{\rho t}f(t)\}={\bar{f}}(s-\rho )$, ${\mathcal {L}}\{t^{\alpha -1}E_{\alpha ,\alpha }(-kt^\alpha )\}=\frac{1}{s^\alpha +k}$ and ${\mathcal {L}}^{-1}\{{\bar{f}}(s){\bar{g}}(s)\}=\int _0^tf(\eta )g(t-\eta )d\eta$ and applying the inverse Laplace transform, we can derive

$$\begin{aligned} M_z(t)&=M_z(0)e^{-\rho t}E_{\alpha }(-k_1t^{\alpha })\nonumber \\&\quad \;+k_2\int _0^te^{-\rho \eta }\eta ^{\alpha -1}E_{\alpha ,\alpha }(-k_1\eta ^{\alpha })d\eta \nonumber \\&=M_z(0)e^{-\rho t}E_{\alpha }\left( -\frac{t^{\alpha }}{T'_1}\right) \nonumber \\&\quad \;+\frac{M_0}{T'_1}\int _0^te^{-\rho \eta }\eta ^{\alpha -1}E_{\alpha ,\alpha }\left( -\frac{\eta ^{\alpha }}{T'_1}\right) d\eta . \end{aligned}$$

(13)

Suppose that $M_{+}(t)=M_x(t)+iM_y(t)$ with $M_{+}(0)=M_x(0)+iM_y(0)$. Combining (11) and (12) gives

$$\begin{aligned} {^{C}_0D_t^{(\alpha ,\rho )}}M_+(t)&=-i\varpi _0M_+(t)-\frac{M_+(t)}{T'_2}. \end{aligned}$$

(14)

Applying the Laplace transform and the inverse Laplace transform to (14), we can obtain

$$\begin{aligned} M_+(t)=M_{+}(0)e^{-\rho t}E_{\alpha }(-k_3t^{\alpha }), \end{aligned}$$

(15)

where $k_3=i\varpi _0+\frac{1}{T'_2}$. As the analytical solutions (13) and (15) involve the Mittag–Leffler function and its integral, which is very challenging to calculate, we need to resort to a numerical solution of the problem.

3.2 Numerical solution

We use the L1 scheme on graded mesh to solve the tempered problem (10)–(12). Define $\mathbf{M}^n$ as the numerical approximation of $\mathbf{M}(t_n)$. At $t=t_n$, we can obtain the following numerical scheme:

$$\begin{aligned} {_0^C{\mathbb {D}}_t^{(\alpha ,\rho )}}M_z^n&=\frac{M_0-M_z^n}{T'_1},\\ {_0^C{\mathbb {D}}_t^{(\alpha ,\rho )}}M_x^n&=\varpi _0M_y^n-\frac{M_x^n}{T'_2},\\ {_0^C{\mathbb {D}}_t^{(\alpha ,\rho )}}M_y^n&=-\varpi _0M_x^n-\frac{M_y^n}{T'_2}. \end{aligned}$$

3.3 Numerical examples

Some numerical results are shown in Figs. 1, 2, 3. Figure 1a shows that the fractional order $\alpha$ enhances the relaxation at a small time scale and then delays the relaxation later compared to the classical case $\alpha =1$. The smaller the fractional index is, the slower it converges to its final asymptotic value. Figure 1b illustrates that the tempered parameter $\rho$ can accelerate the process to reach its asymptotic value. Contrary to the fractional order $\alpha$, the larger $\rho$ is, the more rapid it converges to its maximum value. In addition, with a large $\alpha$ and moderate $\rho$ ($\alpha =0.9$, $\rho =0.2$), we can obtain a similar final state by choosing small $\alpha$ only ($\alpha =0.5$, $\rho =0.0$). This facilitates avoiding choosing too small $\alpha$ since a small fractional index means strong heterogeneity or low regularity of the solution.

Fig. 1

The evolution of $M_z(t)$ with different $\alpha$ (a) and different $\rho$ (b), where the parameters used are $M_z(0)=0$, $M_0=100$, $T_1'=1(ms)^{\alpha }$

Fig. 2

The evolution of $M_y(t)$ with different $\alpha$ and $\rho$, where the parameters used are $T_2'=20(ms)^{\alpha }$ $\tilde{f}_0=\frac{\varpi _0}{2\pi }=160 Hz$, $M_x(0)=0$, $M_y(0)=100$

The effects of $\alpha$ and $\rho$ on $M_{y}(t)$ are presented in Fig. 2. We can see that, for the classical case $\alpha =1$, $\rho =0$, it needs a long relaxation time to decay to zero. The fractional order boosts the decay. The smaller $\alpha$ is, the faster it decays. The tempered parameter $\rho$ promotes the decay further. We can observe that it decays faster with $\alpha =0.9$, $\rho =0.2$ than with $\alpha =0.8$, $\rho =0.0$, which is similar to the discussion above. Figure 3 displays the effects of $\rho$ on the evolution of $M_{x}(t)$ versus $M_{y}(t)$. Similarly, the tempered parameter advances the decay process to zero.

Fig. 3

The evolution of $M_x(t)$ versus $M_y(t)$ with different $\rho$, where the parameters used are $T_2'=20(ms)^{\alpha }$ $\tilde{f}_0=\frac{\varpi _0}{2\pi }=160Hz$, $M_x(0)=0$, $M_y(0)=100$

Fig. 4

Voxel-level data fittings of the classical model (Monoexp), time-fractional model (Fractional) and tempered time-fractional model (Tempered). Here we consider two different voxels (locations (a) (132, 55, 47) and (b) (80, 60, 68)) in different brain regions, of which the data are extracted from [32]. The mean square error (MSE) is defined as $\sum _{j=1}^n(y'_j-y_j)^2/n$, where $y_j$ is the real data, $y'_j$ is the predicted data and n is the number of data

The Bloch equation is related to the MRI signal by $S(t)=A_0\sqrt{M_x^2(t)+M_y^2(t)}+C$, where $A_0$ is the amplitude of the signal and C is a constant accounting for the background noise [31]. Figure 4 shows the voxel-level temporal fitting based on three models, namely the classical mono-exponential model ($A_0\exp (-t/T_2)+C$), the time-fractional model ($A_0E_{\alpha }(-t^{\alpha }/T_2)+C$) [29] and our tempered time-fractional model. The experimental data are extracted from the paper [32]. The ‘lsqcurvefit’ MATLAB routine is used to perform all the data fittings with a maximum number of iteration and a relative tolerance set as $1\times 10^7$ and $1\times 10^{-7}$, respectively. The specific fitting procedure can be referred to [31]. This procedure is generally followed here to determine the initial conditionsfor our model with $A_0$, $T_2$ and C based on the results from the exponential model and $\alpha$ from the fractional model. These parameters are constrained in ranges to $70{-}130\%$ of their starting values. $\omega$ and $\rho$ are new parameter to the other two models, which are initialised to 10 and 0.1 with ranges confined in $0{-}250$ and $0{-}1$. The mean-square error (MSE) is adopted to measure the quality of the fitting. From Fig. 4, we can see that for the signal at position (a), both the classical model and the time-fractional model have a large MSE, while the tempered model has a clear improvement in the accuracy of fitting the voxel-level MRI data with less than a half MSE. For the data at position (b), the classical model does not fit well with an MSE 558.10. Although the time-fractional model improves the fitting little, the MSE is still large with a value of 524.59. In comparison, the MSE for the tempered time-fractional model reduces significantly to 132.11. We conclude that the tempered time-fractional model is effective in fitting the MRI signal.

4 Application II: tempered diffusion problems

In the CTRW model, when a Poissonian waiting time pdf together with a Gaussian jump length pdf is applied, the standard diffusion equation can be derived, describing the Brownian motion. When the finite waiting time pdf is replaced with a divergent long-tailed waiting time pdf, it leads to the time-fractional diffusion equation used to characterise the anomalous diffusion [2]. Different from the standard diffusion with mean-squared displacement $\langle x^{2}(t)\rangle \sim t$, the mean-squared displacement of time-fractional diffusion has the property $\langle x^{2}(t)\rangle \sim t^{\alpha }$, where $0<\alpha <1$ is the fractional order. In this section, we consider a tempered long-tailed waiting time pdf to solve a tempered diffusion problem

$$\begin{aligned}&{^{C}_0D_t^{(\alpha ,\rho )}}u(x,t)=D\frac{\partial ^2u(x,t)}{\partial x^2}+f(x,t),\nonumber \\&\quad x\in (0,l),~ t\in (0,T], \end{aligned}$$

(16)

subject to the initial and boundary conditions

$$\begin{aligned}&u(x,0)=\psi (x),\quad x\in [0,l],\quad u(0,t)=u(l,t)=0,\nonumber \\&\quad t\in (0,T]. \end{aligned}$$

(17)

4.1 The mean-squared displacement

To study the mean-squared displacement of the problem, we suppose $f(x,t)=0$, $\psi (x)=\delta (x)$. Applying the Laplace transform and Fourier transform successively to obtain

$$\begin{aligned} \widehat{{\overline{u}}}(\omega ,s)=\frac{(s+\rho )^{\alpha -1}}{(s+\rho )^{\alpha }+D\omega ^2}. \end{aligned}$$

As $-\frac{\partial ^2 \widehat{{\overline{u}}}(\omega ,s)}{\partial \omega ^2}\big |_{\omega =0}=-\frac{2D}{(s+\rho )^{\alpha +1}}$, we derive

$$\begin{aligned} \langle x^2(t) \rangle =e^{-\rho t}\frac{2Dt^\alpha }{\Gamma (1+\alpha )}. \end{aligned}$$

Different to the mean-squared displacement of the time-fractional diffusion in [2], one tempered factor $e^{-\rho t}$ is added.

4.2 Analytical solution

We first use the finite Fourier transform technique to solve the problem (16)–(17). The associated eigenvalue problem that needs to be solved is

$$\begin{aligned} -\frac{d^2 \varphi }{dx^2}=\lambda ^2\varphi ,\quad \varphi (0)=0,\quad \varphi (l)=0, \end{aligned}$$

to which the solution is $\lambda _n^2=\frac{n^2\pi ^2}{l^2}$ for $n=1,2,\ldots$ and $\varphi _n(x)=\sqrt{\frac{2}{l}}\sin \left( \lambda _nx\right)$. Define the finite Fourier transform ${\widetilde{u}}(\lambda _n,t):=\langle u, \varphi _n\rangle =\int _0^lu(x,t)\varphi _n(x)dx$. Applying it to (16), we obtain

$$\begin{aligned} \left\langle {^{C}_0D_t^{(\alpha ,\rho )}}u(x,t), \varphi _n\right\rangle =D\left\langle \frac{\partial ^2u(x,t)}{\partial x^2} ,\varphi _n\right\rangle +\langle f(x,t),\varphi _n\rangle . \end{aligned}$$

(18)

As $D\left\langle \frac{\partial ^2u(x,t)}{\partial x^2} ,\varphi _n\right\rangle =-D\left\langle \frac{d^2 \varphi _n}{d x^2} ,u\right\rangle =-D\lambda _n^2\langle u, \varphi _n\rangle$, (18) can be recast into

$$\begin{aligned} {^{C}_0D_t^{(\alpha ,\rho )}}{\widetilde{u}}=-D\lambda _n^2{\widetilde{u}}+{\widetilde{f}}. \end{aligned}$$

(19)

Applying the Laplace transform to (19) and denoting $\overline{{\widetilde{u}}}(\lambda _n,s)={\mathcal {L}}\{{\widetilde{u}}(\lambda _n,t)\}$, we can derive

$$\begin{aligned} \overline{{\widetilde{u}}}=\frac{(s+\rho )^{\alpha -1}{\widetilde{u}}(0)}{(s+\rho )^{\alpha }+D\lambda _n^2} +\frac{\overline{{\widetilde{f}}}}{(s+\rho )^{\alpha }+D\lambda _n^2}. \end{aligned}$$

(20)

Imposing the inverse Laplace transform to (20) leads to

$$\begin{aligned}&{\widetilde{u}}(\rho _n,t)={\widetilde{u}}(0)e^{-\rho t}E_{\alpha }(-D\lambda _n^2 t^\alpha )\\&\quad +\int _0^te^{-\rho \eta }\eta ^{\alpha -1}E_{\alpha ,\alpha }\left( -D\lambda _n^2 \eta ^\alpha \right) {\widetilde{f}}(t-\eta )d\eta . \end{aligned}$$

Then the analytical solution can be derived as

$$\begin{aligned} u(x,t)&=\sqrt{\frac{2}{l}}\sum _{n=1}^{\infty }\left[ {\widetilde{u}}(0)e^{-\rho t}E_{\alpha }\left( -D\lambda _n^2 t^\alpha \right) \right. \nonumber \\&\quad \left. +\int _0^te^{-\rho \eta }\eta ^{\alpha -1}E_{\alpha ,\alpha }\right. \nonumber \\&\quad \left. \left( -D\lambda _n^2 \eta ^\alpha \right) {\widetilde{f}}(t-\eta )d\eta \right] \sin \left( \frac{n\pi x}{l}\right) , \end{aligned}$$

(21)

where ${\widetilde{u}}(0)=\langle \psi , \varphi _n\rangle$, ${\widetilde{f}}=\langle f, \varphi _n\rangle$.

4.3 Numerical solution

4.3.1 Numerical scheme

Now we consider the numerical solution to (16)–(17). Firstly, we denote a uniform mesh in space domain. Define $h=\frac{l}{M}$, $x_i=ih$, $i=0,1,2,\ldots ,M$, where M is a positive integer. For the space Laplacian operator, we use the standard second-order central difference scheme to approximate:

$$\begin{aligned} \frac{\partial ^2u(x_i,t_n)}{\partial x^2}&=\frac{u(x_{i-1},t_n)-2u(x_i,t_n)+u(x_{i+1},t_n)}{h^2}\\&\quad +{\mathcal {O}}(h^2)\\&:=\delta _x^2u(x_i,t_n)+{\mathcal {O}}(h^2). \end{aligned}$$

Let $u_i^n$ be the numerical solution to $u(x_i,t_n)$. Then, the numerical scheme to (16)–(17) is

$$\begin{aligned}&{_0^C{\mathbb {D}}_t^{(\alpha ,\rho )}}u_i^n=D\delta _x^2u_i^n+f_i^n,\nonumber \\&\quad 1\le i\le M-1,~1\le n\le N, \end{aligned}$$

(22)

$$\begin{aligned}&u_0^n=u_M^n=0,\quad \mathrm{{for}}~0<n\le N, \end{aligned}$$

(23)

$$\begin{aligned}&u_i^0=\psi (x_i), \quad \mathrm{{for}}~0\le i\le M. \end{aligned}$$

(24)

And the fast numerical scheme to (16)–(17) is

$$\begin{aligned}&{^{FC}_0{\mathbb {D}}_t^{(\alpha ,\rho )}}u_i^n=D\delta _x^2u_i^n+f_i^n,\nonumber \\&\quad 1\le i\le M-1,~1\le n\le N, \end{aligned}$$

(25)

$$\begin{aligned}&u_0^n=u_M^n=0,\quad \mathrm{{for}}~0<n\le N, \end{aligned}$$

(26)

$$\begin{aligned}&u_i^0=\psi (x_i), \quad \mathrm{{for}}~0\le i\le M. \end{aligned}$$

(27)

4.3.2 Stability

Theorem 1

The numerical scheme (22) is unconditional stable, and it holds that

$$\begin{aligned} ||u^n||_{\infty }&\le ||u^0||_{\infty }+\tau _n^\alpha \Gamma (2-\alpha )\sum _{j=1}^n\theta _{n,j}||f^j||_{\infty }, \end{aligned}$$

where $\theta _{n,n}=1$, $\theta _{n,j}=\sum _{k=1}^{n-j}\tau _{n-k}^{\alpha }(d_{n,k}-d_{n,k+1})\theta _{n-k,j}$.

Proof

According to Lemma 4.2 in [22], it is straightforward to derive

$$\begin{aligned}&||e^{\rho t_n} u^n||_{\infty } \le ||u^0||_{\infty }+\tau _n^\alpha \Gamma (2-\alpha )\\&\sum _{j=1}^n\theta _{n,j}||e^{\rho t_j}f^j||_{\infty }, \end{aligned}$$

which leads to

$$\begin{aligned} ||u^n||_{\infty }&\le e^{-\rho t_n}||u^0||_{\infty }+e^{-\rho t_n}\tau _n^\alpha \Gamma (2-\alpha )\\&\sum _{j=1}^n\theta _{n,j}||e^{\rho t_j}f^j||_{\infty }\\&\le ||u^0||_{\infty }+\tau _n^\alpha \Gamma (2-\alpha )\sum _{j=1}^n\theta _{n,j}||f^j||_{\infty }. \end{aligned}$$

The proof is completed. $\square$

4.3.3 Convergence

Theorem 2

Define $u^n$ and $U^n$ as the exact and numerical solution vector, respectively. Then, there exists a positive constant C independent of $\tau$ and h such that

$$\begin{aligned} ||u^n-U^n||_{\infty }\le C(h^2+T^\alpha N^{-\min \{2-\alpha ,r\alpha \}}). \end{aligned}$$

Proof

The truncation error of (22) at $(x_i,t_n)$ is

$$\begin{aligned} |{\mathfrak {R}}^n|\le C(h^2+T^\alpha n^{-\min \{2-\alpha ,r\alpha \}}). \end{aligned}$$

Invoking Theorem 1, we have

$$\begin{aligned} \max |u(x_i,t_n)-u^n_i|&\le C \tau _n^\alpha \Gamma (2-\alpha )\sum _{j=1}^n\theta _{n,j}||{\mathfrak {R}}^j||_{\infty }\\&\le C\tau _n^\alpha \Gamma (2-\alpha )\sum _{j=1}^n\theta _{n,j} (h^2\\&\quad +T^\alpha j^{-\min \{2-\alpha ,r\alpha \}})\\&\le CT^\alpha (h^2+ N^{-\min \{2-\alpha ,r\alpha \}}), \end{aligned}$$

where the following inequality [22] has been used

$$\begin{aligned} \tau _n^{\alpha }\sum _{j=1}^nj^{-\beta }\theta _{n,j}\le \frac{T^\alpha N^{-\beta }}{1-\alpha }. \end{aligned}$$

The proof is completed. $\square$

Similarly, we can obtain the convergence of the fast numerical scheme (25)–(27).

Theorem 3

Define $u^n$ and $U^n$ as the exact and numerical solution vector, respectively. Then, there exists a positive constant C independent of $\tau$ and h such that

$$\begin{aligned} ||u^n-U^n||_{\infty }\le C(h^2+T^\alpha N^{-\min \{2-\alpha ,r\alpha \}}+\varepsilon ). \end{aligned}$$

4.3.4 Numerical examples

In this section, for problem (16)–(17), we consider the case $f(x,t)=0$, $u(x,0)=\sin (x)$ in the domain $[0,\pi ]$. Using (21), it is straightforward to derive the analytical solution $u(x,t)=e^{-\rho t}E_{\alpha }(-Dt^{\alpha })\sin (x)$. We apply the numerical scheme (22)–(24) and the fast numerical scheme (25)–(27) to solve the problem. All of the computations were conducted using MATLAB R2018a on a DELL desktop with the configuration: Intel(R) Core(TM) i7-6700 CPU3.40GHz and RAM 16.0 GB. Firstly, a comparison between the numerical solution and the exact solution at different times is plotted in Fig. 5, where the parameters chosen are $M=40$, $N=M^{2}$, $D=0.5$, $\alpha =0.8$, $\rho =0.5$. It can be observed that the numerical solution agrees with the exact solution very well.

Fig. 5

The comparison between the numerical solution and the exact solution at different times, where the parameters chosen are $M=40$, $N=M^2$, $D=0.5$, $\alpha =0.8$, $\rho =0.5$

Table 5 The error and convergence order of the L1 scheme (22) and the fast numerical scheme (25) for different $\alpha$ and N at $t=1$, where the parameters used are $M=2^{11}$, $D=2$, $\rho =0.5$, $r=\frac{2(2-\alpha )}{\alpha }$, $\varepsilon =10^{-9}$

Table 6 The error and CPU time comparison between the L1 scheme (22) and the fast numerical scheme (25) for different $\alpha$ and M with $N=M^2$ at $t=1$, where the parameters utilised are $D=2$, $\rho =0.5$, $r=\frac{2(2-\alpha )}{\alpha }$, $\varepsilon =10^{-9}$

Next, we compare the convergence of the two numerical schemes. Table 5 shows the maximum error and convergence order of the two numerical schemes for different $\alpha$ and N at $t=1$, where the parameters used are $M=2^{11}$, $D=2$, $\rho =0.5$, $r=\frac{2(2-\alpha )}{\alpha }$, $\varepsilon =10^{-9}$. It can be seen that the optimal $2-\alpha$ order is obtained, which illustrates the effectiveness of the L1 scheme on the graded mesh to solve the tempered problem. Table 6 presents the error and CPU time comparison between the normal numerical scheme (22) and the fast numerical scheme (25) for different $\alpha$ and M with $N=M^{2}$ at $t=1$, where the parameters utilised are $D=2$, $\rho =0.5$, $r=\frac{2(2-\alpha )}{\alpha }$, $\varepsilon =10^{-9}$. We can observe that the L1 scheme on the graded mesh is time-consuming for a large M and the amount of time for $\alpha =0.4$ and $\alpha =0.8$ is not too much different. In contrast to the normal numerical scheme, the fast numerical scheme can reduce the CPU time significantly without losing accuracy, and there is a clear difference for the total time between the cases $\alpha =0.4$ and $\alpha =0.8$. This demonstrates the strong feasibility and applicability to adapt the fast method to deal with the tempered problem. Furthermore, we analyse the impacts of the fractional index $\alpha$ and tempered parameter $\rho$. Figure 6a shows that for the general anomalous diffusion $\rho =0$, the fractional order can boost the diffusion compared to the classical diffusion ($\alpha =1$). The smaller the fractional order is, the faster it diffuses. Compared to the fractional order $\alpha$, the tempered parameter $\rho$ (Fig. 6b) can accelerate the diffusion further. The larger $\rho$ is, the more rapid it decays, and with a large $\alpha$ and a moderate $\rho$ ($\alpha =0.9$, $\rho =3.0$) it diffuses faster than with a single fractional order $\alpha$ ($\alpha =0.6$), which affirms the discussion aforementioned again.

5 Application III: tempered two-layered problems

It is well known that the diffusion-induced MRI signal attenuation curve diverges from the mono-exponential decay at high b-values for human brain tissues. A stretched exponential model [33] was proposed to describe the diffusion-induced signal attenuation effectively, which proves to be a fundamental extension of the fractional Bloch–Torrey equation [34]. In [35], Zhou and co-workers applied the fractional models to analyse the diffusion images of human brain tissues in vivo, in which there was a clear contrast between the grey matter and white matter for the diffusivity (D) and fractional index ($\beta$). For example, the diffusivity and fractional index for the white matter are $D\approx (0.41\pm 0.008)\times 10^{-3}$ and $\beta \approx 0.64\pm 0.01$, while those for the grey matter are $D\approx (0.66\pm 0.007)\times 10^{-3}$ and $\beta \approx 0.82\pm 0.01$. It means different tissues exhibit distinct memory or heterogeneity in a heterogeneous medium, which motives us to generalise the fractional-order model to consider the anomalous diffusion in a composite material. To date, Zeng et al. [36] used a discrete least squares collocation method to deal with a coupled system of time-fractional diffusion equations with different fractional indices in an irregularly shaped region. Feng et al. [21, 37] proposed a two-dimensional time-fractional model to study the underlying transport phenomena in a binary medium based on the Riemann–Liouville fractional derivative, in which it showed that the generalised transport model could exhibit the correct physical solution behaviour and produce a more accurate overall mass balance. In this section, we focus on a two-layered problem with a tempered operator (see Fig. 7):

$$\begin{aligned} {^{C}_0D_t^{(\alpha _1,\rho _1)}}X_1(x,t)&=D_1\frac{\partial ^2X_1(x,t)}{\partial x^2}+S_{a_1}+S_{b_1}X_1,\nonumber \\&\quad x\in (l_0,l_1), \end{aligned}$$

(28)

$$\begin{aligned} {^{C}_0D_t^{(\alpha _2,\rho _2)}}X_2(x,t)&=D_2\frac{\partial ^2X_2(x,t)}{\partial x^2}+S_{a_2}+S_{b_2}X_2,\nonumber \\&\quad x\in (l_1,l_2), \end{aligned}$$

(29)

subject to the initial and boundary conditions

$$\begin{aligned}&X_1(x,0)=X_{1,0}(x),\quad X_2(x,0)=X_{2,0}(x), \end{aligned}$$

(30)

$$\begin{aligned}&X_1(l_0,t) =f_L(t),\quad X_2(l_2,t)=f_R(t), \end{aligned}$$

(31)

Fig. 6

The impacts of the fractional index $\alpha$ (Figure a) and tempered parameter $\rho$ (Figure b) on the diffusion profile at $t=0.1$ with $D=2.5$, $M=40$, $N=M^2$

where $D_i$ is the diffusivity coefficients, $S_{a_i}$ and $S_{b_i}$ are some constants. To guarantee the flux exchange at the interface $x=l_1$ is consistent with intrinsic physics, we define the following boundary conditions

$$\begin{aligned} X_1(l_1,t)&=X_2(l_1,t),\quad D_{1}\frac{\partial X_1(l_1,t) }{\partial x}\nonumber \\&=D_{2}\frac{\partial X_2(l_1,t) }{\partial x}. \end{aligned}$$

(32)

Fig. 7

An illustration of a two-layered problem

5.1 Semi-analytical solution

Problem (28)–(32) will be solved using the finite Fourier and Laplace transforms. Define $d_i=l_i-l_{i-1},i=1,2$, $\lambda _{i,n}$ be the eigenvalues and $\varphi _{i,n}(x)$ be the corresponding eigenfunctions. Then, the Sturm–Liouville system in each layer reads:

$$\begin{aligned}&-\frac{d^2 \varphi _{1,n}(x)}{dx^2} = \lambda _{1,n}^2 \varphi _{1,n}(x),\\&\quad \varphi _{1,n}(l_{0})=0,\quad \frac{\varphi _{1,n}(l_{1})}{dx}=0,\\&-\frac{d^2 \varphi _{2,n}(x)}{dx^2} = \lambda _{2,n}^2 \varphi _{2,n}(x),\\&\quad \frac{\varphi _{2,n}(l_{1})}{dx}=0,\quad \varphi _{2,n}(l_{2})=0. \end{aligned}$$

It is straightforward to derive that $\lambda _{1,n}=\frac{(2n+1)\pi }{2d_1}$, $\varphi _{1,n}=\sqrt{\frac{2}{d_1}}\sin [\lambda _{1,n}(x-l_{0})]$ and $\lambda _{2,n}=\frac{(2n+1)\pi }{2d_2}$, $\varphi _{2,n}=\sqrt{\frac{2}{d_2}}\sin [\lambda _{2,n}(l_{2}-x)]$. Define the finite Fourier transform ${\widetilde{X}}_i(\lambda _{i,n},t) := \langle X_i,\varphi _{i,n} \rangle = \int _{l_{i-1}}^{l_i} \, X_i(x,t)\, \varphi _{i,n}(x)\, dx$ within the $i^{\mathrm {th}}$ layer. We apply the finite Fourier transform to obtain

$$\begin{aligned} {^{C}_0D_t^{(\alpha _i,\rho _i)}} {\widetilde{X}}_i&=D_{i} \left\langle \frac{\partial ^2 X_i}{\partial x^2}, \varphi _{i,n} \right\rangle + S_{a_i} \left\langle 1,\varphi _{i,n} \right\rangle \nonumber \\&\quad + S_{b_i} \,\widetilde{{X}_i},\quad i=1,2. \end{aligned}$$

(33)

Integrating by parts the first term on the right-hand side yields

$$\begin{aligned} D_{i} \left\langle \frac{\partial ^2 X_i}{\partial x^2}, \varphi _{i,n} \right\rangle&= -D_{i} \left\{ \left\langle -\frac{d^2}{dx^2} \varphi _{i,n}, X_i \right\rangle \right. \nonumber \\&\quad \;+ \left[ X_i(l_{i},t)\frac{d \varphi _{i,n}}{dx}(l_{i})\right. \nonumber \\&\quad \; \left. -X_i(l_{i-1},t) \frac{d \varphi _{i,n}}{dx}(l_{i-1})\right. \nonumber \\&\quad \; - \left. \left. \frac{\partial X_i}{\partial x}(l_i,t) \varphi _{i,n}(l_i)\right. \right. \nonumber \\&\quad \;+\left. \left. \frac{\partial X_i}{\partial x}(l_{i-1},t) \varphi _{i,n}(l_{i-1})\right] \right\} . \end{aligned}$$

(34)

Next, define $D_{1}\frac{\partial X_1(l_1,t) }{\partial x}=D_{2}\frac{\partial X_2(l_1,t) }{\partial x}:= v_{12}(t)$. Substituting (31), (32) and (34) into (33), the transformed layer equations are then given by

$$\begin{aligned}&{^{C}_0D_t^{(\alpha _1,\rho _1)}} {\widetilde{X}}_1 =-D_{1} \lambda _{1,n}^2 {\widetilde{X}}_1 + D_1f_L(t)\frac{d \varphi _{1,n}(l_0)}{dx}\\&+v_{12}(t)\varphi _{1,n}(l_1)+S_{a_1}\langle 1,\varphi _{1,n}\rangle +S_{b_1}\, {\widetilde{X}}_1,\\&{^{C}_0D_t^{(\alpha _2,\rho _2)}} {\widetilde{X}}_2 =-D_{2} \lambda _{2,n}^2 {\widetilde{X}}_2 - v_{12}(t) \varphi _{2,n}(l_1)\\&-D_2f_R(t)\frac{d \varphi _{2,n}(l_2)}{dx}+ S_{a_2} \langle 1,\varphi _{2,n} \rangle + S_{b_2} \, {\widetilde{X}}_2, \end{aligned}$$

together with the transformed initial conditions ${\widetilde{X}}_i(\lambda _{i,n},0) = \langle X_{i,0}(x),\varphi _{i,n} \rangle :={\widetilde{X}}_{i,0}, i=1,2$. We now apply the Laplace transform in time and denote $\overline{{\widetilde{X}}}_i(\lambda _{i,n},s) = {\mathcal {L}}\left\{ {\widetilde{X}}_i(\lambda _{i,n},t) \right\} ,i=1,2$ to obtain

$$\begin{aligned}&(s+\rho _1)^{\alpha _1} \overline{{\widetilde{X}}}_1 -(s+\rho _1)^{\alpha _1-1} {\widetilde{X}}_{1,0} \\&= -D_{1} \lambda _{1,n}^2 \overline{{\widetilde{X}}}_1 + D_1\bar{f_L}(s)\frac{d \varphi _{1,n}(l_0)}{dx}\\&\quad +{\bar{v}}_{12}(s) \varphi _{1,n}(l_1)+\frac{S_{a_1}{{\mathbb {I}}_{1,n}}}{s} + S_{b_1} \overline{{\widetilde{X}}}_1, \\&(s+\rho _2)^{\alpha _2} \overline{{\widetilde{X}}}_2 -(s+\rho _2)^{\alpha _2-1} {\widetilde{X}}_{2,0} \\&= -D_{2} \lambda _{2,n}^2 \overline{{\widetilde{X}}}_2-{\bar{v}}_{12}(s) \varphi _{2,n}(l_1)-D_2{\bar{f}}_{R}(s)\frac{d \varphi _{2,n}(l_2)}{dx}\\&\quad +\frac{S_{a_2}{{\mathbb {I}}_{2,n}}}{s} + S_{b_2} \overline{{\widetilde{X}}}_2, \end{aligned}$$

which can be rearranged into the form

$$\begin{aligned} \overline{{\widetilde{X}}}_1(\lambda _{1,n},s)=&\frac{(s+\rho _1)^{\alpha _1-1}{\widetilde{X}}_{1,0}}{\eta _{1,n}(s)} \nonumber \\&+\frac{1}{\eta _{1,n}(s)}\left[ D_1\bar{f_L}(s)\frac{d \varphi _{1,n}(l_0)}{dx} +{\bar{v}}_{12}(s) \varphi _{1,n}(l_1)\right] \nonumber \\&+\frac{S_{a_1}{\mathbb {I}}_{1,n}}{s\eta _{1,n}(s)}, \end{aligned}$$

(35)

$$\begin{aligned} \overline{{\widetilde{X}}}_2(\lambda _{2,n},s)=&\frac{(s+\rho _2)^{\alpha _2-1}{\widetilde{X}}_{2,0}}{\eta _{2,n}(s)} \nonumber \\&-\frac{1}{\eta _{2,n}(s)} \left[ {\bar{v}}_{12}(s)\varphi _{2,n}(l_1)+D_2{\bar{f}}_{R}(s)\frac{d \varphi _{2,n}(l_2)}{dx}\right] \nonumber \\&+\frac{S_{a_2}{\mathbb {I}}_{2,n}}{s\eta _{2,n}(s)}, \end{aligned}$$

(36)

where ${\mathbb {I}}_{i,n}=\langle 1,\varphi _{i,n} \rangle$ and $\eta _{i,n}(s)=(s+\rho _i)^{\alpha _i}+D_{i} \lambda _{i,n}^{2}-S_{b_i}$, $i=1,2$. In order to calculate the unknown interfacial flux value ${\bar{v}}_{12}(s)$, we need the boundary condition at the interface:

$$\begin{aligned} {\overline{X}}_1(l_1,s) ={\overline{X}}_2(l_1,s). \end{aligned}$$

Noting that ${\overline{X}}_i(x,s)=\sum _{n=0}^{\infty }\overline{{\widetilde{X}}}_i(\lambda _{i,n},s) \varphi _{i,n}(x)$, $i=1,2$, then we have

$$\begin{aligned} \sum _{n=0}^{\infty }\overline{{\widetilde{X}}}_1(\lambda _{1,n},s) \varphi _{1,n}(l_1)= \sum _{n=0}^{\infty }\overline{{\widetilde{X}}}_2(\lambda _{2,n},s) \varphi _{2,n}(l_1). \end{aligned}$$

(37)

Substituting expressions (35) and (36) into (37), (37) can be rearranged into the following form:

$$\begin{aligned}&\sum _{n=0}^\infty \left\{ \frac{\varphi ^2_{1,n}(l_1)}{\eta _{1,n}(s)} + \frac{\varphi ^2_{2,n}(l_1)}{\eta _{2,n}(s)} \right\} {\overline{v}}_{12}(s)\\&\quad =-\sum _{n=0}^\infty \left( \frac{(s+\rho _1)^{\alpha _1-1}{\widetilde{X}}_{1,0}}{\eta _{1,n}(s)}+\frac{D_1\bar{f_L}(s)}{\eta _{1,n}(s)}\frac{d \varphi _{1,n}(l_0)}{dx}\right. \\&\left. +\frac{S_{a_1}{\mathbb {I}}_{1,n}}{s\eta _{1,n}}\right) \varphi _{1,n}(l_1)\\&\qquad +\sum _{n=0}^\infty \left( \frac{(s+\rho _2)^{\alpha _2-1}{\widetilde{X}}_{2,0}}{\eta _{2,n}(s)} -\frac{D_2{\bar{f}}_{R}(s)}{\eta _{2,n}(s)}\frac{d \varphi _{2,n}(l_2)}{dx}\right. \\&\qquad \left. +\frac{S_{a_2}{\mathbb {I}}_{2,n}}{s\eta _{2,n}(s)} \right) \varphi _{2,n}(l_1), \end{aligned}$$

which can be solved for ${\overline{v}}_{12}(s)$ at a given value of s. Applying the inverse Laplace transform and solved the solutions numerically within each layer, we can obtain

$$\begin{aligned} {{\widetilde{X}}}_i(\lambda _{i,n},t)&={\mathcal {L}}^{-1}\left\{ {\overline{{\widetilde{X}}}}_i(\lambda _{i,n},s) \right\} \\&=\frac{1}{2 \pi i} \int _{\Gamma } e^{st} \overline{{{\widetilde{X}}}}_i(\lambda _{i,n},s) \, ds\\&= \frac{1}{2 \pi i} \int _{\Gamma } \frac{e^z}{t} \overline{{{\widetilde{X}}}}_i(\lambda _{i,n},z/t) \, dz \approx \\&\quad -2\Re \left( \sum _{k=1}^{K/2}c_{2k-1} \frac{ \overline{{{\widetilde{X}}}}_i(\lambda _{i,n},z_{2k-1}/t)}{t}\right) , \end{aligned}$$

where $z=st$, $c_{2k-1}$ and $z_{2k-1}$ are the residues and poles of the near-best minimax approximation of $e^z$ on the negative real line by rational functions of type (K, K) as computed by the Carathéodor–Fejér method [38]. The final solution in each layer can be written as

$$\begin{aligned} X_i(x,t)=\sum _{n=0}^{\infty }{\widetilde{X}}(\lambda _{i,n},t) \varphi _{i,n}(x),\quad i=1,2. \end{aligned}$$

Fig. 8

The comparison between the numerical solution and the semi-analytical solution at different times for problem (28)–(32), where the parameters chosen are $\alpha _1=0.9$, $\alpha _2=0.8$, $\rho _1=0.1$, $\rho _2=0.5$, $D_1=0.25$, $D_2=0.5$, $S_{a_1}=S_{a_2}=0.1$, $S_{b_1}=S_{b_2}=-0.1$

Fig. 9

The impacts of the fractional index $\alpha$ (Figure a) and tempered parameter $\rho$ (Figure b) on the diffusion profile at $t=0.1$ with $D_1=0.25$, $D_2=0.5$, $S_{a_1}=S_{a_2}=0.1$, $S_{b_1}=S_{b_2}=-0.1$, a $\alpha _1=1$, $\rho _1=0$, $\rho _2=0.0$, b $\alpha _1=1$, $\rho _1=0$, $\alpha _2=0.9$

5.2 Numerical solution

Firstly, we do the grid partition. Define $t_n=T\left( \frac{n}{N} \right) ^r$, $r=\min \left\{ \frac{2(2-\alpha _1)}{\alpha _1} ,\frac{2(2-\alpha _2)}{\alpha _2} \right\}$, $n=0,1,2,\ldots ,N$, $\tau _n=t_n-t_{n-1}$, $n=1,2,\ldots ,N$. Define $x_i=l_0+ih$, $i=0,1,2,\ldots ,M$, where $h=\frac{l_2-l_0}{M}$ is the uniform spatial step. Divide the grid points into two parts: on layer 1 $\{ x_0,x_1,\ldots ,x_{M_1}\}$ and on layer 2 $\{ x_{M_1}, x_{M_1+1}, \ldots ,x_{M}\}$. Applying the L1 scheme and central difference scheme, we can obtain the numerical solution of the two-layered problem:

$$\begin{aligned}&{_0^C{\mathbb {D}}_t^{(\alpha _1,\rho _1)}}X^n_{1,i}=D_1\delta _x^2X^n_{1,i}+S_{a_1}+S_{b_1}X^n_{1,i},\\&\quad 1\le i\le M_1,~1\le n\le N,\\&{_0^C{\mathbb {D}}_t^{(\alpha _2,\rho _2)}}X^n_{2,i}=D_2\delta _x^2X^n_{2,i}+S_{a_2}+S_{b_2}X^n_{2,i},\\&\quad M_1\le i\le M,~1\le n\le N. \end{aligned}$$

Exploiting the boundary conditions (32) to overlap the value in $x=l_1$, the matrix form in terms of the solution on $[l_0,l_2]$ can be derived, which can be solved using the general iterative method.

5.3 Numerical examples

We consider the two-layered problem on a domain [0, 1] with $l_{1}=0.5$ and initial conditions $X_{1,0}(x)=X_{2,0}(x)=1$ and boundary conditions $f_{L}(t)=f_{R}(t)=0$. Figure 8 shows a comparison between the semi-analytical solution and the numerical solution at different times, from which there is perfect agreement between the semi-analytical solution and the numerical solution. It demonstrates that both methods work well. Next, we consider a scenario where one layer is normal material characterised by $\alpha _{1}=1$, $\rho _{1}=0$, while the other layer exhibits memory ($\alpha _{2}\ne 1$). Figure 9a presents the effects of the fractional index $\alpha _{2}$ on the solution profile. It can be seen that a small $\alpha _{2}$ can accelerate the diffusion. Figure 9b displays the impacts of the tempered parameter $\rho _{2}$ on the solution profile. Similarly, the tempered parameter $\rho _{2}$ can further promote the decay and with a moderate parameter pair $\alpha _{2}=0.9$, $\rho _{2}=0.5$; it decays faster than that only with a small $\alpha _{2}$. We can see that the influence of the tempered parameter on the two-layered problem is similar to that on a homogeneous medium.

6 Conclusions

In this paper, we extend two classical numerical schemes, the L1 scheme on graded mesh and the WSGL formula with correction terms, to deal with the benchmark problem with a tempered operator. Both schemes are effective. In addition, a fast algorithm for the time tempered Caputo derivative is developed to reduce the running time significantly. Furthermore, the tempered operator is applied to different models to investigate the tempered solution behaviour. An important finding is that, compared with the fractional index, the tempered parameter could further accelerate the diffusion, and the tempered model with two parameters $\alpha$ and $\rho$ is more flexible. In the future, we will explore high-dimensional tempered diffusion problems in heterogeneous media.