# A Fast Compact Finite Difference Method for Fractional Cattaneo Equation Based on Caputo–Fabrizio Derivative

Authors: Haili Qiao; Zhengguang Liu; Aijie Cheng
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The Cattaneo equations with Caputo–Fabrizio fractional derivative are investigated. A compact finite difference scheme of Crank–Nicolson type is presented and analyzed, which is proved to have temporal accuracy of second order and spatial accuracy of fourth order. Since this derivative is defined with an integral over the whole passed time, conventional direct solvers generally take computational complexity ofOMN2 and require memory of OMN, with M and N the number of space steps and time steps, respectively. We develop a fast evaluation procedure for the Caputo–Fabrizio fractional derivative, by which the computational cost is reduced to OMN operations; meanwhile, only OM memory is required. In the end, several numerical experiments are carried out to verify the theoretical results and show the applicability of the fast compact difference procedure.


---

## Body

## 1. Introduction

Fractional diffusion equations have become a strong and forceful tool to describe the phenomenon of anomalous diffusion, and more research works have been obtained in the last decades [1–6]. However, since the fractional derivative is nonlocal and has weak singularity, it is impossible to solve fractional diffusion equations analytically in most cases. Instead, seeking numerical solutions is becoming an indispensable tool for research work about fractional equations.Different from the traditional derivative of the integer order, the fractional derivative depends on the total information in the correlative region, and this is the so-called nonlocal properties. Just because of this, it consumes computational time extremely to solve fractional equations. We hope to develop effective numerical schemes, which not only have better stability and higher accuracy but also require less storage memory and save computational cost.About stability and convergence analysis of the numerical schemes for fractional equations, the readers can refers to [7, 8] for spatial fractional order equation, [9–18] for temporal fractional diffusion equations, and [19–22] for space-time-fractional equations. About the complexity, i.e., storage requirement and computation cost of an algorithm, researchers devote themselves to reduce storage requirement and computational time by analyzing the particular structure of coefficient matrices arising from the discretization system or reutilizing the intermediate data skillfully. We call these algorithms fast methods, including fast finite difference methods [23–28], fast finite element methods [29], and fast collocation methods [30, 31]. A fast method for Caputo fractional derivatives is proposed [32, 33]. Lu et al. [34] presented a fast method of approximate inversion for triangular Toeplitz tridiagonal block matrix, which is successfully applied to the fractional diffusion equations. Comparatively, there is less research work about the fast method for temporal fractional derivative than that for spatial fractional operators.A time-fractional Cattaneo equation is considered with the following form:(1)∂ux,t∂t+∂αux,t∂tα=∂2ux,t∂x2+fx,t,x,t∈Ω×0,T,ux,0=ϕx,∂u∂tt=0=ψx,x∈Ω,ux,t=0,x∈∂Ω,t>0,where 1<α<2; Ω=a,b for one-dimensional case, and Ω=a,b×c,d for two-dimensional case; fx,t is the source term; ϕxandψx are the prescribed functions for initial conditions; and ∂αu/∂tα is a new Caputo fractional derivative without singular kernel, which is defined in the next section.Our purpose is to establish a fast finite difference scheme of high order for this equation. We will extract the recursive relation between thek+1 time step and the k time step of the finite difference solution. The computational work is significantly reduced from OMN2 to OMN, and the memory requirement from OMN to OM, where M and N are the total numbers of points for space steps and time steps, respectively. For improving the accuracy, a compact finite difference scheme is established. Theoretical analysis shows that the fast compact difference scheme has spatial accuracy of fourth order and temporal accuracy of second order. Several numerical experiments are implemented, which verify the effectiveness, applicability, and convergence rate of the proposed scheme.This paper is organized as follows: some definitions and notations are prepared in Section2. The compact finite difference scheme is described and then the stability and convergence rates are rigorously analyzed for the scheme in Section 3. The compact finite difference scheme is extended to the case of two space dimensions in Section 4. Fast evaluation and efficient storage are established skillfully in Section 5. Some numerical experiments are carried out in Section 6. In the end, we summarize the major contribution of this paper in Section 7.


## 2. Some Notations and Definitions

We provide some definitions which will be used in the following analysis.First, let us recall the usual Caputo fractional derivative of orderα with respect to time variable t, which is given by(2)DtαaCux,t=1Γn−α∑att−sn−α−1unx,sds,n−1<α<n.By replacing the kernel functiont−s−α with the exponential function exp−αt−s/1−α and 1/Γ1−α with Mα/1−α, Caputo and Fabrizio [35] proposed the following definition of fractional time derivative.Definition 1 (see [35]).
Letu⋅,t∈H1a,b, b>a, α∈0,1, then the new Caputo derivative of the fractional order is defined as(3)Dtαux,t=Mα1−α∫atu′x,sexp−αt−s1−αds,where Mα is a normalization function satisfying M0=M1=1. When the function u does not belong to H1a,b, then this derivative can be reformulated as(4)Dtαux,t=αMα1−α∫atux,t−ux,sexp−αt−s1−αds.Definition 2 (see [36]).
The above new Caputo derivative of order0<α<1 can also be reformulated as(5)Dtα0CFux,t=11−α∑0tu′x,sexp−αt−s1−αds.Definition 3 (see [35, 36]).
Letu⋅,t∈H1a,b, if n≥1, and α∈0,1, the fractional time derivative Dtα0CFux,t of order n+α is defined by(6)Dtn+α0CFux,t=Dtα0CFDtn0CFux,t.
Particularly, for1<α<2, we have(7)Dtα0CFux,t=Mα2−α∑0tu″x,sexp1−αt−s2−αds.Remark 1.
The Caputo–Fabrizio (CF) operator was proposed with a nonsingular kernel for describing material heterogeneities that do not exhibit power-law behavior [35].Remark 2.
An open discussion is ongoing about the mathematical construction of the CF operator. Ortigueira and Tenreiro Machado [37] indicated that the CF fractional derivative is neither a fractional operator nor a derivative operator, the authors of [38, 39] showed that this operator cannot describe dynamic memory, and Giusti [40] indicated that this operator can be expressed as an infinite linear combination of Riemann–Liouville integrals with integer powers. As responses to these criticisms, Atangana and Gómez-Aguilar [41] pointed out the need to account for a fractional calculus approach without an imposed index law and with nonsingular kernels. Furthermore, Hristov [42] indicated that the CF operator is not applicable for explaining the physical examples in [37, 40]; instead, he suggested that the CF operator can be used for the analysis of materials that do not follow a power-law behavior. The authors of [43] believe that models with CF operators produce a better representation of physical behaviors than do integer-order models, providing a way to model the intermediate (between elliptic and parabolic or between parabolic and hyperbolic) behaviors.
To obtain the accuracy of the fourth order in spatial directions, the following lemma is necessary.Lemma 1 (see [44]).
Denoteθs=1−s35−31−s2. If fx∈C6a,b,h=b−a/M,xi=a+ih0≤i≤M, then, it holds that(8)112f″xi−1+10f″xi+f″xi+1=1h2fxi−1−2fxi+fxi+1+h4360∫01f6xi−sh+f6xi+shθsds,1≤i≤M−1.


## 3. Compact Finite Difference Scheme for One-Dimensional Fractional Cattaneo Equation

In order to construct the finite difference schemes, the intervala,b is divided into subintervals with xi=a+ih,0≤i≤M, and 0,T is discretized with tk=kΔt,0≤k≤N, where h=b−a/M and Δt=T/N are the spatial grid size and temporal step size, respectively. Denote Ωh=xi:0≤i≤M, ΩΔt=tk:0≤k≤N, then Ωh×ΩΔt becomes a discretization of the practical computational domain a,b×0,T. The values of the function u at the grid points are denoted as ujk=uxj,tk, and the approximate solution at the point xj,tk is denoted as ujk.DenoteVh=vv=v0,v1,…,vM. We also introduce the following notations for any mesh function v∈Vh:(9)δxvj−1/2=vj−vj−1h,δx2vj=vj+1−2vj+vj−1h2,1≤j≤M−1,and define the average operator(10)Avj=112vj+1+10vj+vj−1,1≤j≤M−1,vj,j=0,M.It is easy to see that(11)Avj=I+h212δx2vj,where I is the identical operator. We also denote Av=Av1,Av2,…,AvM for vector v=v1,v2,…,vM, and Au,v=Au,v.For any two gird functionsu,v∈Vh0=vv∈Vh,v0=vM=0, the discrete inner products and norms are defined as(12)u,v=h∑j=1M−1ujvj,u,v=h∑j=1M−1δxuj−1/2δxvj−1/2,u22=u,u,u12=u,u,u∞=max1≤j≤M−1uj.By summation by parts, it is easy to see that(13)δx2u,v=−δxu,δxv=−u,v=u,δx2v.For the average operatorA, define(14)Av,v≜Av,v=vA2.Additionally, letVΔt=vv=v0,v1,⋯,vN be the space of grid function defined on ΩΔt. For any function v∈VΔt, a difference operator is introduced as follows:(15)δtvk=vk−vk−1Δt.


### 3.1. Compact Finite Difference Scheme

We will consider the time-fractional Cattaneo equation equipped with the Caputo–Fabrizio derivative. Vivas-Cruz et al. [43] gave the theoretical analysis of a model of fluid flow in a reservoir with the Caputo–Fabrizio operator. They proved that this model cannot be used to describe nonlocal processes since it can be represented as an equivalent differential equation with a finite number of integer-order derivatives.The finite difference methods usually lead to stencils through the whole history passed by the solution which consume too much computational work. In this paper, we will establish a high-order finite difference scheme and propose a procedure to reduce the computational cost. In [43], the authors proposed a recurrence formula of discretized CF operator and obtained an algorithm which can be considered a stencil with a one-step expression without the need of integrals over the whole history. It seems that the procedure in our paper and the algorithm in [43] are different in approach but equally satisfactory in result.For obtaining effective approximation with high order, we introduce the numerical discretization for the fractional Cattaneo equation by means of compact finite difference methods.At the nodexi,tk+1/2, the differential equation takes the following form:(16)∂u∂txi,tk+1/2+∂αu∂tαxi,tk+1/2=∂2u∂x2xi,tk+1/2+fxi,tk+1/2,1≤i≤M−1,1≤k≤N−1.The approximation of the fractional derivative is given by [45](17)Dtα0CFuxi,tk+1/2=1α−1ΔtM0δtuik+1−∑n=1kMk−n−Mk−n+1δtuij−Mkψi+Rik+1/2,with truncation error Rik+1/2=OΔt2, and(18)Mn=exp1−α2−αΔtn−exp1−α2−αΔtn+1.Furthermore, by Lemma1 and equation (15), the space and time derivative are approximated by(19)∂2u∂x2xi,tk+1/2=12δx2uik+1+δx2uikA+OΔt2+h4,(20)∂u∂txi,tk+1/2=δtuik+1+OΔt2.Substituting (17) and (19)∼(20) into (16), we get(21)Aδtuik+1+1α−1ΔtAM0δtuik+1−∑n=1kMk−n−Mk−n+1δtuin−Mkψi=12δx2uik+1+δx2uik+Afik+1/2+Rik+1/2,1≤i≤M−1,0≤k≤N−1,and there exists a constant C, depending on the function u and its derivatives such that(22)Rik+1/2≤CΔt2+h4.By the initial and boundary value conditions, we have(23)ui0=ϕi,1≤i≤M−1,u0k=uMk=0,0≤k≤N.A compact finite difference scheme can be established by omitting the truncation termRik+1/2 and replacing the exact solution uik in equation (21) with numerical solution uik:(24)α−1ΔtAδtuik+1+M0Aδtuik+1−α−1Δt2δx2Uk+1=α−1Δt2δx2Uk+∑n=1kMk−n−Mk−n+1Aδtuin+MkAψi+α−1ΔtAfik+1/2,1≤i≤M−1,1≤k≤N−1,ui0=ϕi,1≤i≤M−1,u0k=uMk=0,0≤k≤N.


### 3.2. Stability Analysis and Optimal Error Estimates

#### 3.2.1. Stability Analysis

The following Lemma aboutMn is useful for the analysis of stability.Lemma 2 (see [45]).
For the definition ofMn, Mn>0 and Mn+1<Mn, ∀n≤k, are held.Multiplyinghδtuik+1 on both sides of equation (24) and summing up with respect to i from 1 to M−1, the following equation is obtained:(25)α−1ΔtδtUk+1A2+M0δtUk+1A2−α−12δx2Uk+1,Uk+1+α−12δx2Uk+1,Uk=α−12δx2Uk,Uk+1−α−12δx2Uk,Uk+∑n=1kMk−n−Mk−n+1AδtUn,δtUk+1+MkAψ,δtUk+1+α−1ΔtAfk+1/2,δtUk+1.Observing equation (13), we have(26)δx2Uk+1,Uk+1=−δxUk+1,δxUk+1=−δxUk+112≤0,δx2Uk,Uk=−δxUk,δxUk=−δxUk12≤0,δx2Uk+1,Uk=−δxUk+1,δxUk=−Uk+1,Uk=δx2Uk,Uk+1.By the triangle inequality and Lemma2, we obtain(27)∑n=1kMk−n−Mk−n+1AδtUn,δtUk+1+MkAψ,δtUk+1≤∑n=1k12Mk−n−Mk−n+1δtUnA2+δtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−n−Mk−n+1δtUnA2+12M0−MkδtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−nδtUnA2−∑n=1k12Mk−n+1δtUnA2+12M0δtUk+1A2+12MkψA2.Combining equation (25) with (26)∼(27), we get(28)∑n=1k+1Mk−n+1δtUnA2+α−1δxUk+112≤∑n=1kMk−nδtUnA2+α−1δxUk12+MkψA2+α−1Δtfk+1/2A2.Let(29)QUk=∑n=1kMk−nδtUnA2+α−1δxUk12.Summing up with respect tok from 0 to N−1 leads to(30)QUN≤QU0+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.The initial conditionU0=ϕ implies that QU0=α−1δxϕ12, and then(31)QUN≤α−1δxϕ12+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.Theorem 1.
For scheme (24), we have the following stable conclusion:(32)QUm≤α−1δxϕ12+∑k=0m−1MkψA2+α−1Δt∑k=0m−1fk+1/2A2,∀0≤m≤N.


#### 3.2.2. Optimal Error Estimates

Combining equations (21) and (23) with (24), we get an error equation as follows:(33)α−1ΔtAδteik+1+M0Aδteik+1−α−1Δt2δx2ek+1=α−1Δt2δx2ek+∑n=1kMk−n−Mk−n+1Aδtein+α−1ΔtRik+1/2,where Rik+1/2=OΔt2+h4 and eik=uik−uik, ∀k≥0.Multiplyinghδteik+1 on both sides of equation (33) and summing up with respect to i from 1 to M−1, we get(34)α−1Δtδtek+1A2+M0δtek+1A2+α−12δxek+112=α−12δxek12+∑n=1kMk−n−Mk−n+1Aδten,δtek+1+α−1ΔtRk+1/2,δtek+1.By the triangle inequality and Lemma 2, we obtain(35)∑n=1kMk−n−Mk−n+1Aδten,δtek+1≤∑n=1k12Mk−n−Mk−n+1δtenA2+δtek+1A2=∑n=1k12Mk−n−Mk−n+1δtenA2+12M0−Mkδtek+1A2≤∑n=1k12Mk−nδtenA2−∑n=1k12Mk−n+1δtenA2+12M0δtek+1A2.Combining equation (34) with (35), we have(36)∑n=1k+1Mk−n+1δtenA2+α−1δxek+112≤∑n=1kMk−nδtenA2+α−1δxek12+α−1ΔtRk+1/22.By the definition ofQ in stability analysis, the inequality (36) can be rearranged as(37)Qek+1≤Qek+α−1ΔtRk+1/22.Summing up with respect tok from 0 to N−1, we get(38)QeN≤Qe0+CΔt2+h42.Observing that the initial errore0=0 implies Qe0=0. Then, we have(39)α−1δxeN12≤QeN≤CΔt2+h42,∑n=1NMN−nδtenA2≤QeN≤CΔt2+h42.Theorem 2.
Suppose that the exact solution of the fractional Cattaneo equation is smooth sufficiently, then there exists a positive constant C, independent ofh,k, and Δt such that(40)ek1≤CΔt2+h4,∀0≤k≤N,where eik=uik−uik.


## 3.1. Compact Finite Difference Scheme

We will consider the time-fractional Cattaneo equation equipped with the Caputo–Fabrizio derivative. Vivas-Cruz et al. [43] gave the theoretical analysis of a model of fluid flow in a reservoir with the Caputo–Fabrizio operator. They proved that this model cannot be used to describe nonlocal processes since it can be represented as an equivalent differential equation with a finite number of integer-order derivatives.The finite difference methods usually lead to stencils through the whole history passed by the solution which consume too much computational work. In this paper, we will establish a high-order finite difference scheme and propose a procedure to reduce the computational cost. In [43], the authors proposed a recurrence formula of discretized CF operator and obtained an algorithm which can be considered a stencil with a one-step expression without the need of integrals over the whole history. It seems that the procedure in our paper and the algorithm in [43] are different in approach but equally satisfactory in result.For obtaining effective approximation with high order, we introduce the numerical discretization for the fractional Cattaneo equation by means of compact finite difference methods.At the nodexi,tk+1/2, the differential equation takes the following form:(16)∂u∂txi,tk+1/2+∂αu∂tαxi,tk+1/2=∂2u∂x2xi,tk+1/2+fxi,tk+1/2,1≤i≤M−1,1≤k≤N−1.The approximation of the fractional derivative is given by [45](17)Dtα0CFuxi,tk+1/2=1α−1ΔtM0δtuik+1−∑n=1kMk−n−Mk−n+1δtuij−Mkψi+Rik+1/2,with truncation error Rik+1/2=OΔt2, and(18)Mn=exp1−α2−αΔtn−exp1−α2−αΔtn+1.Furthermore, by Lemma1 and equation (15), the space and time derivative are approximated by(19)∂2u∂x2xi,tk+1/2=12δx2uik+1+δx2uikA+OΔt2+h4,(20)∂u∂txi,tk+1/2=δtuik+1+OΔt2.Substituting (17) and (19)∼(20) into (16), we get(21)Aδtuik+1+1α−1ΔtAM0δtuik+1−∑n=1kMk−n−Mk−n+1δtuin−Mkψi=12δx2uik+1+δx2uik+Afik+1/2+Rik+1/2,1≤i≤M−1,0≤k≤N−1,and there exists a constant C, depending on the function u and its derivatives such that(22)Rik+1/2≤CΔt2+h4.By the initial and boundary value conditions, we have(23)ui0=ϕi,1≤i≤M−1,u0k=uMk=0,0≤k≤N.A compact finite difference scheme can be established by omitting the truncation termRik+1/2 and replacing the exact solution uik in equation (21) with numerical solution uik:(24)α−1ΔtAδtuik+1+M0Aδtuik+1−α−1Δt2δx2Uk+1=α−1Δt2δx2Uk+∑n=1kMk−n−Mk−n+1Aδtuin+MkAψi+α−1ΔtAfik+1/2,1≤i≤M−1,1≤k≤N−1,ui0=ϕi,1≤i≤M−1,u0k=uMk=0,0≤k≤N.


## 3.2. Stability Analysis and Optimal Error Estimates

### 3.2.1. Stability Analysis

The following Lemma aboutMn is useful for the analysis of stability.Lemma 2 (see [45]).
For the definition ofMn, Mn>0 and Mn+1<Mn, ∀n≤k, are held.Multiplyinghδtuik+1 on both sides of equation (24) and summing up with respect to i from 1 to M−1, the following equation is obtained:(25)α−1ΔtδtUk+1A2+M0δtUk+1A2−α−12δx2Uk+1,Uk+1+α−12δx2Uk+1,Uk=α−12δx2Uk,Uk+1−α−12δx2Uk,Uk+∑n=1kMk−n−Mk−n+1AδtUn,δtUk+1+MkAψ,δtUk+1+α−1ΔtAfk+1/2,δtUk+1.Observing equation (13), we have(26)δx2Uk+1,Uk+1=−δxUk+1,δxUk+1=−δxUk+112≤0,δx2Uk,Uk=−δxUk,δxUk=−δxUk12≤0,δx2Uk+1,Uk=−δxUk+1,δxUk=−Uk+1,Uk=δx2Uk,Uk+1.By the triangle inequality and Lemma2, we obtain(27)∑n=1kMk−n−Mk−n+1AδtUn,δtUk+1+MkAψ,δtUk+1≤∑n=1k12Mk−n−Mk−n+1δtUnA2+δtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−n−Mk−n+1δtUnA2+12M0−MkδtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−nδtUnA2−∑n=1k12Mk−n+1δtUnA2+12M0δtUk+1A2+12MkψA2.Combining equation (25) with (26)∼(27), we get(28)∑n=1k+1Mk−n+1δtUnA2+α−1δxUk+112≤∑n=1kMk−nδtUnA2+α−1δxUk12+MkψA2+α−1Δtfk+1/2A2.Let(29)QUk=∑n=1kMk−nδtUnA2+α−1δxUk12.Summing up with respect tok from 0 to N−1 leads to(30)QUN≤QU0+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.The initial conditionU0=ϕ implies that QU0=α−1δxϕ12, and then(31)QUN≤α−1δxϕ12+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.Theorem 1.
For scheme (24), we have the following stable conclusion:(32)QUm≤α−1δxϕ12+∑k=0m−1MkψA2+α−1Δt∑k=0m−1fk+1/2A2,∀0≤m≤N.


### 3.2.2. Optimal Error Estimates

Combining equations (21) and (23) with (24), we get an error equation as follows:(33)α−1ΔtAδteik+1+M0Aδteik+1−α−1Δt2δx2ek+1=α−1Δt2δx2ek+∑n=1kMk−n−Mk−n+1Aδtein+α−1ΔtRik+1/2,where Rik+1/2=OΔt2+h4 and eik=uik−uik, ∀k≥0.Multiplyinghδteik+1 on both sides of equation (33) and summing up with respect to i from 1 to M−1, we get(34)α−1Δtδtek+1A2+M0δtek+1A2+α−12δxek+112=α−12δxek12+∑n=1kMk−n−Mk−n+1Aδten,δtek+1+α−1ΔtRk+1/2,δtek+1.By the triangle inequality and Lemma 2, we obtain(35)∑n=1kMk−n−Mk−n+1Aδten,δtek+1≤∑n=1k12Mk−n−Mk−n+1δtenA2+δtek+1A2=∑n=1k12Mk−n−Mk−n+1δtenA2+12M0−Mkδtek+1A2≤∑n=1k12Mk−nδtenA2−∑n=1k12Mk−n+1δtenA2+12M0δtek+1A2.Combining equation (34) with (35), we have(36)∑n=1k+1Mk−n+1δtenA2+α−1δxek+112≤∑n=1kMk−nδtenA2+α−1δxek12+α−1ΔtRk+1/22.By the definition ofQ in stability analysis, the inequality (36) can be rearranged as(37)Qek+1≤Qek+α−1ΔtRk+1/22.Summing up with respect tok from 0 to N−1, we get(38)QeN≤Qe0+CΔt2+h42.Observing that the initial errore0=0 implies Qe0=0. Then, we have(39)α−1δxeN12≤QeN≤CΔt2+h42,∑n=1NMN−nδtenA2≤QeN≤CΔt2+h42.Theorem 2.
Suppose that the exact solution of the fractional Cattaneo equation is smooth sufficiently, then there exists a positive constant C, independent ofh,k, and Δt such that(40)ek1≤CΔt2+h4,∀0≤k≤N,where eik=uik−uik.


## 3.2.1. Stability Analysis

The following Lemma aboutMn is useful for the analysis of stability.Lemma 2 (see [45]).
For the definition ofMn, Mn>0 and Mn+1<Mn, ∀n≤k, are held.Multiplyinghδtuik+1 on both sides of equation (24) and summing up with respect to i from 1 to M−1, the following equation is obtained:(25)α−1ΔtδtUk+1A2+M0δtUk+1A2−α−12δx2Uk+1,Uk+1+α−12δx2Uk+1,Uk=α−12δx2Uk,Uk+1−α−12δx2Uk,Uk+∑n=1kMk−n−Mk−n+1AδtUn,δtUk+1+MkAψ,δtUk+1+α−1ΔtAfk+1/2,δtUk+1.Observing equation (13), we have(26)δx2Uk+1,Uk+1=−δxUk+1,δxUk+1=−δxUk+112≤0,δx2Uk,Uk=−δxUk,δxUk=−δxUk12≤0,δx2Uk+1,Uk=−δxUk+1,δxUk=−Uk+1,Uk=δx2Uk,Uk+1.By the triangle inequality and Lemma2, we obtain(27)∑n=1kMk−n−Mk−n+1AδtUn,δtUk+1+MkAψ,δtUk+1≤∑n=1k12Mk−n−Mk−n+1δtUnA2+δtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−n−Mk−n+1δtUnA2+12M0−MkδtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−nδtUnA2−∑n=1k12Mk−n+1δtUnA2+12M0δtUk+1A2+12MkψA2.Combining equation (25) with (26)∼(27), we get(28)∑n=1k+1Mk−n+1δtUnA2+α−1δxUk+112≤∑n=1kMk−nδtUnA2+α−1δxUk12+MkψA2+α−1Δtfk+1/2A2.Let(29)QUk=∑n=1kMk−nδtUnA2+α−1δxUk12.Summing up with respect tok from 0 to N−1 leads to(30)QUN≤QU0+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.The initial conditionU0=ϕ implies that QU0=α−1δxϕ12, and then(31)QUN≤α−1δxϕ12+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.Theorem 1.
For scheme (24), we have the following stable conclusion:(32)QUm≤α−1δxϕ12+∑k=0m−1MkψA2+α−1Δt∑k=0m−1fk+1/2A2,∀0≤m≤N.


## 3.2.2. Optimal Error Estimates

Combining equations (21) and (23) with (24), we get an error equation as follows:(33)α−1ΔtAδteik+1+M0Aδteik+1−α−1Δt2δx2ek+1=α−1Δt2δx2ek+∑n=1kMk−n−Mk−n+1Aδtein+α−1ΔtRik+1/2,where Rik+1/2=OΔt2+h4 and eik=uik−uik, ∀k≥0.Multiplyinghδteik+1 on both sides of equation (33) and summing up with respect to i from 1 to M−1, we get(34)α−1Δtδtek+1A2+M0δtek+1A2+α−12δxek+112=α−12δxek12+∑n=1kMk−n−Mk−n+1Aδten,δtek+1+α−1ΔtRk+1/2,δtek+1.By the triangle inequality and Lemma 2, we obtain(35)∑n=1kMk−n−Mk−n+1Aδten,δtek+1≤∑n=1k12Mk−n−Mk−n+1δtenA2+δtek+1A2=∑n=1k12Mk−n−Mk−n+1δtenA2+12M0−Mkδtek+1A2≤∑n=1k12Mk−nδtenA2−∑n=1k12Mk−n+1δtenA2+12M0δtek+1A2.Combining equation (34) with (35), we have(36)∑n=1k+1Mk−n+1δtenA2+α−1δxek+112≤∑n=1kMk−nδtenA2+α−1δxek12+α−1ΔtRk+1/22.By the definition ofQ in stability analysis, the inequality (36) can be rearranged as(37)Qek+1≤Qek+α−1ΔtRk+1/22.Summing up with respect tok from 0 to N−1, we get(38)QeN≤Qe0+CΔt2+h42.Observing that the initial errore0=0 implies Qe0=0. Then, we have(39)α−1δxeN12≤QeN≤CΔt2+h42,∑n=1NMN−nδtenA2≤QeN≤CΔt2+h42.Theorem 2.
Suppose that the exact solution of the fractional Cattaneo equation is smooth sufficiently, then there exists a positive constant C, independent ofh,k, and Δt such that(40)ek1≤CΔt2+h4,∀0≤k≤N,where eik=uik−uik.


## 4. Compact Finite Difference Scheme in Two Dimensions

In this section, the following fractional Cattaneo equation in two dimensions will be considered:(41)∂ux,y,t∂t+∂αux,y,t∂tα=∂2ux,y,t∂x2+∂2ux,y,t∂y2+fx,y,t,x,y,t∈a,b×c,d×0,T,ux,y,0=ϕx,y,∂u∂tt=0=ψx,y,a≤x≤b,c≤y≤d,ua,y,t=ub,y,t=ux,c,t=ux,d,t=0,t>0,where 1<α<2, fx,y,t is the source term, ϕx,yandψx,y are the given functions, and ∂αu/∂tα is defined by the new Caputo fractional derivative without singular kernel.In order to construct the finite difference schemes, the rectanglea,b×c,d is discretized with xi=a+ihx0≤i≤Mx and .yj=c+jhy0≤j≤My, and the time interval 0,T is discretized with tk=kΔt,0≤k≤N, where hx=b−a/Mx, hy=d−c/My, and Δt=T/N are the spatial grid and temporal step sizes, respectively. Denote Ωh=xi,yj:0≤i≤Mx,0≤j≤My and ΩΔt=tk:0≤k≤N, then Ωh×ΩΔt is a discretization of the physical computational domain a,b×c,d×0,T. ui,jk=uxi,yj,tk denotes the values of function u at the grid points, and ui,jk denotes the values of the numerical solution at the point xi,yj,tk.DenoteVh=vv=v0,0,v0,1,…,vMx,My and Vh0=vv∈Vh,v0,j=vMx,j=vi,0=vi,My=0,0≤i≤Mx,0≤j≤My. For any mesh function v∈Vh, we use the following notations:(42)δxvi−1/2,j=vi,j−vi−1,jhx,δx2vi,j=vi+1,j−2vi,j+vi−1,jhx2,δyvi,j−1/2=vi,j−vi,j−1hy,δy2vi,j=vi,j+1−2vi,j+vi,j−1hy2,1≤i≤Mx−1,1≤j≤My−1,and define the average operator(43)Axvi,j=112vi+1,j+10vi,j+vi−1,j,1≤i≤Mx−1,0≤j≤My,vi,j,i=0,Mx,0≤j≤My.Ayvi,j=112vi,j+1+10vi,j+vi,j−1,1≤j≤My−1,0≤i≤Mx,vi,j,j=0,My,0≤i≤Mx.It is clear that(44)Axvi,j=I+hx212δx2vi,j,Ayvi,j=I+hy212δy2vi,j.We also denote AxAyu,v=AxAyu,v. It is easy to see that AxAy=AyAx.For any gird functionu,v∈Vh0, the discrete inner product and norms are defined as follows:(45)u,v=hxhy∑i=1Mx−1∑j=1My−1ui,jvi,j,u22=u,u.For the average operatorAxAy, define(46)AxAyv,v≜AxAyv,v=vA2.


### 4.1. Compact Finite Difference Scheme

At the nodexi,yj,tk+1/2, the differential equation is rewritten as(47)∂u∂txi,yj,tk+1/2+∂αu∂tαxi,yj,tk+1/2=∂2u∂x2xi,yj,tk+1/2+∂2u∂y2xi,yj,tk+1/2+fxi,yj,tk+1/2,1≤i≤Mx−1,1≤j≤My−1,1≤k≤N−1.For the approximation of the time-fractional derivative, we have the following approximation [45]:(48)Dtα0CFuxi,yj,tk+1/2=1α−1ΔtM0δtui,jk+1−∑n=1kMk−n−Mk−n+1δtui,jn−Mkψi,j+Ri,jk+1/2,where the truncation error Ri,jk+1/2=OΔt2 and(49)Mn=exp1−α2−αΔtn−exp1−α2−αΔtn+1.Furthermore, we also have(50)∂2u∂x2xi,yj,tk+1/2=12δx2ui,jk+1+δx2ui,jkAx+OΔt2+hx4,(51)∂2u∂y2xi,yj,tk+1/2=12δy2ui,jk+1+δy2ui,jkAy+OΔt2+hy4,(52)∂u∂txi,yj,tk+1/2=δtui,jk+1+OΔt2.Substituting (48) and (50)∼(52) into (47) leads to(53)AxAyδtui,jk+1+1α−1τAxAyM0δtui,jk+1−∑n=1kMk−n−Mk−n+1δtui,jn−Mkψi,j=12δx2ui,jk+1+δx2ui,jk+12δy2ui,jk+1+δy2ui,jk+AxAyfi,jk+1/2+Ri,jk+1/21≤i≤Mx−1,1≤j≤My−1,0≤k≤N−1,and there exists a constant C, depending on the function u and its derivatives such that(54)Ri,jk+1/2≤CΔt2+hx4+hy4.By the initial and boundary conditions, we have(55)ui,j0=ϕi,j,1≤i≤Mx−1,1≤j≤My−1,u0,jk=uMx,jk=ui,0k=ui,Myk=0,0≤k≤N.Omitting the truncation errorRi,jk+1/2 and replacing the true solution ui,jk with numerical solution ui,jk, a compact finite difference scheme can be obtained as follows:(56)α−1ΔtAxAyδtUi,jk+1+M0AxAyδtUi,jk+1−Ayα−1Δt2δx2Ui,jk+1−Axα−1Δt2δy2Ui,jk+1=Ayα−1Δt2δx2ui,jk+Axα−1Δt2δy2ui,jk+∑n=1kMk−n−Mk−n+1AxAyδtui,jn+MkAxAyψi,j+α−1ΔtAxAyfi,jk+1/2,1≤i≤Mx−1,1≤j≤My−1,1≤k≤N−1.ui,j0=ϕi,j,1≤i≤Mx−1,1≤j≤My−1,u0,jk=uMx,jk=ui,0k=ui,Myk=0,0≤k≤N.


### 4.2. Stability Analysis and Optimal Error Estimates

#### 4.2.1. Stability Analysis

Definition 4 (see [46]).
For any gird functionu∈Vh0, define the norm(57)⦀∇hu⦀A2=δxu2−hy212δyδxu2+δyu2−hx212δxδyu2.
The lemmas below is useful in the subsequent analysis of stability.Lemma 3 (see [46]).
For any gird functionu∈Vh0, the following equation is held:(58)13u2≤uA2≤u2.Lemma 4 (see [46]).
For any gird functionu∈Vh0, the following equation is held:(59)23∇hu2≤⦀∇hu⦀A2≤∇hu2,Ayδx2un−1/2+Axδy2un−1/2,δtun−1/2=−12Δt⦀∇hun⦀A2−⦀∇hun−1⦀A2.where ∇hu2=δxu22+δyu22.MultiplyinghxhyδtUi,jk+1 on both sides of equation (56) and summing up w.r.t. i,j from 1 to Mx−1 and from 1 to My−1, respectively, the following equation is obtained:(60)α−1ΔtδtUk+1A2+M0δtUk+1A2=α−1ΔtAyδx2Uk+1/2+Axδy2Uk+1/2,δtUk+1+∑n=1kMk−n−Mk−n+1AyAxδtUn,δtUk+1+MkAyAxψ,δtUk+1+α−1ΔtAyAxfk+1/2,δtUk+1.Observing Lemma4, we have(61)Ayδx2Uk+1/2+Axδy2Uk+1/2,δtUk+1=−12Δt⦀∇huk+1⦀A2−⦀∇huk⦀A2,(62)α−1ΔtAyAxfk+1/2,δtUk+1≤α−1Δt2fk+1/2A2+δtuk+1A2.By the triangle inequality and Lemma2, we obtain(63)∑n=1kMk−n−Mk−n+1AyAxδtUn,δtUk+1+MkAyAxψ,δtUk+1≤∑n=1k12Mk−n−Mk−n+1δtUnA2+δtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−n−Mk−n+1δtUnA2+12M0−MkδtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−nδtUnA2−∑n=1k12Mk−n+1δtUnA2+12M0δtUk+1A2+12MkψA2.Combining equation (60) with (61)∼(63),we get(64)∑n=1k+1Mk−n+1δtUnA2+α−1⦀∇hUk+1⦀A2≤∑n=1kMk−nδtUnA2+α−1∇hUkA2+MkψA2+α−1Δt⦀fk+1/2⦀A2.Let(65)QUk=∑n=1kMk−nδtUnA2+α−1⦀∇hUk⦀A2.Summing up with respect tok from 0 to N−1, we get(66)QUN≤QU0+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.Noting thatU0=ϕ, we have QU0=α−1∇hϕA2. It follows that(67)QUN≤α−1⦀∇hϕ⦀A2+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.Theorem 3.
For the compact finite difference scheme (56), the following stability inequality holds:(68)QUm≤α−1∇hϕA2+∑k=0m−1MkψA2+α−1Δt∑k=0m−1fk+1/2A2,∀0≤m≤N.Similar to the stability, the convergence can also be analyzed.Theorem 4.
Suppose that the exact solution of the fractional Cattaneo equation is sufficiently smooth, then there exists a positive constantC independent of h,k, and Δt such that(69)ek1≤CΔt2+h4,∀0≤k≤N,where ei,jk=ui,jk−ui,jk and h=maxhx,hy.


## 4.1. Compact Finite Difference Scheme

At the nodexi,yj,tk+1/2, the differential equation is rewritten as(47)∂u∂txi,yj,tk+1/2+∂αu∂tαxi,yj,tk+1/2=∂2u∂x2xi,yj,tk+1/2+∂2u∂y2xi,yj,tk+1/2+fxi,yj,tk+1/2,1≤i≤Mx−1,1≤j≤My−1,1≤k≤N−1.For the approximation of the time-fractional derivative, we have the following approximation [45]:(48)Dtα0CFuxi,yj,tk+1/2=1α−1ΔtM0δtui,jk+1−∑n=1kMk−n−Mk−n+1δtui,jn−Mkψi,j+Ri,jk+1/2,where the truncation error Ri,jk+1/2=OΔt2 and(49)Mn=exp1−α2−αΔtn−exp1−α2−αΔtn+1.Furthermore, we also have(50)∂2u∂x2xi,yj,tk+1/2=12δx2ui,jk+1+δx2ui,jkAx+OΔt2+hx4,(51)∂2u∂y2xi,yj,tk+1/2=12δy2ui,jk+1+δy2ui,jkAy+OΔt2+hy4,(52)∂u∂txi,yj,tk+1/2=δtui,jk+1+OΔt2.Substituting (48) and (50)∼(52) into (47) leads to(53)AxAyδtui,jk+1+1α−1τAxAyM0δtui,jk+1−∑n=1kMk−n−Mk−n+1δtui,jn−Mkψi,j=12δx2ui,jk+1+δx2ui,jk+12δy2ui,jk+1+δy2ui,jk+AxAyfi,jk+1/2+Ri,jk+1/21≤i≤Mx−1,1≤j≤My−1,0≤k≤N−1,and there exists a constant C, depending on the function u and its derivatives such that(54)Ri,jk+1/2≤CΔt2+hx4+hy4.By the initial and boundary conditions, we have(55)ui,j0=ϕi,j,1≤i≤Mx−1,1≤j≤My−1,u0,jk=uMx,jk=ui,0k=ui,Myk=0,0≤k≤N.Omitting the truncation errorRi,jk+1/2 and replacing the true solution ui,jk with numerical solution ui,jk, a compact finite difference scheme can be obtained as follows:(56)α−1ΔtAxAyδtUi,jk+1+M0AxAyδtUi,jk+1−Ayα−1Δt2δx2Ui,jk+1−Axα−1Δt2δy2Ui,jk+1=Ayα−1Δt2δx2ui,jk+Axα−1Δt2δy2ui,jk+∑n=1kMk−n−Mk−n+1AxAyδtui,jn+MkAxAyψi,j+α−1ΔtAxAyfi,jk+1/2,1≤i≤Mx−1,1≤j≤My−1,1≤k≤N−1.ui,j0=ϕi,j,1≤i≤Mx−1,1≤j≤My−1,u0,jk=uMx,jk=ui,0k=ui,Myk=0,0≤k≤N.


## 4.2. Stability Analysis and Optimal Error Estimates

### 4.2.1. Stability Analysis

Definition 4 (see [46]).
For any gird functionu∈Vh0, define the norm(57)⦀∇hu⦀A2=δxu2−hy212δyδxu2+δyu2−hx212δxδyu2.
The lemmas below is useful in the subsequent analysis of stability.Lemma 3 (see [46]).
For any gird functionu∈Vh0, the following equation is held:(58)13u2≤uA2≤u2.Lemma 4 (see [46]).
For any gird functionu∈Vh0, the following equation is held:(59)23∇hu2≤⦀∇hu⦀A2≤∇hu2,Ayδx2un−1/2+Axδy2un−1/2,δtun−1/2=−12Δt⦀∇hun⦀A2−⦀∇hun−1⦀A2.where ∇hu2=δxu22+δyu22.MultiplyinghxhyδtUi,jk+1 on both sides of equation (56) and summing up w.r.t. i,j from 1 to Mx−1 and from 1 to My−1, respectively, the following equation is obtained:(60)α−1ΔtδtUk+1A2+M0δtUk+1A2=α−1ΔtAyδx2Uk+1/2+Axδy2Uk+1/2,δtUk+1+∑n=1kMk−n−Mk−n+1AyAxδtUn,δtUk+1+MkAyAxψ,δtUk+1+α−1ΔtAyAxfk+1/2,δtUk+1.Observing Lemma4, we have(61)Ayδx2Uk+1/2+Axδy2Uk+1/2,δtUk+1=−12Δt⦀∇huk+1⦀A2−⦀∇huk⦀A2,(62)α−1ΔtAyAxfk+1/2,δtUk+1≤α−1Δt2fk+1/2A2+δtuk+1A2.By the triangle inequality and Lemma2, we obtain(63)∑n=1kMk−n−Mk−n+1AyAxδtUn,δtUk+1+MkAyAxψ,δtUk+1≤∑n=1k12Mk−n−Mk−n+1δtUnA2+δtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−n−Mk−n+1δtUnA2+12M0−MkδtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−nδtUnA2−∑n=1k12Mk−n+1δtUnA2+12M0δtUk+1A2+12MkψA2.Combining equation (60) with (61)∼(63),we get(64)∑n=1k+1Mk−n+1δtUnA2+α−1⦀∇hUk+1⦀A2≤∑n=1kMk−nδtUnA2+α−1∇hUkA2+MkψA2+α−1Δt⦀fk+1/2⦀A2.Let(65)QUk=∑n=1kMk−nδtUnA2+α−1⦀∇hUk⦀A2.Summing up with respect tok from 0 to N−1, we get(66)QUN≤QU0+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.Noting thatU0=ϕ, we have QU0=α−1∇hϕA2. It follows that(67)QUN≤α−1⦀∇hϕ⦀A2+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.Theorem 3.
For the compact finite difference scheme (56), the following stability inequality holds:(68)QUm≤α−1∇hϕA2+∑k=0m−1MkψA2+α−1Δt∑k=0m−1fk+1/2A2,∀0≤m≤N.Similar to the stability, the convergence can also be analyzed.Theorem 4.
Suppose that the exact solution of the fractional Cattaneo equation is sufficiently smooth, then there exists a positive constantC independent of h,k, and Δt such that(69)ek1≤CΔt2+h4,∀0≤k≤N,where ei,jk=ui,jk−ui,jk and h=maxhx,hy.


## 4.2.1. Stability Analysis

Definition 4 (see [46]).
For any gird functionu∈Vh0, define the norm(57)⦀∇hu⦀A2=δxu2−hy212δyδxu2+δyu2−hx212δxδyu2.
The lemmas below is useful in the subsequent analysis of stability.Lemma 3 (see [46]).
For any gird functionu∈Vh0, the following equation is held:(58)13u2≤uA2≤u2.Lemma 4 (see [46]).
For any gird functionu∈Vh0, the following equation is held:(59)23∇hu2≤⦀∇hu⦀A2≤∇hu2,Ayδx2un−1/2+Axδy2un−1/2,δtun−1/2=−12Δt⦀∇hun⦀A2−⦀∇hun−1⦀A2.where ∇hu2=δxu22+δyu22.MultiplyinghxhyδtUi,jk+1 on both sides of equation (56) and summing up w.r.t. i,j from 1 to Mx−1 and from 1 to My−1, respectively, the following equation is obtained:(60)α−1ΔtδtUk+1A2+M0δtUk+1A2=α−1ΔtAyδx2Uk+1/2+Axδy2Uk+1/2,δtUk+1+∑n=1kMk−n−Mk−n+1AyAxδtUn,δtUk+1+MkAyAxψ,δtUk+1+α−1ΔtAyAxfk+1/2,δtUk+1.Observing Lemma4, we have(61)Ayδx2Uk+1/2+Axδy2Uk+1/2,δtUk+1=−12Δt⦀∇huk+1⦀A2−⦀∇huk⦀A2,(62)α−1ΔtAyAxfk+1/2,δtUk+1≤α−1Δt2fk+1/2A2+δtuk+1A2.By the triangle inequality and Lemma2, we obtain(63)∑n=1kMk−n−Mk−n+1AyAxδtUn,δtUk+1+MkAyAxψ,δtUk+1≤∑n=1k12Mk−n−Mk−n+1δtUnA2+δtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−n−Mk−n+1δtUnA2+12M0−MkδtUk+1A2+12MkψA2+δtUk+1A2=∑n=1k12Mk−nδtUnA2−∑n=1k12Mk−n+1δtUnA2+12M0δtUk+1A2+12MkψA2.Combining equation (60) with (61)∼(63),we get(64)∑n=1k+1Mk−n+1δtUnA2+α−1⦀∇hUk+1⦀A2≤∑n=1kMk−nδtUnA2+α−1∇hUkA2+MkψA2+α−1Δt⦀fk+1/2⦀A2.Let(65)QUk=∑n=1kMk−nδtUnA2+α−1⦀∇hUk⦀A2.Summing up with respect tok from 0 to N−1, we get(66)QUN≤QU0+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.Noting thatU0=ϕ, we have QU0=α−1∇hϕA2. It follows that(67)QUN≤α−1⦀∇hϕ⦀A2+∑k=0N−1MkψA2+α−1Δt∑k=0N−1fk+1/2A2.Theorem 3.
For the compact finite difference scheme (56), the following stability inequality holds:(68)QUm≤α−1∇hϕA2+∑k=0m−1MkψA2+α−1Δt∑k=0m−1fk+1/2A2,∀0≤m≤N.Similar to the stability, the convergence can also be analyzed.Theorem 4.
Suppose that the exact solution of the fractional Cattaneo equation is sufficiently smooth, then there exists a positive constantC independent of h,k, and Δt such that(69)ek1≤CΔt2+h4,∀0≤k≤N,where ei,jk=ui,jk−ui,jk and h=maxhx,hy.


## 5. Efficient Storage and Fast Evaluation of the Caputo–Fabrizio Fractional Derivative

Since time-fractional derivative operator is nonlocal, the traditional direct method for numerically solving the fractional Cattaneo equations generally requires totalOMN memory units and OMN2 computational complexity, where N and M are the total number of time steps and space steps, respectively.In this section, we develop a fast solution method for the finite difference scheme of the time-fractional Cattaneo equation.Let(70)Nnk=Mk−n−Mk−n+1=exp1−α2−αk−nΔt1−2exp1−α2−αΔt+exp1−α2−α2Δt,1≤n≤k,then(71)Nnk+1=Mk+1−n−Mk+1−n+1=exp1−α2−αk−nΔt1−2exp1−α2−αΔt+exp1−α2−α2Δtexp1−α2−αΔt=exp1−α2−αΔtNnk,Nk+1k+1=M0−M1=1−2exp1−α2−αΔt+exp1−α2−α2Δt.So we have(72)∑n=1k+1Nnk+1Aδtuin=exp1−α2−αΔt∑n=1kNnkAδtuin+Nk+1k+1Aδtuik+1,∑n=1k+1Nnk+1AyAxδtui,jn=exp1−α2−αΔt∑n=1kNnkAyAxδtui,jn+Nk+1k+1AyAxδtUi,jk+1.Remark.
We find that at thek-th level, only O1 operations are needed to compute the k-th level since the k−1-th level is known at that point. Thus, the total operations are reduced from ON2 to ON, and the memory requirement decreases from ON to O1. We conclude that that this fast method significantly reduces the total computational cost from OMN2 to OMN and the memory requirement from OMN to OM.


## 6. Numerical Experiments

In this section, we carry out several numerical experiments to check the effectiveness of the proposed scheme. The convergence rate and CPU consumption are all compared in the simulations. We take the space-time domainΩ=0,1,T=1 for one-dimensional case and Ω=0,1×0,1,T=1 for two-dimensional case. These simulations are implemented in Matlab, and the numerical experiments are run on a computer with 4 GB memory. The time-fractional Cattaneo equation of the following forms is considered.Example 1.
We provide the exact solutionux,t=etsinπx, and for different α, we have different fx,t accordingly:(73)∂ux,t∂t+∂αux,t∂tα=∂2ux,t∂x2+fx,t,ux,0=sinπx,∂u∂tt=0=sinπx,u0,t=u1,t=0.
In Tables1 and 2, we take Δt=h2 and h=Δt to examine the discrete l∞-norm (l2-norm) errors and corresponding spatial and temporal convergence rates, respectively. We list the errors and convergence rates (order) of the proposed compact finite difference (CD) scheme, which is almost OΔt2+h4 for different α. Additionally, Table 3 shows the CPU time (CPU) consumed by direct compact (DCD) scheme and fast compact difference (FCD) scheme, respectively. It is obvious that the FCD scheme has a significantly reduced CPU time over the DCD scheme. For instance, when α=1.5, we choose h=0.1 and Δt=1/50,000 and observe that the FCD scheme consumes only 94 seconds, while the DCD scheme consumes 3692 seconds. We can find that the performance of the FCD scheme will be more conspicuous as the time step size Δt decreases.
In Figure1, we set h=0.1 and α=1.75 and change the total number of time steps N to plot out the CPU time (in seconds) of the FCD scheme and DCD scheme. We can observe that the CPU time increases almost linearly with respect to N for the FCD scheme, while the DCD scheme scales like ON2.
Tables4 and 5 show the discrete l∞ errors and convergence rates of the compact finite difference scheme for Example 1. The space rates are almost Oh4 for fixed Δt=2−13, and the time convergence rates are always OΔt2 for fixed h=0.001. We can conclude that the numerical convergence rates of our scheme approach almost to OΔt2+h4.Table 1
ConsideringΔt=h2, the discrete l∞ error and convergence rates of u with different α for Example 1.
h
α
=
1.25
α
=
1.5
α
=
1.75
e
N
∞
Order
e
N
∞
Order
e
N
∞
Order
2
−
3
1.5320e − 04
—
1.5941e − 04
—
1.8076e − 04
—
2
−
4
9.5109e − 06
4.0097
9.8981e − 06
4.0094
1.1230e − 05
4.0086
2
−
5
5.9343e − 07
4.0024
6.1762e − 07
4.0024
7.0079e − 07
4.0022
2
−
6
3.7074e − 08
4.0006
3.8585e − 08
4.0006
4.3783e − 08
4.0005Table 2
Consideringh=Δt, the discrete l2 error and convergence rates of u with different α for Example 1.
Δ
t
α
=
1.25
α
=
1.5
α
=
1.75
e
N
2
Order
e
N
2
Order
e
N
2
Order
2
−
4
0.0017
—
0.0018
—
0.0021
—
2
−
6
1.0832e − 04
4.0382
1.1271e − 04
4.0373
1.2781e − 04
4.0344
2
−
8
6.7252e − 06
4.0097
6.9990e − 06
4.0094
7.9405e − 06
4.0086
2
−
10
4.1961e − 07
4.0024
4.3672e − 07
4.0023
4.9553e − 07
4.0021
2
−
12
2.6215e − 08
4.0006
2.7284e − 08
4.0005
3.0959e − 08
4.0005Table 3
The CPU time consumption of the fast compact difference scheme and direct difference scheme for Example 1.
Δ
t
α
=
1.25
α
=
1.5
α
=
1.75
DCD
FCD
DCD
FCD
DCD
FCD
1
/
1000
3.4008
1.8876
3.4632
1.9032
3.3540
1.9188
1
/
2500
14.2740
4.6800
14.1960
4.7580
13.8840
4.7424
1
/
5000
47.0967
9.3600
46.3010
9.4692
46.1450
9.3756
1
/
7500
99.4194
13.9464
95.6754
14.0088
96.7206
13.8684
1
/
10,000
169.9318
18.6733
164.1130
18.6421
165.8290
18.4393
1
/
25,000
981.6362
46.6910
950.9820
46.9251
959.1409
46.0046
1
/
50,000
3826.9853
93.4133
3692.3096
94.0062
3754.7412
91.9625Figure 1
The CPU time (a) and the log-log CPU time (b) versus the total number of time stepsN for Example 1
(a)
(b)Table 4
ConsideringΔt=2−13, the discrete l∞ error and convergence rates of u with different α for Example 1.
h
α
=
1.25
α
=
1.5
α
=
1.75
e
N
∞
Order
e
N
∞
Order
e
N
∞
Order
2
−
3
2.2486e − 04
—
2.2733e − 04
—
2.3705e − 04
—
2
−
4
1.3985e − 05
4.0071
1.4140e − 05
4.0069
1.4745e − 05
4.0069
2
−
5
8.6897e − 07
4.0084
8.7882e − 07
4.0081
9.1726e − 07
4.0068
2
−
6
5.0194e − 08
4.1137
5.1024e − 08
4.1063
5.4092e − 08
4.0838Table 5
Consideringh=0.001, the discrete l∞ error and convergence rates of u with different α for Example 1.
Δ
t
α
=
1.25
α
=
1.5
α
=
1.75
e
N
∞
Order
e
N
∞
Order
e
N
∞
Order
2
−
3
0.0046
—
0.0043
—
0.0036
—
2
−
4
0.0011
2.0641
0.0011
1.9668
9.0063e − 04
1.9990
2
−
5
2.8662e − 04
1.9403
2.7171e − 04
2.0174
2.2518e − 04
1.9999
2
−
6
7.1658e − 05
1.9999
6.7928e − 05
2.0000
5.6295e − 05
2.0000
2
−
7
1.7915e − 05
2.0000
1.6982e − 05
2.0000
1.4074e − 05
2.0000Example 2.
The example is described by(74)∂ux,t∂t+∂αux,t∂tα=∂2ux,t∂x2+fx,t,ux,0=0,∂u∂tt=0=0,u0,t=t2exp−x02γ,u1,t=t2exp−1−x02γ.
Note that the exact solution of the above problem isux,t=t2exp−x−x02/γ, where x0∈0,1 and γ is a constant, and for different α, we have different fx,t accordingly.
We apply the fast compact difference scheme to discretize the equation. In Figure2, we set γ=0.001,α=1.5,andM=N=100 and plot exact and numerical solutions at time T=1 for Example 2 with different x0. For x0=0.5,α=1.5,andM=N=100, we also plot exact and numerical solutions with the different γ in Figure 3. In Figure 4, for h=0.1 and α=1.5, we vary the total number of time steps N to plot out the CPU time (in seconds) of the FCD scheme and DCD scheme. The numerical experiments verified our theoretical results. In Table 6, by equating Δt=h2 and fixing x0=0.5, we compute the discrete l∞ error and convergence rates with different fractional derivative orders α and different γ. It shows that the compact finite difference scheme has space accuracy of fourth order and temporal accuracy of second order. We set h=Δt and fix γ=0.01, and the discrete l2 error and convergence rates with different α and x0 are displayed in Table 7. These numerical convergence rates are almost approaching OΔt2+h4 as Example 1.Figure 2
Consideringγ=0.001,α=1.5,M=N=100, real and numerical solutions of u at T=1 for Example 2 with the value of x0 equal to (a) 0, (b) 0.2, (c) 0.4, (d) 0.6, (e) 0.8, and (f) 1.
(a)
(b)
(c)
(d)
(e)
(f)Figure 3
Consideringx0=0.5,α=1.5,M=N=100, real and numerical solutions of u at T=1 for Example 2 with the value of γ equal to (a) 1, (b) 0.1, (c) 0.01, and (d) 0.001.
(a)
(b)
(c)
(d)Figure 4
Consideringx0=0.5,γ=0.01,α=1.5,andM=10, the CPU time (a) and the log-log CPU time (b) versus the total number of time steps N for Example 2.
(a)
(b)Table 6
ConsideringΔt=h2andx0=0.5, the discrete l∞ error and convergence rates of u with different α and γ for Example 2.
γ
h
α
=
1.25
α
=
1.5
α
=
1.75
e
N
∞
Order
e
N
∞
Order
e
N
∞
Order
1
2
−
3
7.3224e − 05
—
5.5643e − 05
—
3.8478e − 05
—
2
−
4
4.5707e-06
4.0018
3.4695e − 06
4.0034
2.3900e − 06
4.0090
2
−
5
2.8558e − 07
4.0004
2.1671e − 07
4.0009
1.4913e − 07
4.0024
2
−
6
1.7847e − 08
4.0001
1.3543e − 08
4.0001
9.3169e − 09
4.0006
0.1
2
−
3
0.0012
—
0.0012
—
0.0012
—
2
−
4
7.0462e − 05
4.0900
6.9731e − 05
4.1051
6.9007e − 05
4.1201
2
−
5
4.3571e − 06
4.0154
4.3110e − 06
4.0157
4.2650e − 06
4.0161
2
−
6
2.7160e − 07
4.0038
2.6871e − 07
4.0039
2.6583e − 07
4.0040
0.01
2
−
4
0.0090
—
0.0090
—
0.0090
—
2
−
5
4.9257e − 04
4.1915
4.9235e − 04
4.1922
4.9232e − 04
4.1923
2
−
6
2.9868e − 05
4.0437
2.9855e − 05
4.0436
2.9854e − 05
4.0436
2
−
7
1.8531e − 06
4.0106
1.8523e − 06
4.0106
1.8522e − 06
4.0106Table 7
Consideringh=Δtandγ=0.01, the discrete l2 error and convergence rates of u with different αandx0 for Example 2.
x
0
Δ
t
α
=
1.25
α
=
1.5
α
=
1.75
e
N
2
Order
e
N
2
Order
e
N
2
Order
0
2
−
8
0.0044
—
0.0044
—
0.0042
—
2
−
10
2.4205e − 04
4.1841
2.3721e − 04
4.2133
2.3098e − 04
4.1845
2
−
12
1.4684e − 05
4.0430
1.4390e − 05
4.0430
1.4013e − 05
4.0429
0.25
2
−
8
0.0027
—
0.0027
—
0.0027
—
2
−
10
1.4902e − 04
4.1794
1.4894e − 04
4.1802
1.4892e − 04
4.1803
2
−
12
9.0680e − 06
4.0386
9.0633e − 06
4.0386
9.0616e − 06
4.0386
0.5
2
−
8
0.0027
—
0.0027
—
0.0027
—
2
−
10
1.4848e − 04
4.1846
1.4842e − 04
4.1852
1.4842e − 04
4.1852
2
−
12
9.0323e − 06
4.0390
9.0288e − 06
4.0390
9.0286e − 06
4.0390
0.75
2
−
8
0.0027
—
0.0027
—
0.0027
—
2
−
10
1.4902e − 04
4.1794
1.4894e − 04
4.1802
1.4892e − 04
4.1803
2
−
12
9.0680e − 06
4.0386
9.0633e − 06
4.0386
9.0616e − 06
4.0386
1
2
−
8
0.0044
—
0.0044
—
0.0042
—
2
−
10
2.4205e − 04
4.1841
2.3721e − 04
4.2133
2.3098e − 04
4.1845
2
−
12
1.4684e − 05
4.0430
1.4390e − 05
4.0430
1.4013e − 05
4.0429Example 3.
If the exact solution is given byux,t=etsinπxsinπy, we have different fx,y,t for different α accordingly:(75)∂ux,y,t∂t+∂αux,y,t∂tα=∂2ux,y,t∂x2+∂2ux,y,t∂y2+fx,y,t,ux,0=sinπxsinπy,∂u∂tt=0=sinπxsinπy,u0,y,t=u1,y,t=ux,0,t=ux,1,t=0.
In Figure5, h=0.1 and α=1.75 are fixed, and the total number of time steps N vary to plot out the CPU time (in seconds) of the FCD procedure and DCD procedure, and it presents an approximately linear computation complexity for FCD procedure. We set Δt=h2 in Table 8, and h=Δt in Table 9, the discrete l∞ error, discrete l2 error, and convergence rates with different derivative orders α are presented. The fourth-order space accuracy and second-order temporal accuracy can be observed clearly.Figure 5
The CPU time (a) and the log-log CPU time (b) versus the total number of time stepsN for Example 3.
(a)
(b)Table 8
ConsideringΔt=h2, the discrete l∞ error and convergence rates of u with different α for Example 3.
h
α
=
1.25
α
=
1.5
α
=
1.75
e
N
∞
Order
e
N
∞
Order
e
N
∞
Order
2
−
2
0.0028
—
0.0028
—
0.0030
—
2
−
3
1.6948e − 04
4.0462
1.7328e − 04
4.0142
1.8411e − 04
4.0263
2
−
4
1.0522e − 05
4.0096
1.0759e − 05
4.0095
1.1434e − 05
4.0092
2
−
5
6.5652e − 07
4.0024
6.7133e − 07
4.0024
7.1346e − 07
4.0024Table 9
Consideringh=Δt, the discrete l2 error and convergence rates of u with different α for Example 3.
Δ
t
α
=
1.25
α
=
1.5
α
=
1.75
e
N
2
Order
e
N
2
Order
e
N
2
Order
2
−
4
0.0014
—
0.0014
—
0.0015
—
2
−
6
8.4738e − 05
4.0463
8.6640e − 05
4.0142
9.2056e − 05
4.0263
2
−
8
5.2610e − 06
4.0096
5.3795e − 06
4.0095
5.7168e − 06
4.0092
2
−
10
3.2826e − 07
4.0024
3.3566e − 07
4.0024
3.5673e − 07
4.0023


## 7. Conclusion

In this paper, we develop and analyze a fast compact finite difference procedure for the Cattaneo equation equipped with time-fractional derivative without singular kernel. The time-fractional derivative is of Caputo–Fabrizio type with the order ofα1<α<2. Compact difference discretization is applied to obtain a high-order approximation for spatial derivatives of integer order in the partial differential equation, and the Caputo–Fabrizio fractional derivative is discretized by means of Crank–Nicolson approximation. It has been proved that the proposed compact finite difference scheme has spatial accuracy of fourth order and temporal accuracy of second order. Since the fractional derivatives are history dependent and nonlocal, huge memory for storage and computational cost are required. This means extreme difficulty especially for a long-time simulation. Enlightened by the treatment for Caputo fractional derivative [32], we develop an effective fast evaluation procedure for the new Caputo–Fabrizio fractional derivative for the compact finite difference scheme. Several numerical experiments have been carried out to show the convergence orders and applicability of the scheme.Inspired by the work [43], the topic about modelling and numerical solutions of porous media flow equipped with fractional derivatives is very interesting and challenging and will be our main research direction in the future.



---

*Source: 3842946-2020-03-19.xml*