 # A Fourth-Order Compact ADI Scheme for Two-Dimensional Riesz Space Fractional Nonlinear Reaction-Diffusion Equation

In this paper, a second-order backward difference formula (abbr. BDF2) is used to approximate first-order time partial derivative, the Riesz fractional derivatives are approximated by fourth-order compact operators, a class of new alternating-direction implicit difference scheme (abbr. ADI) is constructed for two-dimensional Riesz space fractional nonlinear reaction-diffusion equation. Stability and convergence of the numerical method are analyzed. Numerical experiments demonstrate that the proposed method is effective.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In this paper, we consider two-dimensional Riesz space fractional nonlinear reaction-diffusion equation (19; 21; 20; 22; 23; 33; 34)

 ∂u(x,y,t)∂t=κ1∂αu(x,y,t)∂|x|α+κ2∂βu(x,y,t)∂|y|β+g(x,y,t,u(x,y,t)),(x,y,t)∈Ω×(0,T], (1.1)

with the boundary and initial conditions

 u(x,y,t)=0,(x,y,t)∈∂Ω×(0,T], (1.2)
 u(x,y,0)=φ(x,y),(x,y)∈¯¯¯¯Ω=∂Ω∪Ω, (1.3)

where and , the diffusion coefficients , are positive constants, is a known sufficiently smooth function, satisfies the Lipschitz condition

 |g(x,y,t,u)−g(x,y,t,υ)|≤L|u−υ|,∀u,υ∈R, (1.4)

here is Lipschitz constant, and Riesz fractional derivatives and are defined by

where , symbols , , and denote left and right Riemann-Liouville fractional derivative operators, which are defined by (10; 11; 12; 13)

 RxDαbu(x,y,t)=1Γ(2−α)∂2∂x2∫bxu(ξ,y,t)(ξ−x)1−αdξ,
 RcDβyu(x,y,t)=1Γ(2−β)∂2∂y2∫ycu(x,ξ,t)(y−ξ)1−βdξ,
 RyDβdu(x,y,t)=1Γ(2−β)∂2∂y2∫dyu(x,ξ,t)(ξ−y)1−βdξ,

where () is Gamma function.

In this paper, we assumed that the problem (1.1)-(1.3) has a unique solution . We also supposed for the fixed and , , for the fixed and , , where and are defined as

 ˜u(x,⋅,⋅)={u(x,⋅,⋅),x∈[a,b],0,R∖[a,b],˜u(⋅,y,⋅)={u(⋅,y,⋅),y∈[c,d],0,R∖[c,d],

here and are of the form

 C4+γ(R)={v|v∈L1(R),∫+∞−∞(1+|ϖ|)4+γ|ˆv(ϖ)|dϖ<∞},γ=α, β,

where

is represented as the Fourier transformation of

and defined by

 ˆv(ϖ)=∫+∞−∞e−iϖxv(x)dx,i2=−1.

In recent years, two-dimensional Riesz space fractional nonlinear reaction-diffusion equation plays an essential role in describing the propagation of the electrical potential in heterogeneous cardiac tissue (19; 20; 33; 21; 22), it attracts many author’s attention in constructing numerical methods for problems of the form (1.1)-(1.3). For approximation of Riesz derivative, Meerschaert and Tadjeran (5) initially proposed the shifted Grünwald-Letnikov approximation with first-order accuracy for Riemann-Liouville fractional derivative. Based on this approximation, Tian et al. (7) estabilished a second-order weighted and shifted Grünwald-Letnikov approximation for Riemann-Liouville fractional derivative, and the approximation was applied in Riesz space fractional advection-dispersion equations (6). Hao et al. (18) constructed a class of new weighted and shifted Grünwald-Letnikov approximation with second-order accuracy, and it was applied in (8) for fractional Ginzburg-Landau equation. Ortigueira (14)

initially proposed the fractional centered difference method with second-order accuracy for Riesz fractional derivative, and this method was applied in Riesz space fractional partial differential equation

(15; 3; 21; 32; 24; 25; 16). Ding and Li (9) proposed a novel second-order approximation for Riesz derivative via constructing a new generating function, and this second-order approximation was adopted in (2) for two-dimension Riesz space-fractional diffusion equation. Recently, compact difference operator has been focused on the fractional differential equations for increasing the spatial accuracy. Zhou et al. (17) constructed a third-order quasi-compact difference scheme for Riemann-Liouville fractional derivative. Hao et al. (18) and Zhao et al. (1) proposed fourth-order compact difference operators to approximate Riemann-Liouville and Riesz derivatives, respectively, these compact difference operators have a great contribution on promoting algorithm accuracy. During these years, there also has developed some approximations by finite element method (22; 28; 29), spectral method (20; 33; 34) et al.. As we noticed, for the approximation of first-order time derivative, implicit Euler method (31; 30; 19; 28), Crank-Nicolson method (32; 22; 20; 18; 15; 6; 2), implicit midpoint method (8; 26) and BDF2 method (27; 3; 4)

are usually used, Padé approximations which are based on Runge-Kutta method are also used in recent researches

(25; 24). And these methods have their own advantages for time-dependent problems.

There are some researches (19; 21; 20; 22; 33; 23; 34) on problem (1.1)-(1.3). Liu et al. (19; 21) constructed two ADI finite difference schemes, where Riesz space derivatives were discretized by shifted Grünwald-Letnikov formulae and fractional centered difference operators, respectively, implicit Euler method was applied to discretize time partial derivative, two proposed methods were proven to be stable and convergent. Bueno-Orovio et al. (33) used Fourier spectral method to approximate Riesz space fractional derivative, implicit Euler method was adopted to discretize first-order time partial derivative, a semi-implicit Fourier spectral method was developed. Zeng et al. (20) and Bu et al. (22) applied Galerkin-Legendre spectral method and Galerkin finite element method to approximate Riesz fractional derivative, respectively, two Crank-Nicolson ADI methods were established. Lin et al. (34) used a bivariate polynomial based on shifted Gegenbauer polynomials method to approximate Riesz space fractional derivative, Runge-Kutta method of order 3 was applied to discretize the first-order time partial derivative, a Runge-Kutta Gegenbauer spectral method was constructed. Iyiola et al. also discussed several implicit-explicit schemes in (23). Because the computing scale of two dimensional diffusion equation problem is very big, so a more efficient algorithm is needed. So far, there are many high-order algorithms for Riesz fractional derivative, we noticed that the fourth-order fractional compact difference operator in (1) is symmetric positive definite under certain circumstance, it’s helpful for us to analyze the stability and convergence. As we know, before ADI method is applied for solving two-dimensional nonlinear reaction-diffusion equation problem. the nonlinear source term needs to have linearized approximation. Therefore, how to deal with the nonlinear source term via linearized approximations (19; 21; 35; 27) plays an important role in constructing ADI scheme for two-dimensional Riesz space fractional nonlinear reaction-diffusion equation. The objective of this paper is to try to use BDF method and the fourth-order fractional compact difference operator to construct a class of new high accuracy ADI scheme based on the first-order (19) and second-order (35) linearized approximations for nonlinear source term. Stability and convergence analysis are given by energy method.

The outline of this paper is organized as follows. In 2, the numerical method is constructed for problem (1.1)-(1.3). Then in 3, stability and convergence are discussed, respectively. In 4, we use the proposed method and these methods in literatures (19; 21) to solve the test problems. Numerical results show that the proposed method has high accuracy and efficiency.

## 2 Numerical method

Let , , , where and are spatial step sizes, denotes time step size. and are exact solution and numerical solution of the problem (1.1)-(1.3) at , respectively. We also denote ,   and the boundary grid mesh is .

To discretize the Riesz space fractional derivative, we would introduce the centred difference operators which are defined by (14)

 Δαxu(xi,yj,tn)=−1hαxi∑k=i−M1g(α)ku(xi−khx,yj,tn), (2.1)

and

 Δβyu(xi,yj,tn)=−1hβyj∑k=j−M2g(β)ku(xi,yj−khy,tn), (2.2)

where the coefficients are determined by

 g(γ)0=Γ(α+1)Γ2(α/2+1),g(γ)k=(1−α+1α/2+k)g(γ)k−1,g(γ)−k=g(γ)k ,γ=α, β,k=1,2,⋯, (2.3)

then we have following lemma.

###### Lemma 2.0.

(see (1).) If ,  , for the fixied step-sizes and , it holds that

 Bαx∂αu(xi,yj,tn)∂|x|α=Δαxu(xi,yj,tn)+O(h4x),1≤i≤M1−1,1≤j≤M2−1,0≤n≤N, (2.4)
 Bβy∂βu(xi,yj,tn)∂|y|β=Δβyu(xi,yj,tn)+O(h4y),1≤i≤M1−1,1≤j≤M2−1,0≤n≤N, (2.5)

where the Fourth-order compact operators and are defined as follows

 Bαxu(xi,yj,tn)={cα2u(xi−1,yj,tn)+(1−2cα2)u(xi,yj,tn)+cα2u(xi+1,yj,tn),1≤i≤M1−1,0≤j≤M2,u(xi,yj,tn),i={0,M1},0≤j≤M2.

and

 Bβyu(xi,yj,tn)={cβ2u(xi,yj−1,tn)+(1−2cβ2)u(xi,yj,tn)+cβ2u(xi,yj+1,tn),1≤j≤M2−1,0≤i≤M1,u(xi,yj,tn),j={0,M2},0≤i≤M1,

where

Before approximating the first-order partial derivative, we would introduce the properties of BDF operator.

###### Lemma 2.0.

(see (3).) For any positive integer , if , then

 ∂u(xi,yj,tn)∂t=D(2)tu(xi,yj,tn)+rni,j, (2.6)

where

 D(2)tu(xi,yj,tn)=⎧⎪⎨⎪⎩δtu(xi,yj,t12),n=1,32δtu(xi,yj,tn−12)−12δtu(xi,yj,tn−32),n≥2,
 δtu(xi,yj,tn−12)=1τ(u(xi,yj,tn)−u(xi,yj,tn−1)),

and satisfies

 ∣∣rni,j∣∣={O(τ),n=1,O(τ2),n≥2. (2.7)

The nonlinear source term can be treated by following process

Let

At the point () (1.1) becomes

 ∂u(xi,yj,tn)∂t=κ1∂αu(xi,yj,tn)∂|x|α+κ2∂βu(xi,yj,tn)∂|y|β+g(xi,yj,tn,u(xi,yj,tn)),(xi,yj)∈Ωh,1≤n≤N. (2.8)

Multiplying by in (2.8), we obtain from 2.1 that

 BαxBβy∂u(xi,yj,tn)∂t=Bβyδαxu(xi,yj,tn)+Bαxδβyu(xi,yj,tn)+BαxBβyg(xi,yj,tn,u(xi,yj,tn))+O(h4x+h4y), (2.9)

where and .

Substituting (2.6) into (2.9), we obtain

 BαxBβyD(2)tu(xi,yj,tn)=Bβyδαxu(xi,yj,tn)+Bαxδβyu(xi,yj,tn)+BαxBβyg(xi,yj,tn,u(xi,yj,tn))−BαxBβyrni,j+O(h4x+h4y). (2.10)

Adding a small error term on both side of (2.10), we have

 BαxBβyD(2)tu(xi,yj,tn)+τ2σ2nδαxδβyD(2)tu(xi,yj,tn)=Bβyδαxu(xi,yj,tn)+Bαxδβyu(xi,yj,tn)+BαxBβyˆg(xi,yj,tn,u(xi,yj,tn))+Rni,j, (2.11)

where and . And there exists the positive constants and such that

 ∣∣Rni,j∣∣≤{c1(τ+h4x+h4y),n=1,c2(τ2+h4x+h4y),n≥2. (2.12)

Omitting the truncation error , we can obtain the numerical scheme for solving the problem (1.1)-(1.3) as follows

 (2.13)

where

the boundary and initial conditions are

 uni,j=0,(xi,yj)∈∂Ωh,1≤n≤N, (2.14)
 u0i,j=φ(xi,yj),(xi,yj)∈¯¯¯¯Ωh. (2.15)

Multipling by in (2.13), and factorizing it, we have

 (Bαx−τσnδαx)(Bβy−τσnδβy)uni,j=Hun−1i,j+τσnBαxBβygni,j, (2.16)

where , and we can also rewrite it as follows

Introducing an intermediate variable , let , therefore, we constructed a class of D’Yakonov ADI finite difference scheme for solving the problem (1.1)-(1.3) as follows

1. for the fixed , can be calculated by

 (Bαx−τσnδαx)u∗i,j=Hun−1i,j+τσnBαxBβygni,j,1≤n≤N, (2.17)

with the boundary conditions

 u∗0,j=(Bβy−τσnδβy)un0,j,u∗M1,j=(Bβy−τσnδβy)unM1,j,1≤j≤M2−1,1≤n≤N. (2.18)
2. for the fixed , can be obtained by

 (Bβy−τσnδβy)uni,j=u∗i,j,1≤n≤N, (2.19)

the boundary and initial conditions are (2.14)-(2.15).

## 3 Stability and convergence analysis

In order to analyze the stability and convergence of the method, we introduce some notations and lemmas.

Let

 ˆγh={ζn∣∣ζn=(ζn0,0,⋯,ζnM1,0,⋯,ζn0,M2,⋯,ζnM1,M2),  ζni,j=0 if (xi,yj)∈∂Ωh,0≤n≤N}.

For any , we define the following discrete inner product and corresponding norm

 (un,vn)=hxhyM1−1∑i=1M2−1∑j=1uni,jvni,j,∥un∥=√(un,un).

In addition, for any , we denote , it implies , and denote and .

###### Lemma 3.0.

(see (3)

.) For any positive integer n and real vector

, we have

 4τ3n∑k=2vk(D(2)tvk)≥(vn)2−13(vn−1)2−(v1)2+13(v0)2−23(v1−v0)2,n≥2,4τ3n∑k=1vk(D(2)tvk)≥(vn)2−13(vn−1)2−13(v1)2−13(v0)2,n≥1,4τ3v1(D(2)tv1)=23(v1)2−23(v0)2+23(v1−v0)2,n=1.

It is easy to check that and are symmetric positive definite and self-adjoint (1), Following from Lemma 3.11 in (1), one can prove that there exists the fractional symmetric positive definite difference operators and such that and , here, and are also commutable.

###### Lemma 3.0.

(see (1).) For any it holds that

 13∥un∥2≤∥un∥2B≤∥un∥2,

where .

###### Lemma 3.0.

(see (3; 1; 15).) For any it holds that

 (δαxun∗,j,un∗,j):=−(Λxun∗,j,Λxun∗,j)≤0,1≤j≤M2−1,(δβyuni,∗,uni,∗):=−(Λyuni,∗,Λyuni,∗)≤0,1≤i≤M1−1,

where , are represented as fractional symmetric positive definite difference operators such that , .

It is easy to verify that

 (δαxδβyun,un)≥0.

Theorem 4.1 in (15) means that and are symmetric positive definite operators. Therefore, with the help of commutativity of and , we could introduce a semi-norm be similar with (1; 3).

###### Lemma 3.0.

For any it holds that

 (Bβyδαxun,un)≤0,(Bαxδβyun,un)≤0.
###### Proof.

From 3.5, we have

 (Bβyδαxun,un)=(δαxQyun,Qyun)=hxhyM1−1∑i=1M2−1∑j=1(δαxQyuni,j)(Qyuni,j) =hyM2−1∑j=1(δαxQyun∗,j,Qyun∗,j) ≤0.

Similarly, we can obtain

 (Bαxδβyun,un)≤0.

The proof is completed. ∎

###### Lemma 3.0.

(Discrete Bellman Inequality)   Let  are a series of nonnegative real numbers, satisfying

 ϵn≤ρ2+ρ1τn−1∑k=0ϵk,n=1,⋯,ˆN,

then it holds that

 ϵn≤ρ2eρ1nτ,n=1,⋯,ˆN.

Assuming that is the numerical solution for the numerical method (2.13)-(2.15) starting from another initial value. Denote , where , then we have following consequence.

###### Theorem 3.8.

For any positive real number , if , then the numerical scheme (2.13)-(2.15) is stable, i.e.

 ∥En∥≤3νe18νLT∥(I+2318τ2δαxδβy)E0∥,n≥1.
###### Proof.

According to (2.13)-(2.15), we have the following equations

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩BαxBβyD(2)tEni,j+τ2σ2nδαxδβyD(2)tEni,j=BβyδαxEni,j+B