A compact exponential difference method for multi-term time-fractional convection-reaction-diffusion problems with non-smooth solutions

https://doi.org/10.1016/j.amc.2020.125316Get rights and content

Highlights

  • Multi-term time-fractional convection-reaction-diffusion problems with non-smooth solutions are considered.

  • A compact exponential finite difference method is proposed

  • Taking into account the initial weak singularity of the solution, the stability and convergence of the method is rigorously proved.

  • The spatial fourth-order convergence and the temporal optimal convergence are obtained.

  • Numerical results confirm the theoretical convergence result.

Abstract

This paper is concerned with a numerical method for a class of one-dimensional multi-term time-fractional convection-reaction-diffusion problems, where the differential equation contains a sum of the Caputo time-fractional derivatives of different orders between 0 and 1. In general the solutions of such problems typically exhibit a weak singularity at the initial time. A compact exponential finite difference method, using the well-known L1 formula for each time-fractional derivative and a fourth-order compact exponential difference approximation for the spatial discretization, is proposed on a mesh that is generally nonuniform in time and uniform in space. Taking into account the initial weak singularity of the solution, the stability and convergence of the method is proved and the optimal error estimate in the discrete L2-norm is obtained by developing a discrete energy analysis technique which enables us to overcome the difficulties caused by the nonsymmetric discretization matrices. The error estimate shows that the method has the spatial fourth-order convergence, and reveals how to select an appropriate mesh parameter to obtain the temporal optimal convergence. The extension of the method to two-dimensional problems is also discussed. Numerical results confirm the theoretical convergence result, and show the applicability of the method to convection dominated problems.

Introduction

Fractional differential problems have been successfully applied to the modelling of many different processes and systems, see [1], [2], [3], [4]. Among these problems, multi-term fractional differential problems, where the differential equation contains multiple fractional derivatives, provide a more flexible and precise tool for describing some complex physical phenomena in practice, such as the oxygen delivery through a capillary to tissues [5], the total concentration in solute transport [6], and the sub-diffusive motion in velocity fields [7]. A great deal of work has been done in qualitative analysis of multi-term fractional differential problems. Luchko [8] considered a class of the generalized multi-term time fractional diffusion equations and established some a priori solution estimates. Li et al. [9] investigated the well-posedness and the long-time asymptotic behaviour for initial-boundary value problems of multi-term time-fractional diffusion equations. Using the method of separation of variables, Jiang et al. [10] derived analytic solutions of the generalized multi-term time-fractional diffusion-wave/diffusion equations with nonhomogeneous Dirichlet, Neumann and Robin boundary conditions. By using Luchko’s Theorem and the equivalent relationship between the Laplacian operator and the Riesz fractional derivative, Jiang et al. [11] proposed some new analytic techniques to solve three types of multi-term time-space Caputo-Riesz fractional advection diffusion equations with nonhomogeneous Dirichlet boundary conditions. Some numerical algorithms for multi-term fractional differential problems have also been developed. Shiralashetti et al. [12] proposed a Haar wavelet collocation method for solving multi-term fractional differential equations. Dehghan et al. [13] developed a high-order finite difference method and the Galerkin spectral method for the numerical solutions of multi-term time-fractional differential equations. Liu et al. [14] discussed several finite difference schemes for multi-term time fractional advection-dispersion and wave-diffusion equations based on a fractional predictor-corrector method. Wang [15] and Wang et al. [16] presented a high-order compact finite difference method for solving the two-term fractional mobile/immobile convection-diffusion equation and the two-term modified anomalous fractional sub-diffusion/diffusion-wave equation, respectively. Using a fractional predictor-corrector method combining the L1 and L2 formulae, Ye et al. [17] studied a numerical method for a multi-term time-space Riesz-Caputo fractional differential equation. Recently, Xu et al. [18] constructed a finite difference scheme for a multi-term variable-order fractional diffusion equation which involves the Caputo variable-order time-fractional derivative and the Riesz variable-order space-fractional derivative, and Alikhanov [19] investigated the finite difference methods for the Dirichlet and Robin boundary value problems of a multi-term variable-distributed order diffusion equation. We also mention that Wu et al. [20], [21] introduced a concept of short memory and proposed a short memory fractional differential equation where the initial point in the Caputo fractional derivative is varied. Since the fractional order is defined by the use of a piecewise constant function on different time intervals, the short memory fractional differential equation contains multiple fractional derivatives. This kind of new variable-order fractional differential equation holds short memory effects and provides a new tool for fractional modelling (see [20], [21]).

We propose and analyze a simple and yet efficient numerical method for the following multi-term time-fractional convection-reaction-diffusion problem:{r=1mρr0CDtαru(x,t)=Lu(x,t)+c(x,t)u(x,t)+f(x,t),(x,t)D(0,L)×(0,T],u(0,t)=ϕ0(t),u(L,t)=ϕL(t),t(0,T],u(x,0)=ψ(x),x[0,L],where ρr (1 ≤ r ≤ m) are some positive constants, 0<αm<αm1<<α1<1, 0CDtαu(x,t) represents the αth-order Caputo fractional derivative in t, which is defined by0CDtαu(x,t)=1Γ(1α)0tsu(x,s)(ts)αds,0<α<1,and the term Lu(x,t) is the convection-diffusion term given byLu(x,t)=ax2u(x,t)pxu(x,t)with the convection coefficient p0 and the diffusion coefficient a > 0. Without loss of generality we take ρ1=1. Using the technique of [22], one can generalize Theorem 2.1 of [23] (also see [24]) to the present multi-term time-fractional convection-reaction-diffusion problem (1.1) and conclude that, under suitable regularity and compatibility assumptions, the problem (1.1) has a unique solution u in C([0,L]×[0,T]) such that(i)|x6u(x,t)|1,(ii)|tlu(x,t)|1+tα1l,(x,t)[0,L]×(0,T];l=1,2,where the notation  ≲  is rigorously defined in the final paragraph of this section. The bound (1.4)-(ii) indicates that the solution u exhibits a weak singularity at the initial time t=0.

The problem (1.1) gives a general type of one-dimensional multi-term time-fractional convection-reaction-diffusion problem. A classical example for an application of it is the so-called Basset equation which arises in the modelling of the forces that occur when a spherical object sinks in a (relatively less dense) incompressible viscous fluid (see, e.g., [25], [26]). Another example is the fractional mobile/immobile convection-diffusion model for the total concentration in solute transport. This model was originally proposed in Schumer et al. [6] and was thoroughly discussed, e.g., in [27].

A lot of work has been devoted to numerical methods for the multi-term time-fractional problem (1.1) with p=0 and c(x, t) ≡ 0, for example, a high-order space-time spectral method in Zheng et al. [28], the Galerkin finite element method in Jin et al. [29], the temporal second-order difference methods based on a super-convergence interpolation approximation in Gao et al. [30], the finite difference method based on the L1 formula in Ren and Sun [31], the finite element method combing with the Diethelm fractional backward difference method in Zhao et al. [32], and a fully-discrete method using the finite element method in spatial direction and the L1 formula in temporal direction in Zhao et al. [33]. However, the convergence analysis in most of these works was done under the assumption that the solution u is at least twice continuously differentiable in time t on the closed domain [0, L] × [0, T] and therefore it excludes the case of non-smooth solutions with the initial weak singularity (1.4)-(ii). In the very recent work [23], Huang et al. proposed and analyzed a finite difference method for the problem (1.1) in multidimensional space with p=0 and c(x, t) ≡ c(x) ≤ 0 — a multi-term time-fractional reaction-diffusion problem. A sharp convergence analysis of the method was given, taking into account the initial weak singularity (1.4)-(ii). In that method, each time-fractional derivative is approximated by the L1 formula on a temporal graded mesh, while the spatial derivative is discretized by the classical second-order central difference approximation on a spatial uniform mesh. The same multi-term time-fractional reaction-diffusion problem was considered numerically in Huang and Stynes [24], where a standard finite element method was used for the spatial discretization.

In this paper, we continue the work of Huang et al. [23] by studying a numerical method with high-order spatial accuracy for the multi-term time-fractional convection-reaction-diffusion problem (1.1). We consider the general case of the problem (1.1), where p0, the reaction coefficient c ≡ c(x, t) is variable both spatially and temporally, and c(x, t) is not required to be non-positive. We propose and analyze a high-order compact finite difference method while taking into account the initial weak singularity (1.4)-(ii). The proposed method uses a similar compact stencil in space as the classical second-order central difference method, and so its calculation is not more complicated than the method used in Huang et al. [23], but it has higher spatial accuracy. In order to deal with the initial weak singularity (1.4)-(ii) better, we use a more general nonuniform time mesh with a weak assumption, which includes the temporal graded mesh used in Huang et al. [23], Huang and Stynes [24] as a special case. This makes the proposed method more flexible and effective in practical application.

There are many high-order compact finite difference methods in the literature for solving convection-reaction-diffusion problems. We here use the high-order compact exponential finite difference method, since it usually results in excellent approximations for convection dominated problems (see [34], [35]). Due to the convection term, the coefficient matrices of the proposed method are not symmetric. This fact considerably complicates the convergence analysis of the proposed method, and the analysis technique used in Huang et al. [23] is no longer valid. In order to tackle this problem, we formulate the proposed method as a suitable matrix form by decomposing the coefficient matrices carefully, and then apply the discrete energy analysis method to the matrix form. This analysis technique enables us to overcome the difficulties caused by the nonsymmetric discretization matrices. It is also very different from the one used in Cui [34], where a compact exponential finite difference method was applied to solve a single-term time-fractional convection-reaction-diffusion problem, but the initial weak singularity of the solution was not considered.

The outline of the paper is as follows. In Section 2, we propose a high-order compact exponential finite difference method for the multi-term time-fractional problem (1.1), based on the L1 formula on a general nonuniform time mesh. A full analysis of the stability and convergence of the method is carried out in Section 3 under the regularity assumption (1.4). Section 4 is devoted to discussing the extension of the proposed method to two-dimensional problems. In Section 5, we present some numerical results to confirm the theoretical convergence result and show the applicability of the method to convection dominated problems. We also make some numerical comparisons of the method with the standard central difference method given in Huang et al. [23]. The concluding section contains a brief conclusion.

Notation. Throughout the paper, the notation C, with or without subscript, denotes a positive constant which may differ at different occurrences, but it is always independent of the spatial mesh size, the time step size and the time levels. We use A ≲ B to mean A ≤ CB, and A ≃ B to mean that A ≲ B and B ≲ A. When we write, e.g., tntn1 for all n, we mean that there is a single fixed C such that tnCtn1 for all n.

Section snippets

Compact exponential difference method

For a given positive integer M, we set h=L/M and partition [0, L] into a mesh by the mesh points xi=ih (0 ≤ i ≤ M). Let N be a positive integer. In order to deal with the initial weak singularity (1.4)-(ii) and obtain an effective numerical method for the problem (1.1), we consider a general nonuniform time mesh 0=t0<t1<<tN=T with the assumption: there are two constants γ ≥ 1 and ρ > 0 such that{(i)τn:=tntn1N1tn11γ(1nN),tntn1(2nN),(ii)τn/τn+1ρ(1nN1).

The assumption (2.1)-(i)

Stability and convergence

We carry out the stability and convergence analysis of the compact exponential difference scheme (2.22) under the regularity assumption (1.4) and the mesh assumption (2.1).

Let Sh={v|v=(v1,v2,,vM1)}. For grid functions v,wSh, we define(v,w)=hi=1M1viwi,v=(v,v)12.

For any vSh, we define the vector vRM1 by v=(v1,v2,,vM1)T. It is clear that for any v,wSh, (v,w)=hvTw and v2=hvTv.

Lemma 3.1

The matrix P1 is symmetric and positive definite. Moreover, λmin(P1)23 and λmax(P1) ≤ 1, where λmin(P1) and λ

Extension to two-dimensional problems

In this section, we extend the compact exponential finite difference scheme (2.22) to two-dimensional problems. Let Ω=(0,Lx)×(0,Ly) and Ω¯=[0,Lx]×[0,Ly]. We consider the following two-dimensional form of the problem (1.1):{r=1mρr0CDtαru(x,y,t)=(Lx+Ly)u(x,y,t)+c(x,y,t)u(x,y,t)+f(x,y,t),(x,y,t)DΩ×(0,T],u(x,y,t)=ϕ(x,y,t),(x,y,t)Ω×(0,T],u(x,y,0)=ψ(x,y),(x,y)Ω¯,where ∂Ω is the boundary of Ω andLxu(x,y,t)=ax2u(x,y,t)pxu(x,y,t),Lyu(x,y,t)=by2u(x,y,t)qyu(x,y,t)with a > 0, b > 0, p0 and q0.

Numerical results

We first present numerical results for three test examples of the form (1.1). The exact solution u of each example is explicitly known and is used to compare with the numerical solution un=(u1n,u2n,,uM1n) of the compact exponential difference scheme (2.22) to confirm the theoretical convergence result of Theorem 3.3 and show the applicability of the scheme (2.22) to convection dominated problems. We also make some numerical comparisons of the compact scheme (2.22) with the standard central

Concluding remarks

We have given a numerical treatment for a class of one-dimensional multi-term time-fractional convection-reaction-diffusion problems, where the differential equation contains a sum of the Caputo time-fractional derivatives of different orders between 0 and 1. In general the solutions of such problems typically exhibit a weak singularity at the initial time. However, this initial weak singularity was ignored in most existing numerical works, where the convergence analysis was done under the

Acknowledgements

The authors would like to thank the referees for their valuable comments and suggestions which improved the presentation of the paper.

References (41)

  • J.C. Ren et al.

    Efficient and stable numerical methods for multi-term time fractional sub-diffusion equations

    East Asian J. Appl. Math.

    (2014)
  • J. Zhao et al.

    Stability and convergence of an effective finite element method for multiterm fractional partial differential equations

    Abstr. Appl. Anal.

    (2013)
  • H.-L. Liao et al.

    Sharp error estimate of the nonuniform L1 formula for linear reaction-subdiffusion equations

    SIAM J. Numer. Anal.

    (2018)
  • N. Kopteva

    Error analysis of the L1 method on graded and uniform meshes for a fractional-derivative problem in two and three dimensions

    Math. Comput.

    (2019)
  • H. Brunner

    Collocation Methods for Volterra Integral and Related Functional Differential Equations

    (2004)
  • I. Podlubny

    Fractional Differential Equations

    (1999)
  • K.B. Oldham et al.

    The Fractional Calculus

    (1974)
  • D. Baleanu et al.

    Fractional Calculus: Models and Numerical Methods

    (2012)
  • M.H. Srivastava et al.

    An efficient analytical technique for fractional model of vibration equation

    Appl. Math. Model.

    (2017)
  • R. Schumer et al.

    Fractal mobile/immobile solute transport

    Water Resour. Res.

    (2003)
  • Cited by (9)

    View all citing articles on Scopus

    This work was supported in part by Science and Technology Commission of Shanghai Municipality (STCSM) (No. 18dz2271000).

    View full text