#### Abstract

In this paper, we propose a numerical method for solving the time fractional Richards’ equation. We first approximate the time fractional derivative of the mentioned equations by a scheme of order *O*(*τ*^{2−α}), 0 < a<1; then, we use the finite point method to approximate the spatial derivatives. Before the discrete spatial derivatives, we introduced the basic principles of the finite point method. We solve the one- and two-dimensional versions of these equations using the proposed method. Moreover, the stability properties of the discretized scheme related to time are theoretically analyzed. Numerical results showed the efficiency of the method presented in this paper.

#### 1. Introduction

In recent years, the theory of fractional calculus and fractional differential equation has been widely used in many fields, such as mechanics, physics, biomathematics, engineering, automatic control, fractal, and so on. Richards’ equation is a basic model for describing soil water movement. The multiscale heterogeneity of soil makes the nature of water diffusion process not consistent with the preconditions of applying Fick’s law, which often reflects the abnormal diffusion phenomenon of soil water infiltration. In response to the above phenomenon, the fractional Richards’ equation is proposed to describe the process of water movement in unsaturated soil [1, 2]. After years of research and practice, numerical simulation has become an effective technical means to simulate soil water flow and infiltration [3–5]. Pachepsky et al. [1] used the finite difference method to solve the time fractional Richard’s equation. Chen et al. [6] used the finite difference method and Kansa method to discretize the time fractional derivatives and the space fractional derivatives, respectively, in solving the time fractional diffusion equation. Freitas et al. [7] proposed the modified fractional integral Richard’s equation to predict the anomalous diffusion process of horizontal infiltration in unsaturated media.

Meshless methods have become very popular in physics and engineering for solving partial differential equations because they do not require mesh reorganization in solving large deformation and many discontinuous problems [8–10]. Meanwhile, a meshless method is not restricted by grid, and its basic idea is to arrange nodes according to certain rules in the solution domain, use shape function to represent unknown field variables in the local region, and finally form stiffness matrix equations to solve it [11]. Onate et al. [12, 13] proposed the finite point method, which first added the stability term to the equation, and then constructed the shape function by moving least squares and discretized the partial differential equation with the stability term by collocation method. Tiwari and Kuhnert [14] developed a finite pointset method, which uses Taylor expansion and weighted least squares to approximate spatial derivatives and uses collocation schemes to discretize equations. Shojaei et al. [15] used the meshless finite point method to solve elastodynamic problems through an explicit velocity. In the same year, Kamranian et al. [16] discussed the two-dimensional initial boundary value problem associated with the sine-Gordon equation using the finite point method. Li and Qin [17] used the meshless finite point method to solve the time fractional convection-diffusion equations, and the numerical results show that this method has higher computational accuracy than the finite difference method. Previously, their work [17] used the finite point method to solve the time fractional linear convection-diffusion equations and obtained better results. Currently, the finite point method is used to solve the time fractional nonlinear soil water movement equation, which belongs to the convection-diffusion equation. The finite point method first applies a stability term to the governing equation and then uses the meshless collocation method to discretize the governing equation. The numerical results show that the calculation accuracy is better when the stability term is applied; that is, applying a stability term can reduce the calculation error and make the simulation result more accurate.

In this paper, we consider the time fractional Richards’ equation as the following: for one-dimensional,subject to the following general initial and boundary conditions:and for two-dimensional,subject to the following general initial and boundary conditions:where , , is soil water content, is diffusivity, and is water conductivity. Assume that the nonlinear term satisfies for positive *m*, *M*. satisfies and , where *L* is the different positive constant.

In equations (1) and (3), is the Caputo fractional derivative of order *α* (0<*α* < 1), which is defined as

Equations (1) and (3) belong to the part of soil water and salt transport model. Aiming at the convection and dispersion phenomenon in the soil water and salt transport model, the finite element method (FEM) and finite difference method (FDM) are easy to produce numerical oscillation when convection is dominant. Numerical methods such as characteristic finite element and characteristic finite difference can eliminate numerical oscillation to a certain extent, but these numerical methods need to rely on the grid in the calculation process, which increases the calculation amount. The meshless method defines the shape function at the nodes in the local domain, which overcomes the dependence on the grid in numerical calculation. Meanwhile, by adding nodes in the computational domain, the computational accuracy in the local region can be improved. The finite point method uses the moving least squares (MLS) method to construct the shape function and the collocation method to discretize the governing equation, which belongs to the meshless method. For equations (1) and (3), the time fractional Caputo derivatives were approximated by *L*1 interpolation, and the space derivatives were discredited by finite point method.

The paper is organized as follows: The finite point method is introduced in Section 2. The scheme of discrete time fractional soil water movement equations is deduced in Section 3. In Section 4, the stability of the time discretization of the one-dimensional and two-dimensional time fractional soil water movement equations is analyzed. In Section 5, some numerical examples are provided to show the accuracy of the proposed method. Finally, a conclusion is drawn in Section 6.

#### 2. Finite Point Method

The finite point method (FPM) uses the moving least squares (MLS) method to construct the shape function and the collocation method to discretize the governing equations. In order to reduce the error and avoid the numerical oscillation, the stability term is added to the discrete equation [12, 13].

##### 2.1. Moving Least Squares (MLS) Approximation

In (0, L), the unknown function *u*(*x*) can be approximated in the neighborhood of the point *x*:where is the vector of MLS shape functions corresponding to *N* nodes in the support domain of the point *x* and can be expressed bywhere is a complete monomials basis of order *m*. *A*(*x*) and *B*(*x*) are defined by

Let , then shape function is . The partial derivatives of the shape functions are, respectively,

From the expression of the shape functions, we can see that the smoothness of the shape function is determined by the weight function, i.e., if , , then . In the present work, Gaussian weight function with compact supports are considered, which is defined aswhere *β* is the shape parameter of the Gaussian function and *r* is the relative distance:where is the distance between the calculation point and the node , and is the support domain size of weight function. The support domain size of node is usually determined by , where scale is the impact factor and is used to control the size of the support domain. in the two-dimensional case.

##### 2.2. Adding Stability Term

The one-dimensional time fractional soil water movement equation is shown in equation (1). LetUsing Taylor’s expansion and flow balance principles, the stability term iswhere *h* is the characteristic length ([12, 13]).

After adding the stability term, equation (1) can be converted into the following form:

Similarly, equation (3) can be converted into the following form:wheregradient ∇*r*the coefficient of convection term , and *h* is the characteristic length.

##### 2.3. Collocation Method

Assume that the equation is as follows:

The Dirichlet boundary condition is

The Neumann boundary condition iswhere *L* and *B* are the differential operators, *u* is the unknown function, *u*_{p} is the prescribed value of *u* in , the boundary , and in Ω.

Assume the approximate expression of *u* is ; then, the approximate function will inevitably produce residuals in the domain and on the boundary, respectively. By using the weighted residual method, equations (18)–(20) are replaced by the following formula:the weight functions *W*_{1}, *W*_{2}, and *W*_{3} have different definitions. Let , where is the Dirac *δ* function, from which the point equation can be obtained aswhere *n*_{u} is the number of nodes on the Dirichlet boundary, *n*_{t} is the number of nodes on the Neumann boundary, *n*_{d} is the number of nodes which satisfy *n*_{d} = *N*-*n*_{t}-*n*_{u} in Ω. The system of equations (22)–(24) leads to a system of algebraic equations of the form

If the differential operators *L* and *B* are nonlinear, it is a nonlinear system of equations, and the approximate value of nodes should be solved by the iterative method [18–21]. Monotone iterative ADI (alternating direction implicit) was used for solving coupled systems of nonlinear parabolic equations, and Boglaev [19] stopped the algorithm and the number of iterations by verifying the convergence order. The inexact Hermitian/Skew-Hermitian Splitting (IHSS) was presented for solving system of nonlinear equations [20]. The stopped criterion for the outer iterations in INHSS algorithm were , where is a prescribed accuracy. A novel biparametric six-order iterative scheme for solving nonlinear systems was presented by Bahl et al. [21]. The number of iterations needed to converge to the solution such that the stopping criterion was satisfied [21].

#### 3. Algorithm Construction of Time Fractional Soil Water Movement Equation

##### 3.1. In the One-Dimensional Case

###### 3.1.1. Time Discretization

Let , where is time step. The time fractional derivative at can be approximated by the following scheme:whereand the truncation error is

So, the *L*1 approximation of Caputo derivative is

Let be the numerical approximation to and ; then, equation (1) can be discretized as the following scheme:

###### 3.1.2. Spatial Discretization

Equation (1) can be transformed into

After adding the stability term to equation (31), we can obtain

The *N* nodes in the domain (a, b) are used for the distribution, and the approximate function at the node *z*_{p} iswhere *φ*_{pi} represents the shape function formed by *z*_{p} as the support domain center and the point *z*_{i}, and *N* is the number of nodes in the local support domain with *z*_{p} as the support domain center.

From equation (33), the first derivative, the second derivative, and the third derivative of the approximation function at the node *z*_{p} can be expressed as

Substititing equations (33) and (34) in equation (32), we can obtain the spatial discrete equation.

Let be the numerical approximation to , from equations (29) and (33)–(35), the fully discrete scheme of equation (1) at iswhere

Let

Equation (35) is as follows:where

##### 3.2. In the Two-Dimensional Case

###### 3.2.1. Time Discretization

Let , where is the time step, the approximation of the time fractional derivative at is shown in equation (29). Let be the numerical approximation to and ; then, equation (3) can be discretized as the following scheme:

###### 3.2.2. Spatial Discretization

The governing equation (3) can be rewritten as

From Section 2.2, it can be seen that the two-dimensional fractional soil-water movement equation after adding the stable term iswhere are shown in Section 2.2. Then, equation (3) can be converted into

The approximate function at iswhere represents the shape function formed by as the center of the support domain and .

The approximate expressions of the unknown function, the first derivative, the second derivative, and the third derivative of the unknown function at can be obtained by constructing the shape function. So, the full discrete scheme of equation (3) at is as follows:

#### 4. Stability

In order to analyze the stability of one-dimensional and two-dimensional time fractional soil-water movement equations, we firstly define the functional spaces endowed with standard norms and inner products:where *L*^{2}(Ω) is the space of measurable functions whose square is Lebesgue integrable in Ω and the inner product form of *L*^{2}(Ω) isand the norm in is

Lemma 1 (see [22]). *Let then*(1)*, when ;*(2)*.*

Lemma 1 gives the magnitude relationship of the coefficients in the time-discrete format.

Lemma 2 (see [23]). *Let { T_{h}}, 0 < h ≤ 1, denote a quasiuniform family of subdivisions of a polyhedral domain . Let be a reference finite element space such that is a finite-dimensional space of functions on is a basis for , where 1 ≤ p ≤ ∞, 1 ≤ q ≤ ∞, and 0 ≤ m ≤ l. For K ∈ T_{h}, let (K, P_{K}, N_{K}) be the affine equivalent element, and V_{h} = {: is measurable and }. Then, there exists C = C (l, p, q) such thatfor all .*

Lemma 3 (see [24]). *Let Ω be a given region with boundary , T_{h} be a uniform partition on Ω, h is the size of unit e in T_{h}, and the finite element space on T_{h}, where P_{k} is a polynomial space whose number of times does not exceed k; then, for the function in the finite element space, there exists a constant C independent of h satisfyingwhere l ≤ k, q ≤ p, and n is the dimension.*

Lemmas 2 and 3 refer to the inverse inequality of finite element, but their forms are different. In the discussion of stability, we use the classical finite element inverse estimation inequality; that is, let be a given region with boundary , *T*_{h} be a uniform partition on , is the size of unit e in *T*_{h}, and the finite element space on *T*_{h}, where *P*_{k} is a polynomial space whose number of times does not exceed *k*; then for the function in the finite element space, there exists a constant *C* independent of *h* satisfying:

Lemma 4 (see [25]). *If V_{h} is a two-dimensional linear finite element space, then when (C_{0} = 12), there iswhere m = 0, p = 2 in , and for any C < C_{0}, there is u_{h} ∈ V_{h} such that [52] does not hold.*

Based on Lemmas 2 and 3, a special case of the finite element inverse estimation inequality is given Lemma 4; that is, when *V*_{h} is a linear finite element space, the constant *C*_{0} in the inverse estimation inequality is a specific number.

Lemma 5 (see [26]). *Given , thanks to the integral form of Gronwall’s inequality, it is derived that*

##### 4.1. In the One-Dimensional Case

The time-discrete scheme of equation (1) is shown in equation (30). When discussing the stability of equation (30), we used the finite element inverse estimation inequality, and the finite element inverse estimation inequality was discussed in the finite element space (see Lemmas 2 and 3). Therefore, we choose the finite element space. From Lemmas 2 and 3, the finite element space can be defined as , where

Theorem 1. *Suppose . is the numerical solution of equation (1); thenwhere C is a positive constant.*

*Proof. *We will prove the above result by mathematical induction.

For *m* = 1, we havewhere .

Multiplying equation (56) by and integrating on Ω, and then using integration by parts, we haveDue to , , and ,From Lemma 1, we haveUsing Schwarz’s inequality and the inverse inequality in Lemmas 2 and 3, we obtainhence,when , thenSuppose that now we have provenMultiplying equation (30) by and integrating on Ω, and then using integration by parts, we haveBy a derivation process similar to *m* = 0, we obtainDue to ,and the proof is finished.

##### 4.2. In the Two-Dimensional Case

The time-discrete scheme of equation (3) is shown in equation (41). When discussing the stability of equation (41), we choose the finite element space, which is defined as equation (56).

Theorem 2. *Suppose . is the numerical solution of equation (3); thenwhere C and ρ are positive constants.*

*Proof. *Multiplying equation (41) by and integrating on Ω, and then using Green’s formula, we haveBy , , , and Schwarz inequality, we obtainFor the derivative part on the right side of the inequality, using the inverse inequality in Lemmas 2 and 3, there isHence,when , we haveAssume , there isi.e.,