Implicit-Explicit difference schemes for nonlinear fractional differential equations with non-smooth solutions

# Implicit-Explicit difference schemes for nonlinear fractional differential equations with non-smooth solutions

Wanrong Cao111Department of Mathematics, Southeast University, Nanjing 210096, P.R.China. (wrcao@seu.edu.cn), Fanhai Zeng222 Division of Applied Mathematics, Brown University, Providence RI, 02912, USA. (fanhai_zeng@brown.edu), Zhongqiang Zhang333Department of Mathematical Sciences, Worcester Polytechnic Institute, Worcester, MA, 01609, USA. (zzhang7@wpi.edu). , and George Em Karniadakis 444 Division of Applied Mathematics, Brown University, Providence RI, 02912, USA. (george_karniadakis@brown.edu).
###### Abstract

We propose second-order implicit-explicit (IMEX) time-stepping schemes for nonlinear fractional differential equations with fractional order . From the known structure of the non-smooth solution and by introducing corresponding correction terms, we can obtain uniformly second-order accuracy from these schemes. We prove the convergence and linear stability of the proposed schemes. Numerical examples illustrate the flexibility and efficiency of the IMEX schemes and show that they are effective for nonlinear and multi-rate fractional differential systems as well as multi-term fractional differential systems with non-smooth solutions.

Keywords time-fractional derivatives, IMEX schemes, low regularity, multi-rate systems, multi-term fractional derivatives

## 1 Introduction

We aim at constructing efficient finite difference schemes for fractional ordinary differential equations (FODEs) with non-smooth solutions. In recent decades, due to the increasing interest in problems with anomalous transport dynamics, fractional differential equations have become significant mathematical models in many fields of science and engineering, such as viscoelastic models in blood flow [33], underground transport [20], options pricing model in financial markets [41], etc.

Though some fractional differential equations (FDEs) with special form, e.g., linear equations, can be solved by analytical methods, e.g., the Fourier transform method or the Laplace transform method [34], the analytical solutions of many generalized FDEs (e.g. nonlinear FDEs and multi-term FDEs) are rather difficult to obtain. This encourages us to develop effective numerical methods for solutions of these FDEs. Up to now, a number of finite difference methods have been established for FDEs. One way is to transform the considered FDEs into their integral forms, then numerical methods for the fractional integral operator are developed and the corresponding difference schemes are derived; see [4, 5, 10, 11, 17, 18, 15, 29, 30, 37, 48]. Another approach is to approximate the fractional derivative operators in the considered FDEs directly; see [16, 27, 36, 38, 54]. Besides finite difference methods, there exist also other numerical methods for FDEs, e.g., finite element methods [22, 40], spectral methods [26, 39, 45], matrix methods [32, 35], etc.

For nonlinear FDEs, most of the aforementioned finite difference methods are implicit, and the nonlinear system needs to be solved using the iteration method, which is costly. To avoid extra computational costs and instability caused by iteration, various numerical methods have been proposed, e.g., predictor-corrector methods; see [6, 7, 10, 11, 17, 25, 42, 53] and time-splitting schemes [1].

Implicit-explicit (IMEX) schemes (also known as semi-implicit schemes, linearly-implicit schemes) play a major rule in the numerical treatment of nonlinear/stiff differential equations [12, 21, 24]. Due to approximating the nonlinear/nonstiff part explicitly, we avoid solving nonlinear equations/systems in every time step, and hence the computational cost can be reduced significantly. Moreover, the IMEX schemes generally have better stability than explicit schemes, for the linear/stiff parts are treated implicitly. To the best of authors’ knowledge, there are very limited works on IMEX schemes for FDEs. An IMEX scheme has been proposed for nonlinear fractional anomalous diffusion equations with smooth solutions in [47] and an explicit and implicit finite difference scheme for fractional Cattaneo equation was developed in [19]. Some semi-implicit methods for space-fractional differential equations can be found in the literature, e.g., see [3, 28, 43].

The main contribution of this work is to develop second-order IMEX schemes for nonlinear/stiff FODEs with solutions that have a weak singularity at the origin. Most numerical methods for differential equations are generally intended for problems with solutions of high regularity. However, solutions to FDEs usually have weak singularities at the origin even when the forcing term is smooth, see e.g. [9, 11, 23]. When solving these FDEs, the singularity requires special attention to obtain the expected high accuracy. Several approaches have been proposed to deal with the weak singularity, such as using adaptive grids (nonuniform grids) to keep errors small near the singularity [18, 44, 51], or employing non-polynomial basis functions to include the correct singularity index [2, 14, 52], or using the correction terms to remedy the loss of accuracy and recover high-order schemes [18, 30, 46, 49, 50].

In this paper, we follow the last approach to develop IMEX methods with uniformly second-order accuracy for FODEs, whose solutions are non-smooth and have known structure. To deal with the singularity, we apply correction terms in the proposed schemes, so that the resulting fractional quadratures are either exact or sufficiently accurate for the weakly singular parts of the solutions. The idea of adding correction terms was firstly proposed for approximating fractional-order integrals by the linear multi-step method in [30]. Very recently, the same strategy has been adopted to enhance the accuracy of numerical schemes for FDEs [49, 50]. To treat the nonlinear part explicitly, we further use extrapolation and Taylor expansion to approximate the nonlinear part, where appropriate correction terms are also used to obtain high accuracy for non-smooth solutions. We propose two IMEX schemes, which can work for nonlinear/stiff FODEs with uniformly second-order convergence, even when the solutions have weak singularity at the origin. We also prove the convergence and linear stability of the proposed schemes.

We organize this work as follows. In Section 2, we formulate the IMEX method based on the fractional linear multistep method (FLMM), and then derive two IMEX schemes by applying the extrapolation and Taylor expansion for the nonlinear terms. Moreover, we provide the strategies for introducing correction terms to the schemes to make them uniformly second-order for FODEs with non-smooth solutions. We also present convergence rates of the proposed schemes, the proofs of which are given in Section 6. In Section 3, we discuss the linear stability of the proposed schemes. We present more details on the proposed IMEX schemes in Section 4 to show that they can be applied to stiff systems and multi-term nonlinear FODEs. In Section 5, we present numerical examples to illustrate the computational flexibility and verify our error estimates. We conclude in Section 7 and discuss the performance of the proposed schemes.

## 2 Second-order IMEX schemes

We consider the following nonlinear FODE

 (CDβ0u)(t) = λu(t)+fu(t),t∈(0,T],u(0)=u0, (2.1)

where , , is the Caputo derivative defined by

 (CDβag)(t)=(I1−βag′)(t),(Iβag)(t)=1Γ(β)∫tag(τ)(t−τ)1−βdτ, t>a. (2.2)

We first transform (2.1) into its integral form as

 u(t) = u(0)+λ(Iβ0u)(t)+(Iβ0fu)(t),0≤t≤T. (2.3)

Eq. (2.3) is readily obtained by applying the operator on both sides of (2.1) and using the identity , see, e.g. [11].

In order to obtain second-order schemes, we need to approximate the fractional integrals in (2.3) with second-order quadrature rules. However, the solutions to (2.1) usually have singularity at . The analytical solution to (2.1) can be written as the summation of regular and singular parts, as given in the following lemma. See also [8, 9, 11] for more discussions.

###### Lemma 2.1 ([11])

Suppose , where is a suitable region of variable . Then there exists a function and some and , such that the solution of (2.1) is of the form

 u(t)=ψ(t)+^ν∑ν=1cνtνβ+^η∑η=1dηt1+ηβ, (2.4)

where , .

The solution of (2.1) is usually non-smooth, even if is smooth. Consequently, many existing numerical methods (see e.g.[10, 18, 47]) for (2.1) would produce less accurate numerical solutions when they are directly applied. Next, we will adopt a second-order FLMM developed in [30] to (2.3) to construct our IMEX schemes. Take a uniform partition of time interval , i.e., with . The second-order FLMM used in the present work reads

 Iβ0u(tn)=hβn∑j=0ω(β)n−ju(tj)+hβm∑j=0W(β)n,ju(tj)+O(h2), (2.5)

where are coefficients of the Taylor expansion of the following generating function

 ω(β)(z)=(121+z1−z)β=∞∑j=0ω(β)jzj, (2.6)

and are the starting weights that recover second-order accuracy. If we drop the correction terms in (2.5), then we would lose the second-order accuracy, unless satisfies some special conditions, i.e., is smooth and .

The following lemma states the convergence of (2.5) with no correction terms.

###### Lemma 2.2 ([30, 48])

If , then for (2.5) with ,

 (Iβ0u)(tn)=hβn∑k=0ω(β)n−ku(tk)+O(h2tν+β−2n)+O(h1+νtβ−1n), (2.7)

where are defined by (2.6).

From Lemma 2.1, we see that the analytical solution of FODE (2.1) has the form of (2.4) when satisfies some suitable conditions. In the following, we will always assume that can be expressed in the following form for convenience,

 u(t)−u(0)=m+1∑r=1crtσr+ξ(t)tσm+2,0<σr<σr+1, (2.8)

where is a uniformly continuous function over the interval and are constants. The sequence is uniquely determined by the considered equation. For example, when in (2.1) is smooth, ’s are of the form , see [34, Chapter 5]. Another example is from Lemma 2.1, ’s are of the form (2.4).

Given a sequence of positive numbers , we define the operator as

 Iβ,n,mh,θg=hβn∑k=0ω(β)n−kg(tk)+hβm∑k=1W(β,θ)n,kg(tk)+hβBθng(t0), (2.9)

where satisfies (2.6), and and are given by

 m∑k=1W(β,θ)n,kkθr = Γ(θr+1)Γ(θr+1+β)nθr+β−n∑k=0ω(β)n−kkθr,1≤r≤m, (2.10) Bθn = nβΓ(1+β)−n∑k=0ω(β)n−k−m∑k=1W(β,θ)n,k. (2.11)

Here in (2.9) are called the starting weights that are chosen such that when , which leads to (2.10).

###### Remark 2.3

The linear system (2.10) is ill-conditioned when is large [8, 30, 50]. The large condition number of the Vandermonde-type matrix in (2.10) may lead to big roundoff errors of the starting weights . However, we do not need many correction terms to get satisfactory numerical solutions in computation as observed in [50]. With this observation, we only need to solve the system (2.10) with moderately large condition number. Thus we can obtain reasonable accuracy of the staring weights and hence the numerical solutions, see [8, 50]. We present residuals of the system (2.10) and its corresponding condition numbers in Example 5.1 and Example 5.2.

If we apply (2.9) to approximate , where satisfies (2.8), then by (2.7) we have

 (Iβ0u)(tn) = (Iβ0(u−u(0)))(tn)+tβnu(0)Γ(1+β) (2.12) = hβn∑k=0ω(β)n−k(u(tk)−u0)+hβmu∑k=1W(β,σ)n,k(u(tk)−u0)+tβnu0Γ(1+β)+Rnu = Iβ,n,muh,σu+Rnu,

where is defined by (2.9) and is defined by

 Rnu = O(h2tσmu+1+β−2n)+O(h1+σmu+1tβ−1n). (2.13)

From (2.1), we have

 fu(t)=f(t,u(t))=(CDβ0u)(t)−λu(t). (2.14)

Hence, the regularity of is related to the regularity of . In fact, based on the smoothness assumption of (see Eq. (2.8)), we obtain that has the form

 f(t,u(t))−f(0,u(0))= −λm+1∑r=1crtσr+m+1∑r=1crΓ(σr+1)Γ(σr+1−β)tσr−β+⋯ (2.15) = l+1∑r=1drtδr+ζ(t)tδl+2,

where is uniformly bounded on , , and with ’s from (2.8).

Similar to (2.12), we have

 (Iβt0fu)(tn) = Iβ,n,mfh,δfu+Rnf, (2.16)

where is defined by (2.9),

and the truncation error is given by

 Rnf = O(h2tδmf+1+β−2n)+O(h1+δmf+1tβ−1n). (2.17)

From (2.12) and (2.16), we can derive the following implicit discretization for (2.3)

 u(tn)=u0+λIβ,n,muh,σu+Iβ,n,mfh,δfu+Rnu+Rnf, (2.18)

where and are defined by (2.9), and are defined by (2.13) and (2.17), respectively.

Let be the approximate solution of . Dropping the truncation errors and in (2.18) and replacing with , we derive the fully implicit method for (2.3): to find for such that

 Un= u0+λIβ,n,muh,σU+Iβ,n,mfh,δF, (2.19)

where , .

Given , we need to solve a nonlinear system (2.19) at each time step to get . Next, we further use the extrapolation and Taylor expansion with correction terms to approximate in (2.18), which leads to linear systems and also preserves high-order accuracy.

If , then we have from the Taylor expansion that

 u(tn)=2u(tn−1)−u(tn−2)+O(h2tσ−2n),n≥2, (2.20) u(tn)=u(tn−1)+hu′(tn−1)+O(h2tσ−2n),n≥2. (2.21)

(1) By extrapolation with correction terms: It is clear that (2.20) does not preserve globally second-order accuracy when . Hence, is not a second-order approximation of when is not sufficiently smooth, see (2.15). By adding correction terms to (2.20), we can obtain

 f(tn,u(tn)) = 2f(tn−1,u(tn−1))−f(tn−2,u(tn−2)) (2.22) +˜mf∑k=1ˆW(f)n,k(f(tk,u(tk))−f(t0,u0))+˜Rnf,n≥2,

where are chosen such that the above equation (2.22) is exact, i.e., , for , i.e., satisfy

 ˜mf∑k=1ˆW(f)n,kkδr=nδr−2(n−1)δr+(n−2)δr,r=1,⋯,˜mf. (2.23)

The truncation error in (2.22) satisfies

 ˜Rnf=O(h2tδ˜mf+1−2n) (2.24)

when satisfies (2.15).

Inserting (2.22) into (2.18) yields

 u(tn)= u0+λIβ,n,muh,σu+Iβ,n,mfh,δfu (2.25) +hβω(β)0[−f(tn,u(tn))+2f(tn−1,u(tn−1))−f(tn−2,u(tn−2)) +˜mf∑k=1ˆW(f)n,k(f(tk,u(tk))−f(t0,u0))]+RnE,

where , , , and are defined by (2.13), (2.17), and (2.24), respectively.

From (2.25), we obtain the IMEX method based on the extrapolation technique (abbreviated as IMEX-E) as: given , to find such that

 Un= U0+λIβ,n,muh,σU+Iβ,n,mfh,δF−hβω(β)0Fn (2.26) +hβω(β)0[2Fn−1−Fn−2+˜mf∑k=1ˆW(f)n,k(Fk−F0)],

where , , and are defined by (2.9), and is given by (2.23).

###### Remark 2.4

Given , Eq. (2.26) is a linear equation of . In fact, contains , and it can be eliminated by the following term in the scheme.

(2) By Taylor expansion with correction terms: From the Taylor expansion and (2.21), we have

 f(tn,u(tn))= f(tn−1,u(tn−1))+hf′(tn−1,u(tn−1)) (2.27) +˜mf∑k=1˜W(f)n,k(f(tk,u(tk))−f(t0,u0))+O(h2tδ˜mf+1−2n),n≥2,

where , and the starting weights are chosen such that Eq. (2.27) is exact for , i.e., satisfy

 ˜mf∑k=1˜W(f)n,kkδr=nδr−(n−1)δr−δr(n−1)δr−1,1≤r≤˜mf. (2.28)

Next, we approximate with correction terms, which is given by

 f′(tn−1,u(tn−1)) = ∂uf(tn−1,u(tn−1))∂tu(tn−1)+∂tf(tn−1,u(tn−1)) (2.29) = ∂uf(tn−1,u(tn−1))[u(tn)−u(tn−1)h+1h˜mu∑k=1˜W(u)n,k(u(tk)−u0) +O(htσ˜mu+1−2n)]+∂tf(tn−1,u(tn−1)),n≥2,

where are chosen such that

 u(tn)−u(tn−1)h+1h˜mu∑k=1˜W(u)n,k(u(tk)−u0)=u′(tn−1)

for some , i.e., satisfy

 ˜mu∑k=1˜W(u)n,kkσr=σr(n−1)σr−1−(nσr−(n−1)σr),1≤r≤˜mu. (2.30)

Combining (2.27) and (2.29) yields

 f(tn,u(tn)) = f(tn−1,u(tn−1))+h∂tf(tn−1,u(tn−1)) (2.31) +∂uf(tn−1,u(tn−1))[u(tn)−u(tn−1)+˜mu∑k=1˜W(u)n,k(u(tk)−u0)] +˜mf∑k=1˜W(f)n,k(f(tk,u(tk))−f(t0,u0))+˜Rnu,f,

where the truncation error satisfies

 ˜Rnu,f= O(h2tδ˜mf+1−2n)+O(h2tσ˜mu+1−2n). (2.32)

Inserting (2.31) into (2.18) leads to

 u(tn)= u0+λIβ,n,muh,σu+Iβ,n,mfh,δfu (2.33) +hβω(β)0[−f(tn,u(tn))+f(tn−1,u(tn−1))+h∂tf(tn−1,u(tn−1)) +∂uf(tn−1,u(tn−1))(u(tn)−u(tn−1)+˜mu∑k=1˜W(u)n,k(u(tk)−u0)) +˜mf∑k=1˜W(f)n,k(f(tk,u(tk))−f(t0,u0))]+RnT,

where , , , and are defined by (2.13), (2.17), and (2.32), respectively.

From (2.33), we obtain the IMEX method based on the Taylor expansion technique (abbreviated as IMEX-T) as: given , to find such that

 Un =U0+λIβ,n,muh,σU+Iβ,n,mfh,δF+hβω(β)0[−Fn+Fn−1+h∂tf(tn−1,Un−1) (2.34) +∂uf(tn−1,Un−1)(Un−Un−1+˜mu∑k=1˜W(u)n,k(Uk−U0))+˜mf∑k=1˜W(f)n,k(Fk−F0)],

where , , and are defined by (2.9), and are given by (2.28) and (2.30), respectively.

Next, we present the convergence results for the two schemes (2.26) and (2.34), the proofs of which will be given in Section 6.

###### Theorem 2.5 (Convergence of IMEX-E)

Suppose that is the solution to (2.1) that satisfies (2.8) and is the solution to (2.26), and satisfies the Lipschitz condition with respect to the second argument . If and , then there exists a positive constant independent of and such that

 |u(tn)−Un|≤C(m∑k=1|u(tk)−Uk|+hq), (2.35)

where and

###### Theorem 2.6 (Convergence of IMEX-T)

Suppose that is the solution to (2.1) that satisfies (2.8) and is the solution to (2.34), and satisfies the Lipschitz condition with respect to the second argument . If and , then there exists a positive constant independent of and such that

 |u(tn)−Un|≤C(m∑k=1|u(tk)−Uk|+hq), (2.36)

where and

 q=min{2,σmu+1+β,δmf+1+β,δ˜mf+1+β,σ˜mu+1+β}.

## 3 Linear stability of IMEX schemes

In this section, we discuss the linear stability of the proposed IMEX schemes for the scalar equation

 (CDβ0u)(t)=λu(t)+ρu(t),t≥0,λ,ρ∈C. (3.1)

We recall the definition of stability for the linear equation (3.1).

###### Theorem 3.1 ([29, 30])

Let . The steady-state solution of Eq. (3.1) is stable if and only if , where .

###### Definition 3.2

A numerical method is said to be -stable if its stability region for (3.1) contains the whole sector .

The following theorem is useful to determine stability regions of the numerical schemes.

###### Theorem 3.3 ([17, 18, 31])

Let . Assume that the sequence is convergent and that the quadrature weights () satisfy

 wn=nβ−1Γ(β+1)+vn,∞∑n=1|vn|<∞, (3.2)

then the stability region of the convolution quadrature is

 ΣNumβ={ξ∈C∣∣1−ξwβ(z)≠0:|z|≤1},wβ(z)=∞∑n=0wnzn,

where or is some function of .

We first consider the linear stability of the IMEX-E scheme for the test equation (3.1). From the IMEX-E scheme (2.26), we get

 Un = U0+(λ+ρ)hβ[n∑k=0ω(β)n−kUk+mu∑k=1W(β,σ)n,kUk+BσnU0] +ρhβω(β)0⎡⎣−Un+2Un−1−Un−2+˜mf∑k=1ˆW(f)n,k(Uk−U0)⎤⎦ = U0+(λ+ρ)hβn∑k=0ω(β)n−kUk−ρhβω(β)0(Un−2Un−1+Un−2)+hβm∑k=0ωn,kUk,

where . By comparing coefficients on both sides of the above identity, we can get . Here we do not give the exact expression of , since it does not affect the stability analysis.

Denote , . Then we have from (3) that

 ∞∑n=2Unzn= U0∞∑n=2zn+(λ+ρ)hβ∞∑n=2(n∑k=0ω(β)n−kUk)zn −ρhβω(β)0∞∑n=2(Un−2U