High-order time-stepping methods for two-dimensional Riesz fractional nonlinear reaction–diffusion equations

Highlights

• Two 4th order ETD methods are presented: L-stable method ETDRK04 and A-stable method ETDRK22.

• Allen Cahn equation with cubic nonlinearity and smooth as well as nonsmooth initial data.

• Fisher’s equation with quadratic nonlinearity, Enzyme Kinetics equation with rational nonlinearity.

• Diffusion profiles confirm the efficiency of the L-stable method for non-smooth initial data.

• Computational efficiency of these methods is also demonstrated.

Abstract

Anomalous diffusion and long-range spatial interactions in anisotropic media could be captured by considering the Riesz fractional derivatives rather than the classical Laplacian. While analytical solutions of the resulting fractional reaction–diffusion models may not be available, the numerical methods for approximating them are challenging. In this paper, fourth-order $L$-stable (ETDRK04) and $A$-stable (ETDRK22) linearly implicit predictor-corrector type time-stepping methods are presented. These methods are implemented on two-dimensional Riesz fractional nonlinear reaction–diffusion equations with smooth and non-smooth initial data. Three types of nonlinear reaction–diffusion models are considered: Allen–Cahn equation with cubic nonlinearity, Fisher’s equation with quadratic nonlinearity, and Enzyme Kinetics equation with rational nonlinearity. Fourth-order temporal convergence rate of the methods is proved analytically and computed numerically. Profiles of the numerical solutions corresponding to different orders and rates of diffusion are included. Computational efficiency of the ETDRK04 and ETDRK22 methods over the well known Cox–Matthews ETDRK4 is presented. The superiority of the provided methods, in terms of computational accuracy, efficiency, and reliability, is demonstrated through the numerical experiments.

Introduction

In recent years, considerable interest is shown to model the phenomena of anomalous diffusion due to long-range spatial interaction. Phenomena of anomalous diffusion are observed in many fields of science and engineering. For example, modeling control theory [1], mechanical properties of materials [2], groundwater flow [3], [4], anomalous diffusion in porous media [5], [6]

Although analytic solutions of most of such models are not available or too complicated, yet many analytical techniques for solving linear fractional differential equations are available in the literature, for example, Green’s function and transform methods. However, the analytical solutions are presented in terms of some special functions which are difficult to compute directly.

Several numerical methods have been developed for solving space-fractional partial differential equations. Aceto and Novati [7] developed a numerical method by using rational approximation of the fractional Laplacian operator in reaction–diffusion problems. Khaliq et al. [8] have developed second-order linearly implicit predictor–corrector methods for space-fractional reaction–diffusion equations with nonsmooth initial data. It is shown that the methods based on Matrix Transfer Technique for spatial discretization are unconditionally stable. Furati et al. [9] have developed fourth-order methods for one-dimensional space-fractional reaction–diffusion equations with nonsmooth initial data.

In this paper we consider the following initial boundary value problem

$$\frac{\partial u}{\partial t}=K_x\frac{{\partial }^{\alpha }u}{\partial {\left|x\right|}^{\alpha }}+K_y\frac{{\partial }^{\beta }u}{\partial {\left|y\right|}^{\beta }}+F\left(u\right)+f(x,y,t),(x,y)\in \Omega ,t>0,$$$$u(x,y,0)=\varphi (x,y),(x,y)\in \Omega ,$$$$u(x,y,t)=0,(x,y)\in \partial \Omega ,t>0,$$

where $\Omega =(a,b)\times (c,d)$, $1<\alpha ,\beta \le 2,K_x>0,K_y>0$. The functions $F$ and $f$ are sufficiently smooth functions.

The spatial fractional derivatives in (1.1) are the Riesz derivatives,

$$\frac{{\partial }^{\alpha }}{\partial {\left|x\right|}^{\alpha }}u(x,y,t)=-\frac{1}{2cos\frac{\pi \alpha }{2}}\left[{}_a{D}_x^{\alpha }{+}_x{D}_b^{\alpha }\right]u(x,y,t),$$$$\frac{{\partial }^{\beta }}{\partial {\left|y\right|}^{\beta }}u(x,y,t)=-\frac{1}{2cos\frac{\pi \beta }{2}}\left[{}_c{D}_y^{\beta }{+}_y{D}_d^{\beta }\right]u(x,y,t),$$

where the left Riemann–Liouville fractional derivatives ${}_a{D}_x^{\alpha }$, ${}_c{D}_y^{\beta }$ and the right Riemann–Liouville fractional derivatives ${}_x{D}_b^{\alpha }$, ${}_y{D}_d^{\beta }$ are defined by,

$${}_a{D}_x^{\alpha }u(x,y,t)=\frac{1}{\Gamma (2-\alpha )}\frac{{\partial }^2}{\partial x^2}{\int }_a^x{(x-\xi )}^{1-\alpha }u(\xi ,y,t)d\xi ,$$$${}_c{D}_y^{\beta }u(x,y,t)=\frac{1}{\Gamma (2-\beta )}\frac{{\partial }^2}{\partial y^2}{\int }_c^y{(y-\eta )}^{1-\beta }u(x,\eta ,t)d\eta ,$$$${}_x{D}_b^{\alpha }u(x,y,t)=\frac{1}{\Gamma (2-\alpha )}\frac{{\partial }^2}{\partial x^2}{\int }_x^b{(\xi -x)}^{1-\alpha }u(\xi ,y,t)d\xi ,$$

and

$${}_y{D}_d^{\beta }u(x,y,t)=\frac{1}{\Gamma (2-\beta )}\frac{{\partial }^2}{\partial y^2}{\int }_y^d{(\eta -y)}^{1-\beta }u(x,\eta ,t)d\xi .$$

The advantage of model (1.1) is that it allows for different diffusion rates and orders with respect to each space variable. The Riesz fractional derivative is introduced to describe the generalized diffusive fluxes that do not satisfy Fick’s law, Fourier law, or Newtonian constitutive equation. These nonlocal derivatives can enhance the predictivity of some models since they couple the local state of the field with the state of the entire field. More details on Riesz fractional derivative and fractional Laplacian can be found in [10].

Several numerical methods for two-dimensional reaction–diffusion problems with Riesz space-fractional derivative have been developed. Fawang [11] proposed an alternating direction implicit (ADI) method using first-order backward difference method in time. Zeng [12] developed an ADI Galerkin–Legendre spectral method and implemented it to solve two-dimensional Riesz space-fractional nonlinear reaction–diffusion equation using Crank–Nicolson method. Yang [13] presented a fully discrete scheme by applying Galerkin finite element method in space and first order backward difference method in time. Zhu [14] used second-order finite element method in space with Crank–Nicolson method in time. The method was implemented to solve two-dimensional nonlinear Fisher’s equation with Riesz fractional derivatives in space. Cheng [15] developed a compact ADI scheme for two-dimensional Riesz space-fractional nonlinear reaction–diffusion equations. They showed fourth-order convergence in space and second-order convergence in time. The method was implemented to solve two-dimensional Allen–Cahn and Fisher’s equation. Xing [16] developed a fourth-order difference scheme for two-dimensional Riesz space-fractional diffusion equations. Fourth-order in space and second-order convergence in time was shown. All the methods mentioned above are either first-order or second-order accurate in time.

Cox and Matthews [17] developed a class of Exponential Time Differencing Runge–Kutta (ETDRK) methods for stiff systems of ordinary differential equations. Kassam and Trefethen [18] showed that severe cancellation of errors can occur which affect the accuracy and stability of these methods. Cox–Matthews methods require computing matrix exponential functions and inverting matrix polynomials. Computing matrix exponential functions is computationally expensive because it will not be a sparse matrix even if the original matrix is sparse. Moreover, these methods suffer from numerical instability when inverting the polynomials of matrices. Yousuf et al. [19] resolved these issues by using Padé approximations of the matrix exponential functions. They developed a fourth order $L$-stable method (ETDRK04) by modifying the Cox–Matthews fourth order (ETDRK4) method using the fourth order (0,4)-Padé approximation to the matrix exponential functions. Partial fractions decomposition technique is applied to develop an efficient version of the method which reduces the computational cost by almost 50%.

Recently some methods have been developed based on ETD approach for solving space-fractional reaction–diffusion equations. Khaliq et al. [20] developed a fourth-order method for systems of space-fractional nonlinear Shrödinger equations. The method combines a fourth-order discretization for the Riesz derivative developed by Ding et al. [21] with the fourth-order ETDRK22 method developed by modifying ETDRK method using (2,2)-Padé approximation. Ding and Zhang [22] also used (2,2)-Padé approximation to compute the matrix exponential function and developed a fourth-order time-stepping method for Riesz fractional partial differential equations in one space dimension only. For one-dimensional space-fractional reaction-diffusion equations, Furati et al. [9] developed two fourth-order methods: an $L$-stable ETDRK04 method using (0,4)-Padé and an $A$-stable ETDRK22 method using (2,2)-Padé approximations.

In this paper, we present fourth-order methods in time to solve the two-dimensional Riesz fractional nonlinear reaction–diffusion equations. Classical convergence results typically rely on smoothness assumptions on the underlying data. The nonsmooth data leads to serious degradation in the convergence of the numerical methods unless certain criteria are used. The presented ETDRK04 method works well and no unwanted oscillation are produced even if the initial data is nonsmooth. The second method we present is ETDRK22 which works well when initial data is smooth. We show through examples that these methods yield fourth-order accuracy in time.

The paper is organized as follows. The discretization of Riesz derivative is described in Section 2. In Section 3, the time-stepping methods are explained and implementation algorithms are given. Linear Error analysis of the methods is given in Section 4. Error analysis for the nonlinear reaction–diffusion problem is given in Section 5. Numerical experiments with graphs and convergence tables are provided in Section 6. Finally, a conclusion is included in Section 7.