Numerical analysis of time-fractional Sobolev equation for fluid-driven processes in impermeable rocks

This paper proposes a local meshless radial basis function (RBF) method to obtain the solution of the two-dimensional time-fractional Sobolev equation. The model is formulated with the Caputo fractional derivative. The method uses the RBF to approximate the spatial operator, and a finite-difference algorithm as the time-stepping approach for the solution in time. The stability of the technique is examined by using the matrix method. Finally, two numerical examples are given to verify the numerical performance and efficiency of the method.


Introduction
Extending integer derivatives to fractional orders in differential equations led to an area in mathematics called fractional calculus (FC) [1][2][3][4]. This field, dealing with the calculus of derivatives and integrals of arbitrary order, became popular in the last thirty years, since numerous significant phenomena in signal processing and optics [5][6][7], biological systems [8], quantum mechanics [9], electrochemistry [10], fluid mechanics [11], viscoelasticity [12], and electromagnetics [13] can be described via fractional differential equations (FDEs). Moreover, FC is presently a key tool in the modeling of particle transport occurring in porous heterogeneous media in addition to complex phenomena [1,2]. Therefore, researchers need an effective tool to solve FDEs; however, obtaining exact solutions for such equations is difficult. Hence, it is becoming increasingly important to develop efficient numerical techniques to tackle these problems [14]. Nevertheless, numerical schemes for FDEs lead to mathematical difficulties not faced in integer-order model analysis. The main reasons for this are that the fractional difference operators have a nonlocal nature and their adjoints are not the negatives of themselves [15].
From the perspective of the numerical analysis, there are some fundamental difficulties in numerically solving the fractional derivatives, because some of the good properties of classical approximating operators are lost. During the last decades, the difference method has made some developments for approximating FDEs, [16,17]. The Riemann-Liouville fractional derivative can be discretized by the standard Grünwald-Letnikov formula [18] with only first-order accuracy, but the difference scheme based on the Grünwal-Letnikov formula for time-dependent problems is unstable [16]. To tackle this problem, Meerschaert and Tadjeran [16] presented the shifted Grünwald-Letnikov formula to simulate fractional advection-dispersion flow equations. Sousaa and Li [19] adopted a secondorder discretization for the Riemann-Liouville fractional derivative and established an unconditionally stable weighted average difference method. Ortigueira [20] advanced the fractional centred derivative to solve the Riesz fractional derivative with second-order accuracy. It is worth mentioning that the singularity of the exact solution for time FDEs often occurs near the starting time t = 0 [21]. Therefore, one needs to establish some proper regularity assumptions to develop the numerical schemes [21]. As a result, constructing fast and high-order numerical techniques to solve FDEs is a challenging task.
Hereafter, we study a numerical approach to solve the time-fractional Sobolev equation (TFSE) as with the initial condition (IC) and the boundary condition (BC) where γ and σ denote two positive constants, expresses a continuous domain in R 2 with boundary ∂ , ∇ 2 represents a Laplacian operator in the variable space, T stands for the total time, and u(x, t) is the unknown solution to be determined. The time-fractional derivative in the TFSE denotes the Caputo fractional derivative of order α expressed by Here, the process of modeling fluids is developed using the fractional derivative operator. This generalization includes by means of the Caputo fractional derivative instead of the standard temporal derivative. When α = 1, the TFSE yields the classical pseudoparabolic or Sobolev equation that describes various fluid-mechanics and engineering problems, e.g., quasistationary phenomena in semiconductors, thermal conduction with conductive or thermodynamic temperature, and the flow of fluid through fissured rocks [22]. Indeed, the Sobolev equation is the mathematical model of the vertical nonstationary groundwater flow with dynamic capillary-pressure effect in a porous medium (see [23] and references therein). A noteworthy characteristic of these models consists of their capability of expressing the conservation of some quantities, such as mass, heat, and momentum [22]. The Sobolev model is a special class of the Benjamin-Bona-Mahony-Burgers (BBMB) equation, where the coefficients of the nonlinear term and first-order derivatives equal zero [24].
Some numerical techniques have been advanced to approximate the TFSE. Liu et al. [25] applied a modified reduced-order finite-element (FE) technique. Haq and Hussain [26,27] developed a meshless scheme based on radial basis functions (RBFs). Beshtokov [28] proposed a finite-difference algorithm, while Qin et al. [29] employed a Newton linearized scheme based on the Crank-Nicolson technique and Zhao et al. [30] used a finite-volume element (FVE) approach.
We find in the technical literature three standard numerical approaches, namely the finite element, the finite difference, and the finite volume (abbreviated, FE, FD, and FV) schemes. These approaches have been successfully used to approximate a variety of partial differential equations (PDEs) in diverse subject areas, such as electromagnetism, fluid dynamics, material science, financial markets, and astrophysics. In fact, the behavior of some materials in nature may be modeled via PDEs and solved by traditional numerical methods. Examples include climate and weather modeling performed in geophysics, Navier-Stokes equations arising in fluid dynamics, biphasic modeling in engineering, and Maxwell's equations used in electrodynamics. The FD is common in science and engineering for its simplicity. However, the main shortcoming of the FD is its inability to handle higher-dimensional geometries since the PDE discretization depends on a topological line grid. On the other hand, the FE is the most flexible traditional method with respect to geometry. An alternative strategy is the so-called spectral method. This is an accurate method, but it suffers from geometric restrictions and has predetermined periodic boundary conditions in the Fourier case. The computational domain in the FE is partitioned into smaller subdomains and the solution is constructed in each element using the basis functions. Consequently, the number of dimensions or variables may reach hundreds or thousands in a real-world problem. Hence, the question arises as to whether it is possible to generate grid points in problems with irregular domains and complex geometries in higher dimensions. This problem led to the development of meshless (or mesh-free) methods that are meshless in the sense that they do not require point connectivity in the mesh/grid. Determining the nearest neighbors involves a smaller computational cost than mesh generation in traditional methods.
The main idea behind the RBF-PU is to decompose the original domain into several covering subdomains (patches) and to construct a local RBF approximant over each subdomain. The localized method overcomes the large computational cost posed by the illconditioned and dense matrices of the GRBF method, while maintaining the accuracy resulting from the sparsity of such matrices. Moreover, the RBF-PU also achieves good accuracy with significantly less computational burden than the GRBF. This paper is outlined as follows. Section 2 provides a time-discrete scheme using the finite-difference (FD) formula through the θ -weighted rule. Section 3 discusses the spatial discretization using the RBF-PU. In addition, the stability of the RBF-PU collocation through the matrix method is also analyzed. Section 4 presents two numerical problems that exemplify the accuracy and efficiency of the RBF-PU. Finally, Sect. 5 highlights the main conclusions.

Spatial derivative discretization
This section describes the LRBF-PU collocation scheme. The PU technique determines the numerical solution via a weighted sum of the local approximants created on overlapping subdomains. The RBFs are applied as local approximants and Wendland's compactly supported RBFs [46] are adopted as weight functions. This strategy is based on constructing differentiation matrices that convert a given PDE into an algebraic system. The benefit of RBF-PU comes from its low computational cost owing to its relatively sparse matrices. Let = {x 1 , . . . , x N } ⊆ R d , be a set of scattered data sites, the RBF interpolant S(x) of u(x) to the data points u j = u(x j ), j = 1, . . . , N , takes the following form: in which a j represent the coefficients to be determined from the data, φ j (x, ε) = φ( xx j 2 , ε), j = 1, . . . , N , are the RBFs, r = x -x j 2 stands for the Euclidean norm and x j are centers coinciding with the collocation points. The constant ε denotes the shape parameter, which is responsible for the flatness of the functions. To compute the unknown coefficients λ j , we can impose the interpolation conditions S(x i ) = u i , i = 1, . . . , N ) and as a result we obtain a N × N linear system Thus, the alternative expression for the interpolant (9) can be illustrated as: in which (x) = [ψ 1 (x), . . . , ψ N (x)] T . In view of Eqs. (12), (10), and (9), we deduce the relation between the original basis and the Lagrange radial basis as The nonsingularity of matrix A guarantees that the transformation (13) is valid. Let ⊂ R 2 be an open and bounded domain, and = {x 1 , . . . , x N } ⊆ R d a set of collocation nodes. We partition ⊂ R 2 into M patches or subdomains { j } M j=1 of the domain so that ⊂ M j=1 j . Figure 1 represents a schematic diagram of the square domain along with the related circular subdomains. Also, we define in which the constant C is independent of the number of subdomains [46]. The PU weight functions w j can be established based on Shepard's method [47] as follows: in which ϕ j (x) is a compactly supported function on each subdomain j . It follows that w j (x) = 0 ∀j / ∈ I(x). To ensure the nonnegativity and compact support in subdomain j , we consider in (14) ϕ j (x) = ϕ j x -x j ρ j , j = 1, 2, . . . , M.
Here, ρ j denotes the radius of the subdomain j andx j represents its center node. The function ϕ j denotes one of the compact functions with minimal degree [48]. Hereafter, we adopt the Wendland function ϕ(r) = (1 + 4r) + (1r) 4 , where r = x-x j ρ j [46].
The global approximant is thus constructed in the computational domain as follows where {S j } k j=1 defines a local RBF interpolant on each subdomain j . We approximate a spatial differential operator ∂ |σ | ∂x σ at the interior nodes in order to use the LRBF-PU for the spatial discretization of the PDE as where σ and ϑ ∈ N d 0 . The evaluation of the Laplace operator can be represented as From (16), we express the approximation solution V k by the following expression: We divide the collocation points into two sets, J and I, considered to be sets of boundary and internal points, respectively. Moreover, the overall number of points is equal to N with N = N J + N I , so that N = N I and N = N J denote the total numbers of internal and boundary points, respectively. We split matrix A into 2 matrices, A I and A J , as follows in which Inserting Eq. (19) into Eq. (8) gives a system of N linear equations in the matrix form When k ≥ 1, we obtain where Using system (22), we can calculate the numerical solution at any temporal step n.

Stability analysis
Here, we investigate the stability of the RBF collocation approximation (22) based on the matrix method. Let us define the error e n at the nth time level as where u n and U n are the analytic and approximate solutions at the nth time level, respectively. Then, the error equation of (22) can be written as where P = AA -1 BA -1 is called the amplification matrix. Due to the Lax-Richtmyer definition of stability, Eq. (24) is said to be stable if P ≤ 1, which is equivalent to ρ(P) ≤ 1, where ρ(P), represents the spectral radius of the matrix P.

Numerical examples
This section illustrates the effectiveness of the RBF-PU by means of two numerical examples on the TFSE. To assess the accuracy, we calculate the error norms L ∞ , L 2 and L rms defined as: Example 1 We consider the 2-dim TFSE as The IC and BC can be determined from the exact solution u(x, y, t) = exp(-t) sin(πx) × sin(πy). Table 1 presents the numerical errors L ∞ , L 2 , and L rms and the associated condition number κ(A) when N = 225, M = 25 and T = 1 at various time steps. Table 2 compares the error norms of the LRBF-PU versus the technique described in [26] at different total times when α = 0.5. We can see that the computational accuracy of the RBF-PU is close to [26]. Table 3 presents the values of L ∞ , κ(A) and the CPU running time (in seconds) for a time-step size δt = 1/200 at T = 1. Figure 2 depicts the sparsity structure of the matrix A   for two subdomains, namely M = 81 and M = 49 with N = 400. Figure 3 shows the absolute errors L ∞ at various final times T ∈ {0.1, .5, 1, 2} when α = 0.5.

Example 2 We consider the 2-dim TFSE
The IC and BC can be computed from the exact solution u(x, y, t) = exp(xyt) sin(πx) × sin(πy). Table 4 reports the error norms L ∞ , L 2 , and L rms and the associated condition number κ(A) when N = 225, M = 25 and T = 1 at various time steps. Table 5 compares the error norms of the LRBF-PU with the technique described in [26] at different total times when α = 0.9. We can conclude from Table 5 that the accuracy of the proposed method is slightly   superior to that of [26]. Table 6 illustrates the values of L ∞ , κ(A) and the CPU running time over 1 with a time-step size δt = 1/200 at T = 1. Figure 4 displays the absolute errors L ∞ at various final times T ∈ {0.1, 1, 2, 3, 4, 5} when α = 0.5.

Conclusions
This work adopted an efficient numerical technique following the RBF-PU collocation technique for finding the solution of the TFSE. A major shortcoming of GRBF colloca-  tion methods is the computational burden resulting in dense algebraic systems. The local method deals with the ill-conditioning in GRBF schemes and decreases the associated computational time. The proposed strategy consisted of two stages. First, a FD of order O(δt 2-α ), with the θ -rule, 0 ≤ θ ≤ 1, was implemented to approximate the time dimension. Then, the RBF-PU was provided to discretize the spatial dimension. Numerical experiments illustrated the effectiveness and high accuracy of the RBF-PU.