Numerical study of the nonlinear anomalous reaction–subdiffusion process arising in the electroanalytical chemistry

Highlights

• The paper presents a numerical technique for a nonlinear anomalous reaction–subdiffusion process.

• An implicit time stepping scheme with second-order accuracy is used to discretize the temporal derivatives.

• The RBF-FD method is adopted for approximating the spatial derivatives over local collection nodes.

• The time-discretized scheme is proved to be stable and convergent using the energy method.

• Numerical results illustrate the efficiency and the outperformance of the proposed approach.

Abstract

This paper presents a meshless method based on the finite difference scheme derived from the local radial basis function (RBF-FD). The algorithm is used for finding the approximate solution of nonlinear anomalous reaction–diffusion models. The time discretization procedure is carried out by means of a weighted discrete scheme covering second-order approximation, while the spatial discretization is accomplished using the RBF-FD. The theoretical discussion validates the stability and convergence of the time-discretized formulation which are analyzed in the perspective to the $H^1$-norm. This approach benefits from a local collocation technique to estimate the differential operators using the weighted differences over local collection nodes through the RBF expansion. Two test problems illustrate the computational efficiency of the approach. Numerical simulations highlight the performance of the method that provides accurate solutions on complex domains with any distribution node type.

Introduction

In recent years, research in the field of fractional calculus (FC) attracted an increasing attention due to its wide spectrum of applications [1], [2]. This interest stems from the modeling of physical systems via FC tools, such as in diffusion systems, control theory, heat conduction, and electromagnetics [3], [4], [5], [6]. Several real-world phenomena can be better described and interpreted with fractional mathematical models due to their ability to capture long-term memory and non-local effects. The fractional equations have been studied using different methods. Liu et al. [7] developed a Galerkin mixed finite element (FE) technique combined with a time second-order discrete approach to approximate the nonlinear time fractional diffusion model. Chen et al. [8] proposed a two-grid modified method with characteristics based on the FE for the nonlinear variable-order time-fractional mobile/immobile advection–diffusion model. Wong et al. [9] adopted a Galerkin FE algorithm combined with a second-order finite difference (FD) approach for solving the nonlinear time fractional cable model. Manimaran et al. [10] solved the time-fractional non-local diffusion model by using the Galerkin FE for the space and a FD for the time-fractional derivative. Thakoor and Bhuruth [11] formulated a local radial point interpolation scheme for the time-fractional diffusion and diffusion damped wave problems. Garshasbi et al. [12] presented an iterative procedure based on an implicit FD to approximate a space–time fractional moving boundary problem.

Advanced numerical techniques are of key relevance in the field of fractional diffusion models. Such models arise in the unification of wave propagation and diffusion phenomena, the modeling of anomalous subdiffusive and diffusive systems, and the continuous-time random walk (CTRW) problem [13]. Diffusion processes, such as the Gaussian statistics using Fick’s second law, or the Brownian motion, described through the mean squared displacement (MSD) i.e., $〈x^2\left(t\right)〉\sim {\kappa }_1{t}^{},t\to \infty$, where $k_1$ stands for the diffusion coefficient and the symbol $〈\cdot 〉$ represents an average either over time or over an ensemble of similar particles. These models are not the most adequate for the time duration in different complex systems because of their deviation from the linear time dependence. Also, stochastic processes work poorly because of the central limitation of the Markov theorem. This type of phenomenon is named anomalous diffusion. Anomalous diffusion has been widely adopted since it exhibits a non-linear growth in what regards the MSD during the time period. The diffusion is characterized by a mean-square displacement in a nonlinear growth in time $t$, i.e., $〈x^2\left(t\right)〉\sim {\kappa }_{\alpha }{t}^{\alpha },t\to \infty$, where ${\kappa }_{\alpha }$ denotes the generalized diffusion coefficient and $\alpha$ can have the values $0<\alpha <1$, $\alpha =1$, $1<\alpha <2$ and $\alpha =2$, for the anomalous subdiffusion, normal diffusion, super-diffusion, and ballistic diffusion, respectively. It is worth mentioning that the generalized diffusion coefficient ${\kappa }_{\alpha }$ has the dimension $\left[{\kappa }_{\alpha }\right]={\mathrm{cm}}^2{\mathrm{s}}^{-\alpha }$ [13]. The CTRW models with temporal and spatial memory unify these types of anomalous diffusion [13].

Oldham and co-workers [14], [15], [16], [17] played a pioneer role in electrochemistry where the adoption of fractional-order derivatives and integrals lead to relevant results. The adoption of a semi-integral of the current $i\left(t\right)$, ${}_0{\mathrm{D}}_t^{-\frac{1}{2}}i\left(t\right)$, in the work [15] on semi-integral electroanalysis paved the way to new ideas and the research was later continued with the semi-differential electroanalysis method [18]. Determining the concentration of the electro-active species close to the surface of the electrode is an important topic. The technique proposed by Oldham and Spanier [17] includes a boundary (electrode surface) relationship to replace a diffusion model problem under specific conditions. Based on this concept, Oldham [15] considered the relation

$$m\left(t\right)={}_0{\mathrm{D}}_t^{-\frac{1}{2}}i\left(t\right),$$

that models the electrode surface as the semi-integral of the current $i\left(t\right)$, where $m\left(t\right)$ represents the fractional integral of the faradaic current $i\left(t\right)$ that flows when a complete diffusion control follows a pre-existing equilibrium [19]. Then, the surface concentration $C_s\left(t\right)$ of the electroactive species is given by

$$C_s\left(t\right)=C_0-k{}_0{\mathrm{D}}_t^{-\frac{1}{2}}i\left(t\right),$$

where $k$ is a parameter to be defined below [15] and the constant $C_0$ represents the uniform concentration of the electroactive species. The expression (2) was determined following the normal diffusion model [19]

$$\frac{\partial C(x,t)}{\partial t}={\mathrm{D}}_{\ast }\left[\frac{{\partial }^{2}C(x,t)}{\partial x^2}\right],t>0,0<x<\infty ,$$$$C(x,0)=C_0,C(\infty ,0)=C_0,$$$${\mathrm{D}}_{\ast }{\left[\frac{\partial C(x,t)}{\partial t}\right]}_{x=0}=\frac{i\left(t\right)}{n\mathcal{A}\mathcal{F}},$$

where ${\mathrm{D}}_{\ast }$ denotes the diffusion coefficient, $\mathcal{F}$ is the Faraday’s constant, $n$ represents the number of electrons associated to the reduction of the species, $i\left(t\right)$ stands for the cathodic faradaic current and $\mathcal{A}$ is the electrode area. Therefore, the coefficient $k$ in (2) can be represented as

$$k=\frac{1}{n\mathcal{A}\mathcal{F}\sqrt{{\mathrm{D}}_{\ast }}}\cdot$$

The fractional order diffusion model [2] can be considered as an alternative of the normal diffusion model (3), and can be written as

$$\frac{\partial C(x,t)}{\partial t}={}_0{\mathrm{D}}_t^{1-\alpha }\left[\frac{{\partial }^{2}C(x,t)}{\partial x^2}\right],$$

where ${}_0{\mathrm{D}}_t^{1-\alpha }$ is the Riemann–Liouville (R–L) fractional partial derivative of order $1-\alpha$

$${}_0{\mathrm{D}}_t^{1-\alpha }u(x,t)=\frac{1}{\Gamma \left(\alpha \right)}\frac{\partial }{\partial t}{\int }_0^t\frac{u(x,\varsigma )}{{(t-\varsigma )}^{1-\alpha }}\mathrm{d}\varsigma ,$$

and $\Gamma \left(\cdot \right)$ is the gamma function.

The nonlinear anomalous reaction–diffusion process plays a important and fundamental role not only in chemistry and physics, but also in the modeling of pattern formation that arise in biological, sociological and ecologic systems. A fractional reaction–diffusion model was derived by Henry and Wearne [20] from a CTRW with temporal sources and memory. This description presents a comprehensive model for reaction–diffusion processes accompanied by anomalous diffusion like those exhibited by spatially inhomogeneous environments. A set of continuum fractional diffusion models was proposed by Yuste et al. [21] for describing the reaction front for the $A+B⟶C$ reaction–subdiffusion process. It was shown by Sokolov et al. [22] that we cannot have some close resemblance between the reaction–subdiffusion and the reaction–diffusion models. Moreover, it was considered that the former cannot be derived via a slight change in the diffusion operator of a subdiffusion model. In this article, our target is to study the approximated solution of the nonlinear anomalous reaction–subdiffusion model (ARSDM) in the form:

$$\frac{\partial u(\mathbf{x},t)}{\partial t}={}_0{\mathrm{D}}_t^{1-\alpha }[{\kappa }_1\Delta u(\mathbf{x},t)-{\kappa }_{2}u(\mathbf{x},t)]+\mathcal{G}\left(u(\mathbf{x},t)\right)+f(\mathbf{x},t),(\mathbf{x},t)\in \Omega \times (0,T].$$

The initial and Dirichlet boundary conditions can be considered as

$$u(\mathbf{x},0)=g\left(\mathbf{x}\right),\mathbf{x}\in \Omega \cup \partial \Omega ,$$$$u(\mathbf{x},t)=\psi (\mathbf{x},t),\mathbf{x}\in \partial \Omega ,t>0,$$

in which $\mathbf{x}=(x,y)$, $\alpha \in (0,1)$, ${\kappa }_1$ (the coefficient of anomalous diffusion) and ${\kappa }_2$ are positive constants, $\Omega$ represents a bounded convex polygonal domain with boundary $\partial \Omega$, $\Delta$ denotes the Laplacian operator $\Delta u=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}$, $T$ is final time, and $g\left(\mathbf{x}\right)$ and $\psi (\mathbf{x},t)$ are given functions with sufficient smoothness. The non-linear source term $\mathcal{G}\left(u\right)$ has the second-order continuous partial derivative, satisfies the assumption $\left|\mathcal{G}\left(u\right)\right|\le L\left|u\right|$, and the first order derivative of $\mathcal{G}\left(u\right)$ with respect to $u$ is bounded.

In recent years, several numerical techniques for solving the ARSDM were formulated. Chen et al. [23] proposed a difference approximation scheme to solve the fractional diffusion model describing sub-diffusion and investigated its stability through the Fourier method. Baeumer et al. [24] constructed a technique based on operator splitting to solve the ARSDM. Zhuang et al. [25] derived an implicit FD method for the ARSDM. Yu et al. [26] designed an implicit compact FD scheme, while Li and Ding [27] adopted a high order FD for the ARSDM. Dehghan et al. [28] employed a meshless Galerkin method and Zhu and Xie [29] applied a Galerkin finite element scheme coupled with a first-order implicit–explicit algorithm. Abbaszadeh and Dehghan [30] combined the alternating direction implicit (ADI) technique with the interpolating element free Galerkin (IEFG) method, whereas Jafarabadi and Shivanian [31] implemented the radial point interpolation method (RPIM). Ghehsareh et al. [32] formulated a meshless algorithm based on the RPIM to approximate the ARSDM.

Researchers have been exploring meshless techniques as an alternative for mesh generation, where a group of scattered nodes is employed. The radial basis function (RBF) is a powerful method for interpolating scattered data [33], [34], [35]. The RBF involves a function that includes an independent variable for the distance between the center point and the calculation point. This function is not related with the dimension of the space. In addition, the RBF has a simple symmetrical form and the multivariate problem can be easily transformed into a one-variable problem. Among the advantages of the RBF we can mention the spectral convergence, dimension insensitivity, needless for node connectivity, and easy implementation. Some progress has taken place in the last decade in the application of the RBF for the solution of different types of partial differential equations (PDE). In general, the RBF is superior to the FE method [36], because: (i) it is truly meshless since one can select the collocation points freely (unlike in the FE, the points are not required to be connected) and, therefore, one can avoid complex meshing exercises; (ii) it is an independent spatial coordinate that can be readily developed in order to solve multi-dimensional problems; (iii) does not requires grid generation saving, therefore, pre-processing time.

During the recent years, meshless RBFs have been widely used as a potential alternative for solving PDEs problems in various applications [37], [38], [39]. The meshless RBF can be classified into globalized and localized approaches. The interpolation matrix of the global versions of the RBF (GRBF) is always dense and leads to high condition numbers for an increasing number of nodes (also known as supporting centers). Therefore, we obtain an ill-conditioned and large linear system that can lead to uncertain results. Conversely, these methods have sparse matrices with better condition number. In many cases, the situation can be improved by introducing localization approaches with better conditioned matrices, see e.g., [40], [41], [42], [43], [44]. The first strategy, which is the one adopted in this paper, is the natural generalization of the classical grid-based FD, consisting of a combination the FD and RBF that is called RBF-FD [43], [44]. The RBF-FD is drawing attention due to its limited computational complexity in comparison with other RBF methods, besides its high-order accuracy and parallel nature. The second strategy is to adopt the RBF based on a partition of unity method, which was proposed by Babuška and Melenk [45]. A partition based formulation is also well suited for parallel implementation [42], [46]. Interested readers can see further works for the parallelization of localized RBF techniques in [40], [46], [47], [48] and tackling the ill-conditioning of the matrices by the RBF-QR technique in [49], [50].

The main objective of this paper is to develop a hybrid algorithm using the local RBF-FD collocation technique for obtaining the approximate solution of the nonlinear ARSDM. One interesting issue is put forward in this paper, namely about the accuracy order of the time discretization, due to the importance of the numerical approximation of the time derivative in the fractional diffusion equations. In this line of thought, a second-order time discretization approach is implemented. It is worth to mention that the order of convergence in time of the proposed method is superior to the one of previous schemes found in the literature. In fact, we improve the time accuracy of [25], [28], [29], [30], [31] to second order. The RBF-FD produces a specific accuracy level while requiring a considerably low computational time when compared to global collocation techniques, besides making use of much better conditioned matrices.

This paper is organized into five sections. Section 2 proposes the approximation in time by the FD scheme. In addition, the time-discrete algorithm in terms of unconditional stability and convergence is examined. Section 3 implements the RBF-FD method to discretize the spatial variable. Section 4 dwells with the numerical verification on the applicability and effectiveness of the proposed techniques in the irregular domains. Section 5 summarizes the main aspects and presents the concluding remarks.