Efficient numerical method for multi-term time-fractional diffusion equations with Caputo-Fabrizio derivatives

1.   Introduction

Fractional differential equations have wide applications in various fields of science, including physics, economics, engineering, chemistry, biology and others^[1,2,3,4,5]. There are many kinds of definitions for the fractional derivatives, the most used fractional derivatives are the Riemann-Liouville fractional derivative and the Caputo fractional derivative ^[6,7]. However, both of these operators still present challenges in practical applications. To be more precise, the Riemann-Liouville derivative of a constant is non-zero and the Laplace transform of this derivative contains terms that lack physical significance. The Caputo fractional derivative has successfully addressed both issues, however, its definition involves a singular kernel which poses challenges in analysis and computation. Caputo and Fabrizio ^[8] have proposed a novel definition of the fractional derivative with a smooth kernel, referred to as the Caputo and Fabrizio (CF) derivatives, which present distinct representations for the temporal and spatial variables. The representation in time variable is suitable to use the Laplace transform, and the spatial representation is more convenient to use the Fourier transform. Although there is ongoing debate regarding the mathematical properties of fractional derivatives with non-singular kernels ^[9,10], numerous scholars remain interested in studying differential equations involving such derivatives due to their nice performance in various applications. Considering the CF derivative offers two primary advantages: 1) The utilization of a regular kernel in non-local systems is motivated by its potential to accurately depict material heterogeneities and fluctuations of various scales, which cannot be adequately captured by classical local theories or fractional models with singular kernels, see, e.g., ^[8,11]; 2) CF derivatives have numerical advantages. As we know, the truncation error of the numerical calculation for fractional operators with singular kernels is typically dependent on the order $\alpha$. For instance, in the case of the Caputo fractional derivative, employing classical L1 discretization results in an error of order $2-\alpha$, which becomes highly unfavorable when $\alpha\approx1$. In order to enhance actuarial accuracy, the utilization of higher-order methods will lead to an increase in computational complexity, particularly for problems with high dimensions. However, within the same approximation framework, the CF derivative has a higher truncation error, see Remark 2.2. Further properties and diverse applications of this fractional derivative can be found in various references, such as ^{[11,12,13,14,15,16,17]}.

Let $T, L > 0$, $\Lambda: = (0, S)$. In this paper, we are concerned with the numerical approximation of the multi-term time-fractional diffusion equation

with the initial conditions

and the boundary condition

where

$0 < \alpha_1\leq\alpha_2\leq\cdots\leq\alpha_{n-1}\leq\alpha_n < 1$, and $d_i\geq0, \; i = 1, 2, \cdots, n, n\in \mathbb{N}$. $\varphi(x)$ and $f(x, t)$ are given sufficiently smooth functions in their respective domains. In addition, ${}_0^{\mathrm{CF}}\!D_t^{\alpha_i}u(x, t)$ is the Caputo-Fabrizio derivative operator of order $\alpha_i$ ^[8,11] defined as

If $n = 1$, then $(1.1)$–$(1.3)$ reduces to the single-term time-fractional diffusion equation. The model of $(1.1)$–$(1.3)$, which describes the temporal flow of water within a leaky aquifer at various scales ^[12,18], as well as the electro-magneto-hydrodynamic flow of non-Newtonian biofluids with heat transfer ^[19], etc. For the well-posedness of $(1.1)$–$(1.3)$, we refer to, e.g., ^[17,20,21].

Many researchers have explored the numerical approximation of both single-term and multi-term time fractional diffusion equations. In ^[22], Liu et al. proposed a finite difference method for solving time-fractional diffusion equations in both space and time domains. Lin and Xu ^[23] utilized a finite difference scheme in time and Legendre spectral methods in space to numerically solve the time-fractional diffusion equations. Subsequently, Li and Xu ^[24] improved upon their previous work by proposing a space-time spectral method for these equations. For the numerical treatment of multi-term time-fractional diffusion equations, ^[25] proposed a fully-discrete schemes for one- and two-dimensional multi-term time fractional sub-diffusion equations. These schemes combine the compact difference method for spatial discretization with L1 approximation for time discretization. The Galerkin finite element method and the spectral method were introduced in ^[26] and ^[27,28], respectively. Zhao et al. ^[29] developed a fully-discrete scheme for a class of two-dimensional multi-term time-fractional diffusion equations with Caputo fractional derivatives, utilizing the finite element method in spatial direction and classical L1 approximation in temporal direction. Akman et al. ^[30] proposed a numerical approximation called the L1-2 formula for the Caputo-Fabrizio derivative using quadratic interpolation. In ^[31], finite difference/spectral approximations for solving two-dimensional time CF fractional diffusion equation were proposed and analyzed. Later, a second order scheme ^[32] was devised for addressing this problem. A compact alternating direction implicit (ADI) difference scheme was proposed by ^[33] for solving the two-dimensional time-fractional diffusion equation.

Simulating models with fractional derivatives presents a challenge due to their non-locality, which significantly impedes algorithm efficiency and necessitates greater memory storage compared to traditional local models. In particular, for fractional models, the computational complexity of obtaining an approximate solution is $O(N_T^2)$ and the required memory storage is $O(N_T)$, which contrasts with local models that have a complexity of $O(N_T)$ and require a memory storage of $O(1)$, where $N_T$ denotes the total number of time steps, see, e.g., ^[23,31,32]. To address this issue, several researchers have proposed efficient algorithms for computing the derivatives of Riemann-Liouville, Caputo, and Riesz fractional operators, see e.g., ^{[34,35,36,37]} and the references therein. Recently, a fast compact finite difference method for quasi-linear time-fractional parabolic equations is presented and analyzed in ^[38]. Then, ^[39] proposed a fast second-order numerical scheme for approximating the Caputo-Fabrizio fractional derivative at node $t = t_{k+1/2}$ with computational complexity of $O(N_T)$ and memory storage of $O(1)$.

Inspired by the above mentioned, we extend finite difference/spectral approximations for the multi-term Caputo-Fabrizio time-fractional diffusion equation (1.1)–(1.3). Firstly, we present a L1 formula for the Caputo-Fabrizio derivative. In this context, we introduce two discrete fractional differential operators, namely $L_t^{\alpha}$ and $F_t^{\alpha}$, which are essentially equivalent. However, $F_t^{\alpha}$ effectively utilizes the exponential kernel and incurs lower storage and computational costs compared to $L_t^{\alpha}$. The idea of this approach is essentially identical to that of reference ^[39], albeit with a slightly different formulation in our case; specifically, the approximation is centered at point $t = t_k$ and presented in a more concise manner. The error bounds associated with these two operators will be examined in detail. Secondly, we develop a semi-discrete scheme based on finite difference method for multi-term time-fractional derivatives, with complete proofs of its unconditional stability and convergence rate. A detailed error analysis is carried out for the semi-discrete problem, showing that the temporal accuracy is second order. Finally, we present the fully-discrete scheme based on the Legendre spectral collocation method for spatial discretization. We will investigate both the convergence order of this method and its implementation efficiency, while providing a rigorous proof of its spectral convergence in this paper.

The rest of this paper is organized as follows. In Section 2, a semi-discrete scheme is proposed for (1.1)–(1.3) based on fast L1 finite difference scheme. The stability and convergence analysis of the semi-discrete scheme is presented. In Section 3, we construct a Legendre spectral collocation method for the spatial discretization of the semi-discrete scheme. Error estimates are provided for the full discrete problem. Some numerical results are reported in Section 4. Finally, the conclusions are given in Section 5.

2.   Semi-discretization

Define $t_k: = k\Delta t$, $k = 0, 1, \cdots, N_T$, where $\Delta t: = T/N_T$ is the time step.

2.1.   Fast L1 formula for Caputo-Fabrizio derivative

We first give L1 approximation for fractional Caputo-Fabrizio derivative of function $h(t)$ defined by

In order to simplify the notations, we denote $h(t_k): = h^k$ for $0\leq k\leq M_t$. The L1 formula is obtained by substituting the linear Lagrange interpolation of $h(t)$ into (2.1). Precisely, the linear approximation of the function $h(t)$ on $[t_{j-1}, t_j]$ is written as

and the error in the approximation is

Then we define the discrete fractional differential operator $L_t^\alpha$ by

where $b_{j, k}: = \sigma_{j, k}-\sigma_{j-1, k}$ and

The right hand side of (2.4) involves a sum of all previous solutions $\left\{h^j\right\}_{j = 0}^k$, which reflects the memory effect of the non-local fractional derivative. Thus it requires on average $O\left(N_T\right)$ storage and the total computational cost is $O\left(N_T^2\right)$ with $N_T$ the total number of time steps. This makes both the computation and memory expensive, specially in the case of long time integration. In order to overcome this difficulty, we propose a further approach to the fractional derivative. The idea consists in first splitting the convolution integral in (2.1) into a sum of history part and local part as follows:

Note that a comparable treatment is employed in reference ^[34]. Then the history part $C_h(t_k)$ can be rewritten as

hence we have

Using the simple recurrence relation (2.5), we define the discrete fractional differential operator $F_t^{\alpha}$ by

It is not difficult to see that $F_t^{\alpha}h^k = L_t^{\alpha}h^k$ for $1\leq k\leq N_T$. Comparing $L_t^\alpha h^k$ in (2.4) with $F_t^\alpha h^k$ in (2.6) and (2.7), the former requires all the previous time step values of $h(t)$ while the latter only needs $h^{k-1}$, $h^{k}$ and $F_t^{\alpha}h^{k-1}$. This implies that approximating ${}_0^{\mathrm{CF}}\!D_t^{\alpha}h^k$ by $F_t^\alpha h^k$ considerably reduces the storage and computational costs as compared to $L_t^\alpha h^k$. Roughly speaking, replacing $L_t^\alpha h^k$ by $F_t^\alpha h^k$ allows to reduce the storage cost from $O\left(N_T\right)$ to $O\left(1\right)$, and the computational cost from $O\left(N_T^2\right)$ to $O\left(N_T\right)$.

Remark 2.1. The fast algorithm of Caputo derivative in ^[34] should be noted to retain an additional truncation error $\varepsilon$, whereas the fast algorithm of CF derivative does not introduce this error. Furthermore, it is worth mentioning that other algorithms, such as parallel computational methods ^[40], result in an augmented spatial complexity.

The following lemma provides an error bound for approximation $F_t^\alpha h^k$.

Lemma 2.1. Suppose that $h(t)\in \mathbb{C}^2[0, T]$. For any $0 < \alpha < 1$, let

Then

Proof. We consider proving the following estimate by mathematical induction:

First we have

where $\eta_1\in(0, t_1)$. Therefore, (2.8) holds for $j = 1$. Now suppose that (2.8) holds for $j = k-1$, we need to prove that it holds also for $j = k$. Similar to the proof of $\left|R_1\right|$, we can easily get that

By combining (2.5) and (2.7), we obtain

The estimate (2.8) is proved. Hence

which prove the conclusion of the lemma. □

Remark 2.2. The second rate of convergence of L1 formula has been proven in ^[30] by different methods, here we obtained identical results herein. Note that the rate of convergence of L1 formula for classical Caputo fractional derivative with order $\alpha$ is $2-\alpha$, this result seems reasonable since Caputo-Fabrizio derivative has smooth kernel.

2.2.   Discretization in time

We denote $u^k: = u(x, t_k)$ and $f^k(x): = f(x, t_k)$. Then from (2.6), (2.7) and Lemma 2.1, the time fractional derivative (1.5) at $t = t_k$ can be approximated by

where

Then Eq (1.1) can be rewritten as

where

with $C_{u, \alpha_i} > 0$ for $i = 1, 2, \cdots, n$. Notice that

we denote

the above equations are recast into

Let $u^k$ be the approximation for $u^k(x)$, and $f^k: = f^k(x)$. Then the semi-discrete problem of Eq (1.1) can be written as

where

Moreover, by utilizing relation

we can easily derive an alternative formulation of (2.15)–(2.18) as follows

where

Remark 2.3. Since $L_t^{\alpha_i} = F_t^{\alpha_i}$, (2.19)–(2.22) can be also obtained by using $L_t^{\alpha_i}$ in Eq (1.1). It is noteworthy that Eqs (2.15)–(2.18) offer computational advantages over Eqs (2.19)–(2.22). This is primarily attributed to the straightforward recurrence relation presented in Eqs (2.6) and (2.7). However, (2.19)–(2.22) is more appropriate for our analysis than (2.15)–(2.18), hence it play a crucial role in the subsequent sections.

Theorem 2.1. Let $\widetilde{R}_{k}$ be defined by (2.12), then there exists a constant $\widetilde{c} > 0$ such that

Proof. Without loss of generality, we assume that $\Delta t\in(0, 1)$. By the definition of $R_{k}$ and the inequalities of (2.10), we have

On the other hand, since

for $i = 1, 2, \cdots, n$, we get

This implies that

Therefore, there exists a constant $\widetilde{c} > 0$ such that

which prove the conclusion of the lemma. □

Lemma 2.2. Let the coefficients $b_{j, k}^{(\alpha_i)}$ be defined by (2.9), then for every $i$,

Proof. $b_{j, k}^{(\alpha_i)}\in(0, 1)$ can be easily obtained by the definition of $\sigma_{j, k}^{(\alpha_i)}$ and the monotone property of the function $g(x) = \exp(x)$. Finally, note that

Using the above equalities and the fact $\exp\left(-\frac{\alpha_i\Delta t}{1-\alpha_i}\right)\in(0, 1)$ completes the proof of the lemma.□

Remark 2.4. (2.26) gives a easy way to compute all the coefficients $b_{j, k}^{(\alpha_i)}$.

Lemma 2.3. Let the coefficients $\zeta_{j, k}$ be defined by (2.23), then

Proof. By Lemma 2.2, and the definition of $\zeta_{j, k}$, we can readily arrive at these conclusions.□

2.3.   Stability and convergence analysis of the semi-discrete scheme

To discuss the stability and convergence of the semi-discrete scheme, we introduce functional spaces equipped with standard norms and inner products that will be utilized subsequently. Let $L^2(\Lambda)$ is the space of measurable functions whose square is Lebesgue integrable in $\Lambda$. Then

The inner products of $L^2(\Lambda)$ and $H^1(\Lambda)$ are defined, respectively, by

and the corresponding norms by

The norm $\|\cdot\|_m$ of the space $H^m(\Omega)$ $(m\in\mathrm{N})$ is defined by

In this paper, instead of using the above standard $H^1$-norm, we prefer to define $\|\cdot\|_1$ by

It is widely acknowledged that the standard $H^1$-norm and the norm defined by (2.27) are equivalent; therefore, we will adopt the latter in subsequent discussions.

The variational (weak) formulation of the Eqs (2.15) and (2.16)/(2.20), subject to the boundary condition (2.18), can be expressed as finding $u^{k}\in H_0^1(\Lambda)$ such that for $\forall v\in H_0^1(\Lambda)$

where $u^0: = u(x, 0)$.

For the semi-discretized problems (2.28) and (2.29), we can establish a stability result as follows.

Theorem 2.2. The semi-discretized problems (2.28) and (2.29) is unconditionally stable in the sense that for all $\Delta t > 0$, it holds

Proof. By mathematical induction. First of all, when $k = 1$, we have

Notice that $\|v\|_0\leq\|v\|_1$, taking $v = u^1$ and using the Cauchy-Schwarz inequality, we obtain immediately

Now, suppose

Taking $v = u^{k}$ in (2.29) gives

Hence, by using (2.30) and Lemma 2.3, we have

Thus, the proof is completed.□

Remark 2.5. In the proof of the following theorem, we will demonstrate that $\zeta_{1, k}^{-1}$ is bounded. As shown in Eq (2.25), $|\kappa| = O(1)$. Therefore, it follows that $\kappa\zeta_{1, k}^{-1}$ is also bounded.

We now conduct an error analysis for the solution of the semi-discretized problem.

Theorem 2.3. Assuming $n\neq 1$ and $0 < \alpha_1\leq\alpha_2\leq\cdots\leq\alpha_{n-1}\leq\alpha_n < 1$. Let $u^k(x)$ be the exact solution of (1.1)–(1.3), $\{u^k\}_{k = 0}^{N_T}$ be the solution of semi-discretized problems (2.28) and (2.29) with initial condition $u^0: = u^0(x)$, then the following error estimates hold:

where the constant $\widetilde{c}$ is defined in (2.24) and $S$ is the length of $\Lambda$.

Proof. We shall establish the following estimate through a process of mathematical induction:

Let $\bar{e}^k = u^k(x)-u^k$, $k = 1, 2, \cdots, N_t$. By combining (2.13) and (2.28), the error $\bar{e}^1$ satisfies

Taking $v = \bar{e}^1$ yields $\|\bar{e}^1\|_1^2\leq\|\widetilde{R}_{1}\|_0\|\bar{e}^1\|_0$. This, in conjunction with (2.24), yields

Therefore, (2.32) holds for $k = 1$. Assuming that (2.32) holds for all $k = 1, 2, \cdots, l-1$, it is necessary to demonstrate its validity for $k = l$. By combining (2.13), (2.14) and (2.29), for $\forall v\in H_0^1(\Lambda)$, we have

Let $v = \bar{e}^{l}$ in the above equation, then

Using the induction assumption and Lemma 2.3, we derive

Next we show that $\zeta_{1, k}^{-1}$ is bounded. Considering that function $h(\alpha) = :-\frac{\alpha}{1-\alpha}$ is decreasing on $(0, 1)$, we have

by combining Eq (2.26). Therefore,

Consequently we obtain, for all $k$ such that $k\Delta t\leq T$,

□

3.   Full discretization

3.1.   A shifted Legendre collocation method in space

We shall begin by providing a comprehensive overview of fundamental definitions and properties pertaining to Legendre Gauss-type quadratures. Let $\mathbb{P}_N(\Lambda)$ denote the space of algebraic polynomials of degree less than or equal to $N$ with respect to variable $x$, and $L_N(x)$ be the Legendre polynomial of degree $N$ on the interval $[-1, 1]$. Then the discrete space, denoted by $\mathbb{P}_N^0(\Lambda): = \mathbb{P}_N(\Lambda)\cap H_0^1(\Lambda)$.

Let $\pi_N^1$ be the $H^1$-orthogonal projection operator from $H_0^1(\Lambda)$ into $\mathbb{P}_N^0(\Lambda)$, associated to the norm $\|\cdot\|_1$ defined in (2.27), that is, for all $\psi\in H_0^1(\Lambda)$, define $\pi_N^1\psi\in \mathbb{P}_N^0(\Lambda)$, such that, $\forall v_N\in \mathbb{P}_N^0(\Lambda)$,

From ^[41], the following estimate of projection holds:

Define the Legendre-Gauss-Lobatto nodes and weights as $\xi_p$ and $\omega_p$, $p = 0, 1, \cdots, N$, $N\geq 1$, where $\{\xi_k\}_{k = 0}^{N}$ are the zeroes of $(1-x^2)L'_N(x)$, and

Moreover, the following quadrature holds

The discrete inner product and norm defined as follow, for any continuous functions $\phi, \psi\in \mathbb{C}\left([-1, 1]\right)$,

From ^[42], the discrete norm $\|\cdot\|_N$ is equivalent to the standard $L^2$-norm in $\mathbb{P}_N([-1, 1])$. If we denote $\{\hat{\xi}_{p}\}_{p = 0}^N$ and $\{\hat{\omega}_{p}\}_{p = 0}^N$ as the nodes and weights of shifted Legendre-Gauss-Lobatto quadratures on $\overline{\Lambda}$, then one can easily show that

Thus, we define the discrete inner product and norm on $\overline{\Lambda}$ as follows

It is not difficult to obtain that

and

We introduce the operator of interpolation at the $N+1$ shifted Legendre-Gauss-Lobatto nodes, denoted by $I_N$, i.e., $\forall \psi\in C(\overline{\Lambda})$, $I_N\psi\in \mathbb{P}_N(\Lambda)$, such that

The interpolation error estimate (see ^[41]) is

Now consider the spectral discretization to the problems $(2.28)$ and $(2.29)$ as follows: find $u_N^{k}\in \mathbb{P}_N^0(\Lambda)$, such that

where

For $\{u_N^j\}_{j = 0}^{k-1}$ given, the well-posedness of the problem (3.7) is guaranteed by the well-known Lax-Milgram Lemma.

3.2.   Convergence analysis of the full discretization scheme

To simplify matters, we present the semi-discretized problems (2.28) and (2.29) in a compact form: find $u^{k}\in H_0^1(\Lambda)$, such that

where

We denote by $\|\cdot\|_{1, \widehat{N}}$ the norm associated to the bilinear form $A_{\widehat{N}}(\cdot, \cdot)$:

It follows from (3.3) that for all $v_N\in \mathbb{P}_N(\Lambda)$ the discrete norm $\|\cdot\|_{1, \widehat{N}}$ is equivalent to the norm $\|\cdot\|_1$ defined in (2.27).

Theorem 3.1. Assuming $n\neq 1$ and $0 < \alpha_1\leq\alpha_2\leq\cdots\leq\alpha_{n-1}\leq\alpha_n < 1$. Let $\{u_N^k\}_{k = 0}^{N_T}$ is the solution of the problem (3.7) with the initial condition $u_N^0$ taken to be $I_Nu^0(x)$, $\{u^k\}_{k = 0}^{N_T}$ the solution of the semi-discretized problems (2.28) and (2.29). Suppose that $u^k\in H^m(\Lambda)\cap H_0^1(\Lambda)$ with $m > 1$, for $k = 1, 2, \cdots, N_T$, then there exists a constant $\overline{c} > 0$ such that

Proof. For any $\forall v_{N-1}\in \mathbb{P}_{N-1}^0(\Lambda)$, denote $\rho_N^{k}: = u_N^{k}-v_{N-1}$. It is direct to check that

By virtue of (3.4) gives

hence

For the last term, by definition, we have

and

It is known that the following result holds (see e.g., ^[41,42]): $\forall g\in H^m(\Lambda)$, $m\geq 1$, $\forall \delta_N\in \mathbb{P}_N(\Lambda)$,

Thus for $\forall g\in H^m(\Lambda)$, $m\geq 1$, $\forall g_N, \rho_N\in \mathbb{P}_N(\Lambda)$ we have

Applying the above results to (3.9) and (3.10), we obtain

and

Let $\varepsilon_N^j: = u^j-u_N^j$, using (3.8) and the norm equivalence, for $\forall v_{N-1}\in \mathbb{P}_{N-1}^0(\Lambda)$, we have

and

By triangular inequality

for $\forall v_{N-1}\in \mathbb{P}_{N-1}^0(\Lambda)$, we obtain

and

The above estimate specially holds for $v_{N-1} = \pi_{N-1}^1u^{k}$, which implies

and

Similar to the proof of Theorem 2.2, we can immediately get the following conclusions:

Notice that

and the boundedness of $\kappa$ and $\zeta_{1, k}^{-1}$, then there exists a constant $\overline{c} > 0$ such that

□

4.   Numerical validation

4.1.   Implementation

We provide a comprehensive account of the implementation of problem (3.7) using the shifted Legendre collocation method.

Considering problem (3.7), we express the function $u_N^{k+1}$ in terms of the Lagrangian interpolants based on the shifted Legendre-Gauss-Lobatto points $\hat{\xi}_{i}, \; i = 0, 1, \cdots, N$, i.e.,

where $c_{i}^{k}: = u_N^{k}(\hat{\xi}_i)$, unknowns of the discrete solution. $h_i(x)$ is the Lagrangian polynomials defined in $\Lambda$, which satisfies

where $\delta_{ij}$ is the Kronecker symbols. Taking (4.1) into (3.7), and notice that the homogeneous Dirichlet boundary condition (1.3), then choosing each test function $v_N$ to be $h_l(x)$ $(l = 1, 2, \cdots, N-1)$, we have

Define the matrices

Then, we obtain the matrix representation of the above equation in the following form:

The linear system (4.2) can be solved in particular by the LU factorization or other related computational techniques.

Finally, we discuss about the calculation of $Q^{k}$. When $k = 0$, the initial condition $u_N^0$ taken to be

It is not difficult to see that

which implies that $u_N^0$ satisfies interpolation condition (3.5). When $k\geq 1$, suppose that

then

Furthermore,

In a word, we can easily obtain $Q^{k}$ at each iteration of time-step.

4.2.   Numerical results

We present a series of numerical results to validate our theoretical propositions.

Firstly, to investigate the computational performance of two discrete fractional differential operators $L_t^{\alpha}$ and $F_t^{\alpha}$, we test three examples from ^[30]. Denote $\beta: = \frac{\alpha}{1-\alpha}$.

Example 4.1. Consider the function $h(t) = t^m$ $(m\geq 1, \; m\in \mathbb{N})$, the Caputo-Fabrizio fractional derivative of order $\alpha$ with $0 < \alpha < 1$ of $h(t)$ is written as

Example 4.2. Consider the function $h(t) = \cos(\omega t)$, the Caputo-Fabrizio fractional derivative of order $\alpha$ with $0 < \alpha < 1$ of $h(t)$ is written as

Example 4.3. Consider the function $h(t) = \exp(\omega t)$, the Caputo-Fabrizio fractional derivative of order $\alpha$ with $0 < \alpha < 1$ of $h(t)$ is written as

The proofs based on the method of integration by parts can be found in ^[8,30].

We choose $m = 4$ in Example 4.1 and $\omega = 5$ in Examples 4.2 and 4.3, and set $\alpha = 0.5$, $t\in[0, 2]$. Define the errors

for $L_t^{\alpha}$ and $F_t^{\alpha}$ operators, respectively, where $N_T$ is the last time step. Tables 1–3 give the numerical results of approximation error and CPU time with three examples. Here CPU time represents the total computation time, that is, the whole time for computing the approximations of Caputo-Fabrizio fractional derivatives at every time step. The convergence rates in Tables are given by

Tables 1–3 demonstrate that the errors of both $L_t^{\alpha}$ and $F_t^{\alpha}$ approximations are virtually identical, as a result of their equivalence ($F_t^{\alpha}h^k = L_t^{\alpha}h^k$). Moreover, both approximations have achieved second-order convergence of error, as stated in Lemma 2.1. However, we observe that the CPU time of $F_t^{\alpha}$ approximation increases linearly with respect to $N_T$, while the $L_t^{\alpha}$ approximation increases almost quadratically. This suggests that the $F_t^{\alpha}$ operator holds promise as it requires less storage and incurs lower computational costs than the $L_t^{\alpha}$ operator during computation.

Secondly, we provide preliminary computational findings to demonstrate the efficacy of the finite difference/shifted Legendre collocation method (abbreviated as FCM).

Example 4.4. Consider the following three-term time-fractional diffusion equations:

where $\Lambda = (0, \pi)$, $0 < \alpha_1\leq\alpha_2\leq\alpha_3 < 1$, and

with $\beta_i = \frac{\alpha_i}{1-\alpha_i}$, $i = 1, 2, 3$. The exact solution of the Eq (4.3) is $u(x, t) = (1+t^2)\sin x$, which is sufficiently smooth. In our experiments, we set the parameters $d_1 = 1, d_2 = 2, d_3 = 3$.

The following error norms have been used as the error indicator:

We test Example 4.4 with four cases:

Case 1. $\alpha_1 = \alpha_2 = \alpha_3 = 0.5$, then (4.3) reduces to the single-term time-fractional diffusion equation;

Case 2. $\alpha_1 = 0.3$, $\alpha_2 = 0.5$, $\alpha_3 = 0.7$;

Case 3. $\alpha_1 = 0.2$, $\alpha_2 = 0.3$, $\alpha_3 = 0.4$;

Case 4. $\alpha_1 = 0.6$, $\alpha_2 = 0.7$, $\alpha_3 = 0.8$.

In time discretization, we use $F_t^{\alpha_i}$ operator. Table 4 shows the errors and temporal accuracy of FCM with polynomial degree $N = 20$ at $t_k = 1$ for different cases of $\alpha_i$. Here convergence rates are given by

It can be observed that the FCM exhibits a second-order temporal convergence rate, which is consistent with our theoretical analysis.

Next, we check the spatial accuracy with respect to the polynomial degree $N$. In order to avoid the contamination of temporal error, we need fix the time step $\Delta t$ sufficiently small. Here we take $\Delta t = 10^{-6}$, and terminate computing at $t_k = 0.01$ for saving time. Figure 1 shows the errors with respect to polynomial degree $N$ at $t_k = 0.01$ in semi-log scale. Evidently, the spatial discretization exhibits exponential convergence as demonstrated by the nearly linear curves depicted in this figure. The aforementioned is known as spectral accuracy as expected since the exact solution is a sufficiently smooth function with respect to the space variable.

To further verify the numerical validity, we finally test a two-dimensional problem.

Example 4.5. Consider the following three-term time-fractional diffusion equations:

where $\Omega = [0, 1]\times[0, 1]$, $0 < \alpha_1\leq\alpha_2\leq\alpha_3\leq1$, and

The exact solution of the Eq (4.4) is $u(x, y, t) = \sin(2\pi xy)\exp(t)$, which is sufficiently smooth. In our experiments, we set the parameters $d_1 = 1, d_2 = 2, d_3 = 3$.

In this case, we denote $\{\hat{\xi}_{pq}\}_{p, q = 0}^N: = \{(\hat{\xi}_p, \hat{\xi}_q)\}_{p, q = 0}^N$ and $\{\hat{\omega}_{pq}\}_{p, q = 0}^N: = \{\hat{\omega}_p\hat{\omega}_q\}_{p, q = 0}^N$ as the nodes and weights of shifted Legendre-Gauss-Lobatto quadratures on $\overline{\Omega}$. Then we express the function $u_N^{k+1}$ in terms of the two-dimensional Lagrangian interpolants based on the shifted Legendre-Gauss-Lobatto points $\hat{\xi}_{ij}, \; i, j = 0, 1, \cdots, N$,

where $c_{ij}^{k+1}: = u_N^{k+1}(\hat{\xi}_i, \hat{\xi}_j)$, unknowns of the discrete solution. $h_i(x)$ and $h_j(y)$ are the Lagrangian polynomials defined in $I_x: = [0, 1]$ and $I_y: = [0, 1]$, i.e.,

where $\delta_{il}$ and $\delta_{js}$ are the Kronecker symbols. A linear system such as (4.2) can be readily derived. Here we take $\Delta t = 10^{-6}$. Figure 2 shows the errors with respect to polynomial degree $N$ in semi-log scale. Thanks to the fast scheme (2.7), a small time step does not significantly escalate the computational burden in the time direction, thereby the proposed method is effective even for handling high-dimensional problems.

5.   Conclusions

In this work, we have developed a fully discrete scheme for the multi-term time-fractional diffusion equations with Caputo-Fabrizio derivatives. The proposed approach utilizes the finite difference method to approximate multi-term fractional derivatives in time and employs the Legendre spectral collocation method for spatial discretization. Specifically, we use the exponential property of Caputo-Fabrizio derivative to give a recursive difference calculation scheme, which offers benefits in terms of computational complexity and storage capacity. The proposed scheme has been proved to be unconditionally stable and convergent with order $O\left(\left(\Delta t\right)^{2}+N^{-m}\right)$. Numerical results show good agreement with the theoretical analysis. Due to its high resolution feature in spectral approximation, the proposed method can be extended to handle multi-term time-fractional diffusion equations in higher spatial dimensions.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the Natural Science Foundation of Fujian Province (Grant No. 2022J05198), the Foundation of Fujian Provincial Department of Education (Grant No. JAT210293), and the Natural Science Foundation of Fujian University of Technology (Grant No. GY-Z21069).

We are thankful for the anonymous reviewers for their valuable comments and suggestions.

Conflict of interest

The authors declare no conflicts of interest in this manuscript.