# A Finite Point Method for Solving the Time Fractional Richards’ Equation

Authors: Xinqiang Qin; Xin Yang
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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 orderO(τ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.


---

## Body

## 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,(1)0CDtαθ=∂∂zDθ∂θ∂z+∂Kθ∂z+f,z∈Ω1,t>0,subject to the following general initial and boundary conditions:(2)θz,0=θ0z,z∈Ω1,θz,t=φz,t,z∈∂Ω1,t>0,and for two-dimensional,(3)0CDtαθ=∂∂xDθ∂θ∂x+∂∂zDθ∂θ∂z+∂Kθ∂z+f,z∈Ω1,t>0,subject to the following general initial and boundary conditions:(4)θx,z,0=θ0x,z,z∈Ω2,θx,z,t=ψx,z,t,z∈∂Ω2,t>0,where Ω1=0,L, Ω2=0,L×0,L, θ is soil water content, Dθ is diffusivity, and Kθ is water conductivity. Assume that the nonlinear term Dθ satisfies 0<m≤Dθ≤M for positive m, M. Kθ satisfies Kθ≤Lθ and K′θ≤L, where L is the different positive constant.In equations (1) and (3), 0CDtαθ is the Caputo fractional derivative of order α (0<α < 1), which is defined as(5)aCDtαθt=1Γ1−α∫αtθ′st−sαds.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 L1 interpolation, and the space derivatives were discredited by finite point method.The paper is organized as follows: The finite point method is introduced in Section2. 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 functionu(x) can be approximated in the neighborhood Ωx of the point x:(6)ux≈uhx=∑i=1Nφixui=ΦxU,x∈Ω¯,U=u1,u2,…,uNT,where Φx is the vector of MLS shape functions corresponding to N nodes in the support domain of the point x and can be expressed by(7)Φx=pTxA−1xBx,where pTx=p1x,p2x,…,pmx is a complete monomials basis of order m. A(x) and B(x) are defined by(8)Ax=∑i=1NWixpxipTxi,Bx=W1xpx1W1xpx1,…,WNxpxN.Letγ=A−1p, then shape function is ΦT=γTB. The partial derivatives of the shape functions are, respectively,(9)Φ,iT=γ,iTB+γTB,i,Φ,ijT=γ,ijTB+γ,iTB,j+γTB,ij,Φ,ijkT=γ,ijkTB+γ,ijTB,k+γ,ikTB,j+γ,jkTB,i+γ,iTB,jk+γ,jTB,ik+γ,kTB,ij+γTB,ijk.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., ifwx−xi∈CkΩ, i=1,2,…,n, then ϕix∈CkΩ. In the present work, Gaussian weight function with compact supports are considered, which is defined as(10)wx=e−r2β2−e−β21−e−β2,r≤1,0,r>1,where β is the shape parameter of the Gaussian function and r is the relative distance:(11)r=didmi=x−xidmi,where di is the distance between the calculation point x and the node xi, and dmi is the support domain size of weight function. The support domain size dmi of node xi is usually determined by dmi=scale⋅dc, where scale is the impact factor and dc is used to control the size of the support domain. di=x−xi2+y−yi2 in the two-dimensional case.


### 2.2. Adding Stability Term

The one-dimensional time fractional soil water movement equation is shown in equation (1). Let(12)r=∂Kθ∂θ⋅∂θ∂z+Dθ∂2θ∂z2+∂Dθ∂θ∂θ∂z2+f.Using Taylor’s expansion and flow balance principles, the stability term is(13)R=h2⋅∂r∂z,where h is the characteristic length ([12, 13]).After adding the stability term, equation (1) can be converted into the following form:(14)0CDtαθ−r+h2⋅∂r∂z=0.Similarly, equation (3) can be converted into the following form:(15)0CDtαθ−r+h2u⋅u⋅∇r=0,where(16)r=∂Kθ∂θ⋅∂θ∂z+∂Kθ∂θ⋅∂θ∂x+Dθ∂2θ∂x2+∂Dθ∂θ∂θ∂x2+Dθ∂2θ∂z2+∂Dθ∂θ∂θ∂z2+f−∂Kθ∂θ⋅∂θ∂x,gradient ∇r(17)∇r=∂r∂x,∂r∂zT,the coefficient of convection term u=−∂Kθ/∂θ,−∂Kθ/∂θ, and h is the characteristic length.


### 2.3. Collocation Method

Assume that the equation is as follows:(18)Lux=fx,x∈Ω.The Dirichlet boundary condition is(19)u−up=0,x∈Γu.The Neumann boundary condition is(20)Bu=gx,x∈Γt,where L and B are the differential operators, u is the unknown function, up is the prescribed value of u in Γu, the boundary ∂Ω=Γ=Γt∪Γu, and Γt∩Γu=0 in Ω.Assume the approximate expression ofu is u^; 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:(21)∫ΩW1Lu^−fxdΩ+∫ΓtW2Bu^−gxdΓ+∫ΓuW3u^−updΓ,the weight functions W1, W2, and W3 have different definitions. Let W1=W2=W3=δi, where δi is the Dirac δ function, from which the point equation can be obtained as(22)Lu^xi−fxi=0,xi∈Ω,i=1,…,nd,(23)uxi−up=0,xi∈Γu,i=nd+nt+1,…,nd+nt+nu,(24)Bu^xi−gxi=0,xi∈Γt,i=nd+1,…,nd+nt,where nu is the number of nodes on the Dirichlet boundary, nt is the number of nodes on the Neumann boundary, nd is the number of nodes which satisfy nd = N-nt-nu in Ω. The system of equations (22)–(24) leads to a system of algebraic equations of the form(25)Kup=F.If the differential operatorsL 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 Fuk/Fu0≤δ, 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 xk+1−xk+Fxk<10δ was satisfied [21].


## 2.1. Moving Least Squares (MLS) Approximation

In (0, L), the unknown functionu(x) can be approximated in the neighborhood Ωx of the point x:(6)ux≈uhx=∑i=1Nφixui=ΦxU,x∈Ω¯,U=u1,u2,…,uNT,where Φx is the vector of MLS shape functions corresponding to N nodes in the support domain of the point x and can be expressed by(7)Φx=pTxA−1xBx,where pTx=p1x,p2x,…,pmx is a complete monomials basis of order m. A(x) and B(x) are defined by(8)Ax=∑i=1NWixpxipTxi,Bx=W1xpx1W1xpx1,…,WNxpxN.Letγ=A−1p, then shape function is ΦT=γTB. The partial derivatives of the shape functions are, respectively,(9)Φ,iT=γ,iTB+γTB,i,Φ,ijT=γ,ijTB+γ,iTB,j+γTB,ij,Φ,ijkT=γ,ijkTB+γ,ijTB,k+γ,ikTB,j+γ,jkTB,i+γ,iTB,jk+γ,jTB,ik+γ,kTB,ij+γTB,ijk.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., ifwx−xi∈CkΩ, i=1,2,…,n, then ϕix∈CkΩ. In the present work, Gaussian weight function with compact supports are considered, which is defined as(10)wx=e−r2β2−e−β21−e−β2,r≤1,0,r>1,where β is the shape parameter of the Gaussian function and r is the relative distance:(11)r=didmi=x−xidmi,where di is the distance between the calculation point x and the node xi, and dmi is the support domain size of weight function. The support domain size dmi of node xi is usually determined by dmi=scale⋅dc, where scale is the impact factor and dc is used to control the size of the support domain. di=x−xi2+y−yi2 in the two-dimensional case.


## 2.2. Adding Stability Term

The one-dimensional time fractional soil water movement equation is shown in equation (1). Let(12)r=∂Kθ∂θ⋅∂θ∂z+Dθ∂2θ∂z2+∂Dθ∂θ∂θ∂z2+f.Using Taylor’s expansion and flow balance principles, the stability term is(13)R=h2⋅∂r∂z,where h is the characteristic length ([12, 13]).After adding the stability term, equation (1) can be converted into the following form:(14)0CDtαθ−r+h2⋅∂r∂z=0.Similarly, equation (3) can be converted into the following form:(15)0CDtαθ−r+h2u⋅u⋅∇r=0,where(16)r=∂Kθ∂θ⋅∂θ∂z+∂Kθ∂θ⋅∂θ∂x+Dθ∂2θ∂x2+∂Dθ∂θ∂θ∂x2+Dθ∂2θ∂z2+∂Dθ∂θ∂θ∂z2+f−∂Kθ∂θ⋅∂θ∂x,gradient ∇r(17)∇r=∂r∂x,∂r∂zT,the coefficient of convection term u=−∂Kθ/∂θ,−∂Kθ/∂θ, and h is the characteristic length.


## 2.3. Collocation Method

Assume that the equation is as follows:(18)Lux=fx,x∈Ω.The Dirichlet boundary condition is(19)u−up=0,x∈Γu.The Neumann boundary condition is(20)Bu=gx,x∈Γt,where L and B are the differential operators, u is the unknown function, up is the prescribed value of u in Γu, the boundary ∂Ω=Γ=Γt∪Γu, and Γt∩Γu=0 in Ω.Assume the approximate expression ofu is u^; 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:(21)∫ΩW1Lu^−fxdΩ+∫ΓtW2Bu^−gxdΓ+∫ΓuW3u^−updΓ,the weight functions W1, W2, and W3 have different definitions. Let W1=W2=W3=δi, where δi is the Dirac δ function, from which the point equation can be obtained as(22)Lu^xi−fxi=0,xi∈Ω,i=1,…,nd,(23)uxi−up=0,xi∈Γu,i=nd+nt+1,…,nd+nt+nu,(24)Bu^xi−gxi=0,xi∈Γt,i=nd+1,…,nd+nt,where nu is the number of nodes on the Dirichlet boundary, nt is the number of nodes on the Neumann boundary, nd is the number of nodes which satisfy nd = N-nt-nu in Ω. The system of equations (22)–(24) leads to a system of algebraic equations of the form(25)Kup=F.If the differential operatorsL 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 Fuk/Fu0≤δ, 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 xk+1−xk+Fxk<10δ 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

Lettm=mΔt,m=1,2,…,n, where Δt=T/n is time step. The time fractional derivative at t=tm can be approximated by the following scheme:(26)Dtα0Cθtt=tm=1Γ1−α∑k=1m∫tk−1tkθtk−θtk−1Δt⋅1tn−sαds+R˜m=τ−αΓ2−αa0αθtm−∑k=1m−1am−k−1α−am−kαθtk−am−1αθt0+R˜m,where(27)am−k−1α=m−k1−α−m−k−11−α,am−kα=m−k+11−α−m−k1−α,and the truncation error is(28)R˜m≤CΔt2−α.So, theL1 approximation of Caputo derivative is(29)0CDtαθtt=tm≈τ−αΓ2−αa0αθtm−∑k=1m−1am−k−1α−am−kαθtk−am−1αθt0.Letθm=θmx be the numerical approximation to θx,tm and fm=fx,tm; then, equation (1) can be discretized as the following scheme:(30)τ−αΓ2−αa0αθm−∑k=1m−1am−k−1α−am−kαθk−am−1αθ0=∂∂zDθm∂θm∂z+∂Kθm∂z+fm.


#### 3.1.2. Spatial Discretization

Equation (1) can be transformed into(31)0CDtαθ−∂Kθ∂θ⋅∂θ∂z−Dθ∂2θ∂z2−∂Dθ∂θ∂θ∂z2−f=0.After adding the stability term to equation (31), we can obtain(32)0CDtαθ−∂Kθ∂θ⋅∂θ∂z+Dθ∂2θ∂z2+∂Dθ∂θ∂θ∂z2+f−h2−∂2Kθ∂2θ∂θ∂z2+∂Kθ∂θ⋅∂2θ∂z2+3∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂z2+Dθ∂3θ∂z3+∂2Dθ∂θ2∂θ∂z3+∂f∂z=0.TheN nodes z=zii=1,2,…,N in the domain (a, b) are used for the distribution, and the approximate function at the node zp is(33)θzp≅θ^zp=∑i=1Nφpizpθzi,where φpi represents the shape function formed by zp as the support domain center and the point zi, and N is the number of nodes in the local support domain with zp 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 zp can be expressed as(34)∂θz∂zz=zp=∑i=1N∂φpiz∂zθziz=zp,∂2θz∂z2z=zp=∑i=1N∂2φpiz∂z2θziz=zp,∂3θz∂z3z=zp=∑i=1N∂3φpiz∂z3θziz=zp.Substititing equations (33) and (34) in equation (32), we can obtain the spatial discrete equation.Letθpm be the numerical approximation to θzp,tm, from equations (29) and (33)–(35), the fully discrete scheme of equation (1) at θzp,tm is(35)τ−αΓ2−αa0αθpm−∑k=1m−1am−k−1α−am−kαθpk−am−1αθp0=∑i=1Nαpi⋅θim+∂Dθm∂θ−h2∂2Kθm∂θ2∑i=1N∂ϕpi∂z⋅θim2−3h2∂Dθm∂θ∑i=1N∂ϕpi∂z⋅θim∑i=1N∂2ϕpi∂z2⋅θim−h2∂2Dθm∂θ2∑i=1N∂ϕpi∂z⋅θim3−h2∂fm∂zz=zp+fpm,where(36)αpim=∂Kθm∂θ∂φpi∂zz=zp+Dθpm−h2∂Kθm∂θ∂2φpi∂z2z=zp−h2Dθpm∂3φpi∂z3z=zp.Let(37)θm=θ1m,θ2m,…,θpm,…,θNmT,m=1,2,3,….Equation (35) is as follows:(38)Pmθm=τ−αΓ2−α∑k=1m−1am−k−1α−am−kαθk+τ−αΓ2−αam−1αθ0+Fm,where(39)Pm=τ−αΓ2−αa0α⋅E−A−∂Dθm∂θ−h2∂2Kθm∂θ2∑i=1N∂φpi∂z⋅θim⋅B+3h2⋅∂Dθm∂θ∑i=1N∂2φpi∂z2⋅θim⋅C+h2⋅∂2Dθm∂θ2∑i=1N∂φpi∂z⋅θim2⋅B,A=α11α12⋯α1Nα21α22⋯α2N⋮⋮⋮⋮αp1αp2⋯αpN⋮⋮⋮⋮αN1αN2⋯αNN,B=∂φ11∂z∂φ12∂z⋯∂φ1N∂z∂φ21∂z∂φ22∂z⋯∂φ2N∂z⋮⋮⋮⋮∂φp1∂z∂φp2∂z⋯∂φpN∂z⋮⋮⋮⋮∂φN1∂z∂φN2∂z⋯∂φNN∂z,C=∂2φ11∂z2∂2φ12∂z2⋯∂2φ1N∂z2∂2φ21∂z2∂2φ22∂z2⋯∂2φ2N∂z2⋮⋮⋮⋮∂2φp1∂z2∂2φp2∂z2⋯∂2φpN∂z2⋮⋮⋮⋮∂2φN1∂z2∂2φN2∂z2⋯∂2φNN∂z2,F=−h2⋅∂fm∂zz=z1+f1m−h2⋅∂fm∂zz=z2+f2m⋮−h2⋅∂fm∂zz=zp+fpm⋮−h2⋅∂fm∂zz=zN+fNm.


### 3.2. In the Two-Dimensional Case

#### 3.2.1. Time Discretization

Lettm=mΔt,m=1,2,…,n, where Δt=T/n is the time step, the approximation of the time fractional derivative at t=tm is shown in equation (29). Let θm=θmx,z be the numerical approximation to θx,z,tm and fm=fx,z,tm; then, equation (3) can be discretized as the following scheme:(40)τ−αΓ2−αa0αθm−∑k=1m−1am−k−1α−am−kαθk−am−1αθ0=∂∂xDθm∂θm∂x+∂∂zDθm∂θm∂z+∂Kθm∂z+fm.


#### 3.2.2. Spatial Discretization

The governing equation (3) can be rewritten as(41)0CDtαθ−∂Kθ∂θ⋅∂θ∂z−Dθ∂2θ∂x2−∂Dθ∂θ∂θ∂x2−Dθ∂2θ∂z2−∂Dθ∂θ∂θ∂z2−f=0.From Section2.2, it can be seen that the two-dimensional fractional soil-water movement equation after adding the stable term is(42)0CDtαθ−r+h2u⋅u⋅∇r=0,where r,∇r,u are shown in Section 2.2. Then, equation (3) can be converted into(43)0CDtαθ−∂Kθ∂θ⋅∂θ∂z−Dθ∂2θ∂x2−∂Dθ∂θ∂θ∂x2−Dθ∂2θ∂z2−∂Dθ∂θ∂θ∂z2−f+h2u−∂Kθ∂θ∂2Kθ∂θ2⋅∂θ∂x⋅∂θ∂z+∂Kθ∂θ⋅∂2θ∂z∂x+3∂Dθ∂θ⋅∂θ∂x⋅∂2θ∂x2+Dθ∂3θ∂x3+∂2Dθ∂θ2∂θ∂x3+∂Dθ∂θ⋅∂θ∂x⋅∂2θ∂z2+∂2Dθ∂θ2⋅∂θ∂x⋅∂θ∂z2+∂f∂x+Dθ∂3θ∂z2∂x+2∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂z∂x−∂Kθ∂θ∂2Kθ∂θ2∂θ∂z2+∂Kθ∂θ⋅∂2θ∂z2+∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂x2+Dθ∂3θ∂x2∂z+3∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂z2+Dθ∂3θ∂z3+∂2Dθ∂θ2⋅∂θ∂z⋅∂θ∂x2+2∂Dθ∂θ·∂θ∂x⋅∂2θ∂x∂z+∂2Dθ∂θ2⋅∂θ∂z3+∂f∂z=0.The approximate function atθxp,zq is(44)θxp,zq=∑φpq,ijxp,zq⋅θxi,zj,where φpq,ij represents the shape function formed by θxp,zq as the center of the support domain and θxi,zj.The approximate expressions of the unknown function, the first derivative, the second derivative, and the third derivative of the unknown function atθxp,zq can be obtained by constructing the shape function. So, the full discrete scheme of equation (3) at θxp,zq,tm is as follows:(45)τ−αΓ2−αa0αθpqm−∑k=1m−1am−k−1α−am−k−1α−am−1αθpq0=∑ij∂Kθm∂θ∂φpq,ij∂z+Dθpqm∂2φpq,ij∂x2+Dθpqm∂2φpq,ij∂z2+h2u∂Kθm∂θ⋅∂2Kθm∂θ2⋅∂2φpq,ij∂x∂z+h2u∂Kθm∂θ2∂2φpq,ij∂z∂x+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂x3+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂z2∂x+h2u∂Kθm∂θ2∂2φpq,ij∂z2+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂x2∂z+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂z3θ=θxp,zq⋅θijm+∂Dθm∂θ∑ij∂φpq,ij∂xθijm2+∂Dθm∂θ∑ij∂φpq,ij∂zθijm2+h2u∂Kθm∂θ⋅∂2Dθm∂θ2⋅∑ij∂φpq,ij∂zθijm3h2u∂Kθm∂θ3∂Dθm∂θ∑ij∂φpq,ij∂xθijm∑ij∂2φpq,ij∂x2θijm+∂2Dθm∂θ2∑ij∂φpq,ij∂xθijm3+∂Dθm∂θ∑ij∂φpq,ij∂xθijm∑ij∂2φpq,ij∂z2θijm+∂2Dθm∂θ2∑ij∂φpq,ij∂xθijm∑ij∂φpq,ij∂zθijm2+∂Dθm∂θ∑ij∂φpq,ij∂zθijm∑ij∂2φpq,ij∂x2θijm+3∂Dθm∂θ∑ij∂φpq,ij∂zθijm∑ij∂2φpq,ij∂z2θijm+∂2Dθm∂θ2∑ij∂φpq,ij∂zθijm∑ij∂φpq,ij∂xθijm2+2∂Dθm∂θ∑ij∂φpq,ij∂xθijm∑ij∂2φpq,ij∂x∂zθijm+2∂Dθm∂θ∑ij∂φpq,ij∂zθijm∑ij∂φpq,ij∂z∂xθijm+∂2Kθm∂θ2∑ij∂φpq,ij∂Zθijm2+fm+h2u∂Dθm∂θ∂fm∂x+∂fm∂zx=xp,z=zq.


## 3.1. In the One-Dimensional Case

### 3.1.1. Time Discretization

Lettm=mΔt,m=1,2,…,n, where Δt=T/n is time step. The time fractional derivative at t=tm can be approximated by the following scheme:(26)Dtα0Cθtt=tm=1Γ1−α∑k=1m∫tk−1tkθtk−θtk−1Δt⋅1tn−sαds+R˜m=τ−αΓ2−αa0αθtm−∑k=1m−1am−k−1α−am−kαθtk−am−1αθt0+R˜m,where(27)am−k−1α=m−k1−α−m−k−11−α,am−kα=m−k+11−α−m−k1−α,and the truncation error is(28)R˜m≤CΔt2−α.So, theL1 approximation of Caputo derivative is(29)0CDtαθtt=tm≈τ−αΓ2−αa0αθtm−∑k=1m−1am−k−1α−am−kαθtk−am−1αθt0.Letθm=θmx be the numerical approximation to θx,tm and fm=fx,tm; then, equation (1) can be discretized as the following scheme:(30)τ−αΓ2−αa0αθm−∑k=1m−1am−k−1α−am−kαθk−am−1αθ0=∂∂zDθm∂θm∂z+∂Kθm∂z+fm.


### 3.1.2. Spatial Discretization

Equation (1) can be transformed into(31)0CDtαθ−∂Kθ∂θ⋅∂θ∂z−Dθ∂2θ∂z2−∂Dθ∂θ∂θ∂z2−f=0.After adding the stability term to equation (31), we can obtain(32)0CDtαθ−∂Kθ∂θ⋅∂θ∂z+Dθ∂2θ∂z2+∂Dθ∂θ∂θ∂z2+f−h2−∂2Kθ∂2θ∂θ∂z2+∂Kθ∂θ⋅∂2θ∂z2+3∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂z2+Dθ∂3θ∂z3+∂2Dθ∂θ2∂θ∂z3+∂f∂z=0.TheN nodes z=zii=1,2,…,N in the domain (a, b) are used for the distribution, and the approximate function at the node zp is(33)θzp≅θ^zp=∑i=1Nφpizpθzi,where φpi represents the shape function formed by zp as the support domain center and the point zi, and N is the number of nodes in the local support domain with zp 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 zp can be expressed as(34)∂θz∂zz=zp=∑i=1N∂φpiz∂zθziz=zp,∂2θz∂z2z=zp=∑i=1N∂2φpiz∂z2θziz=zp,∂3θz∂z3z=zp=∑i=1N∂3φpiz∂z3θziz=zp.Substititing equations (33) and (34) in equation (32), we can obtain the spatial discrete equation.Letθpm be the numerical approximation to θzp,tm, from equations (29) and (33)–(35), the fully discrete scheme of equation (1) at θzp,tm is(35)τ−αΓ2−αa0αθpm−∑k=1m−1am−k−1α−am−kαθpk−am−1αθp0=∑i=1Nαpi⋅θim+∂Dθm∂θ−h2∂2Kθm∂θ2∑i=1N∂ϕpi∂z⋅θim2−3h2∂Dθm∂θ∑i=1N∂ϕpi∂z⋅θim∑i=1N∂2ϕpi∂z2⋅θim−h2∂2Dθm∂θ2∑i=1N∂ϕpi∂z⋅θim3−h2∂fm∂zz=zp+fpm,where(36)αpim=∂Kθm∂θ∂φpi∂zz=zp+Dθpm−h2∂Kθm∂θ∂2φpi∂z2z=zp−h2Dθpm∂3φpi∂z3z=zp.Let(37)θm=θ1m,θ2m,…,θpm,…,θNmT,m=1,2,3,….Equation (35) is as follows:(38)Pmθm=τ−αΓ2−α∑k=1m−1am−k−1α−am−kαθk+τ−αΓ2−αam−1αθ0+Fm,where(39)Pm=τ−αΓ2−αa0α⋅E−A−∂Dθm∂θ−h2∂2Kθm∂θ2∑i=1N∂φpi∂z⋅θim⋅B+3h2⋅∂Dθm∂θ∑i=1N∂2φpi∂z2⋅θim⋅C+h2⋅∂2Dθm∂θ2∑i=1N∂φpi∂z⋅θim2⋅B,A=α11α12⋯α1Nα21α22⋯α2N⋮⋮⋮⋮αp1αp2⋯αpN⋮⋮⋮⋮αN1αN2⋯αNN,B=∂φ11∂z∂φ12∂z⋯∂φ1N∂z∂φ21∂z∂φ22∂z⋯∂φ2N∂z⋮⋮⋮⋮∂φp1∂z∂φp2∂z⋯∂φpN∂z⋮⋮⋮⋮∂φN1∂z∂φN2∂z⋯∂φNN∂z,C=∂2φ11∂z2∂2φ12∂z2⋯∂2φ1N∂z2∂2φ21∂z2∂2φ22∂z2⋯∂2φ2N∂z2⋮⋮⋮⋮∂2φp1∂z2∂2φp2∂z2⋯∂2φpN∂z2⋮⋮⋮⋮∂2φN1∂z2∂2φN2∂z2⋯∂2φNN∂z2,F=−h2⋅∂fm∂zz=z1+f1m−h2⋅∂fm∂zz=z2+f2m⋮−h2⋅∂fm∂zz=zp+fpm⋮−h2⋅∂fm∂zz=zN+fNm.


## 3.1.1. Time Discretization

Lettm=mΔt,m=1,2,…,n, where Δt=T/n is time step. The time fractional derivative at t=tm can be approximated by the following scheme:(26)Dtα0Cθtt=tm=1Γ1−α∑k=1m∫tk−1tkθtk−θtk−1Δt⋅1tn−sαds+R˜m=τ−αΓ2−αa0αθtm−∑k=1m−1am−k−1α−am−kαθtk−am−1αθt0+R˜m,where(27)am−k−1α=m−k1−α−m−k−11−α,am−kα=m−k+11−α−m−k1−α,and the truncation error is(28)R˜m≤CΔt2−α.So, theL1 approximation of Caputo derivative is(29)0CDtαθtt=tm≈τ−αΓ2−αa0αθtm−∑k=1m−1am−k−1α−am−kαθtk−am−1αθt0.Letθm=θmx be the numerical approximation to θx,tm and fm=fx,tm; then, equation (1) can be discretized as the following scheme:(30)τ−αΓ2−αa0αθm−∑k=1m−1am−k−1α−am−kαθk−am−1αθ0=∂∂zDθm∂θm∂z+∂Kθm∂z+fm.


## 3.1.2. Spatial Discretization

Equation (1) can be transformed into(31)0CDtαθ−∂Kθ∂θ⋅∂θ∂z−Dθ∂2θ∂z2−∂Dθ∂θ∂θ∂z2−f=0.After adding the stability term to equation (31), we can obtain(32)0CDtαθ−∂Kθ∂θ⋅∂θ∂z+Dθ∂2θ∂z2+∂Dθ∂θ∂θ∂z2+f−h2−∂2Kθ∂2θ∂θ∂z2+∂Kθ∂θ⋅∂2θ∂z2+3∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂z2+Dθ∂3θ∂z3+∂2Dθ∂θ2∂θ∂z3+∂f∂z=0.TheN nodes z=zii=1,2,…,N in the domain (a, b) are used for the distribution, and the approximate function at the node zp is(33)θzp≅θ^zp=∑i=1Nφpizpθzi,where φpi represents the shape function formed by zp as the support domain center and the point zi, and N is the number of nodes in the local support domain with zp 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 zp can be expressed as(34)∂θz∂zz=zp=∑i=1N∂φpiz∂zθziz=zp,∂2θz∂z2z=zp=∑i=1N∂2φpiz∂z2θziz=zp,∂3θz∂z3z=zp=∑i=1N∂3φpiz∂z3θziz=zp.Substititing equations (33) and (34) in equation (32), we can obtain the spatial discrete equation.Letθpm be the numerical approximation to θzp,tm, from equations (29) and (33)–(35), the fully discrete scheme of equation (1) at θzp,tm is(35)τ−αΓ2−αa0αθpm−∑k=1m−1am−k−1α−am−kαθpk−am−1αθp0=∑i=1Nαpi⋅θim+∂Dθm∂θ−h2∂2Kθm∂θ2∑i=1N∂ϕpi∂z⋅θim2−3h2∂Dθm∂θ∑i=1N∂ϕpi∂z⋅θim∑i=1N∂2ϕpi∂z2⋅θim−h2∂2Dθm∂θ2∑i=1N∂ϕpi∂z⋅θim3−h2∂fm∂zz=zp+fpm,where(36)αpim=∂Kθm∂θ∂φpi∂zz=zp+Dθpm−h2∂Kθm∂θ∂2φpi∂z2z=zp−h2Dθpm∂3φpi∂z3z=zp.Let(37)θm=θ1m,θ2m,…,θpm,…,θNmT,m=1,2,3,….Equation (35) is as follows:(38)Pmθm=τ−αΓ2−α∑k=1m−1am−k−1α−am−kαθk+τ−αΓ2−αam−1αθ0+Fm,where(39)Pm=τ−αΓ2−αa0α⋅E−A−∂Dθm∂θ−h2∂2Kθm∂θ2∑i=1N∂φpi∂z⋅θim⋅B+3h2⋅∂Dθm∂θ∑i=1N∂2φpi∂z2⋅θim⋅C+h2⋅∂2Dθm∂θ2∑i=1N∂φpi∂z⋅θim2⋅B,A=α11α12⋯α1Nα21α22⋯α2N⋮⋮⋮⋮αp1αp2⋯αpN⋮⋮⋮⋮αN1αN2⋯αNN,B=∂φ11∂z∂φ12∂z⋯∂φ1N∂z∂φ21∂z∂φ22∂z⋯∂φ2N∂z⋮⋮⋮⋮∂φp1∂z∂φp2∂z⋯∂φpN∂z⋮⋮⋮⋮∂φN1∂z∂φN2∂z⋯∂φNN∂z,C=∂2φ11∂z2∂2φ12∂z2⋯∂2φ1N∂z2∂2φ21∂z2∂2φ22∂z2⋯∂2φ2N∂z2⋮⋮⋮⋮∂2φp1∂z2∂2φp2∂z2⋯∂2φpN∂z2⋮⋮⋮⋮∂2φN1∂z2∂2φN2∂z2⋯∂2φNN∂z2,F=−h2⋅∂fm∂zz=z1+f1m−h2⋅∂fm∂zz=z2+f2m⋮−h2⋅∂fm∂zz=zp+fpm⋮−h2⋅∂fm∂zz=zN+fNm.


## 3.2. In the Two-Dimensional Case

### 3.2.1. Time Discretization

Lettm=mΔt,m=1,2,…,n, where Δt=T/n is the time step, the approximation of the time fractional derivative at t=tm is shown in equation (29). Let θm=θmx,z be the numerical approximation to θx,z,tm and fm=fx,z,tm; then, equation (3) can be discretized as the following scheme:(40)τ−αΓ2−αa0αθm−∑k=1m−1am−k−1α−am−kαθk−am−1αθ0=∂∂xDθm∂θm∂x+∂∂zDθm∂θm∂z+∂Kθm∂z+fm.


### 3.2.2. Spatial Discretization

The governing equation (3) can be rewritten as(41)0CDtαθ−∂Kθ∂θ⋅∂θ∂z−Dθ∂2θ∂x2−∂Dθ∂θ∂θ∂x2−Dθ∂2θ∂z2−∂Dθ∂θ∂θ∂z2−f=0.From Section2.2, it can be seen that the two-dimensional fractional soil-water movement equation after adding the stable term is(42)0CDtαθ−r+h2u⋅u⋅∇r=0,where r,∇r,u are shown in Section 2.2. Then, equation (3) can be converted into(43)0CDtαθ−∂Kθ∂θ⋅∂θ∂z−Dθ∂2θ∂x2−∂Dθ∂θ∂θ∂x2−Dθ∂2θ∂z2−∂Dθ∂θ∂θ∂z2−f+h2u−∂Kθ∂θ∂2Kθ∂θ2⋅∂θ∂x⋅∂θ∂z+∂Kθ∂θ⋅∂2θ∂z∂x+3∂Dθ∂θ⋅∂θ∂x⋅∂2θ∂x2+Dθ∂3θ∂x3+∂2Dθ∂θ2∂θ∂x3+∂Dθ∂θ⋅∂θ∂x⋅∂2θ∂z2+∂2Dθ∂θ2⋅∂θ∂x⋅∂θ∂z2+∂f∂x+Dθ∂3θ∂z2∂x+2∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂z∂x−∂Kθ∂θ∂2Kθ∂θ2∂θ∂z2+∂Kθ∂θ⋅∂2θ∂z2+∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂x2+Dθ∂3θ∂x2∂z+3∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂z2+Dθ∂3θ∂z3+∂2Dθ∂θ2⋅∂θ∂z⋅∂θ∂x2+2∂Dθ∂θ·∂θ∂x⋅∂2θ∂x∂z+∂2Dθ∂θ2⋅∂θ∂z3+∂f∂z=0.The approximate function atθxp,zq is(44)θxp,zq=∑φpq,ijxp,zq⋅θxi,zj,where φpq,ij represents the shape function formed by θxp,zq as the center of the support domain and θxi,zj.The approximate expressions of the unknown function, the first derivative, the second derivative, and the third derivative of the unknown function atθxp,zq can be obtained by constructing the shape function. So, the full discrete scheme of equation (3) at θxp,zq,tm is as follows:(45)τ−αΓ2−αa0αθpqm−∑k=1m−1am−k−1α−am−k−1α−am−1αθpq0=∑ij∂Kθm∂θ∂φpq,ij∂z+Dθpqm∂2φpq,ij∂x2+Dθpqm∂2φpq,ij∂z2+h2u∂Kθm∂θ⋅∂2Kθm∂θ2⋅∂2φpq,ij∂x∂z+h2u∂Kθm∂θ2∂2φpq,ij∂z∂x+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂x3+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂z2∂x+h2u∂Kθm∂θ2∂2φpq,ij∂z2+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂x2∂z+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂z3θ=θxp,zq⋅θijm+∂Dθm∂θ∑ij∂φpq,ij∂xθijm2+∂Dθm∂θ∑ij∂φpq,ij∂zθijm2+h2u∂Kθm∂θ⋅∂2Dθm∂θ2⋅∑ij∂φpq,ij∂zθijm3h2u∂Kθm∂θ3∂Dθm∂θ∑ij∂φpq,ij∂xθijm∑ij∂2φpq,ij∂x2θijm+∂2Dθm∂θ2∑ij∂φpq,ij∂xθijm3+∂Dθm∂θ∑ij∂φpq,ij∂xθijm∑ij∂2φpq,ij∂z2θijm+∂2Dθm∂θ2∑ij∂φpq,ij∂xθijm∑ij∂φpq,ij∂zθijm2+∂Dθm∂θ∑ij∂φpq,ij∂zθijm∑ij∂2φpq,ij∂x2θijm+3∂Dθm∂θ∑ij∂φpq,ij∂zθijm∑ij∂2φpq,ij∂z2θijm+∂2Dθm∂θ2∑ij∂φpq,ij∂zθijm∑ij∂φpq,ij∂xθijm2+2∂Dθm∂θ∑ij∂φpq,ij∂xθijm∑ij∂2φpq,ij∂x∂zθijm+2∂Dθm∂θ∑ij∂φpq,ij∂zθijm∑ij∂φpq,ij∂z∂xθijm+∂2Kθm∂θ2∑ij∂φpq,ij∂Zθijm2+fm+h2u∂Dθm∂θ∂fm∂x+∂fm∂zx=xp,z=zq.


## 3.2.1. Time Discretization

Lettm=mΔt,m=1,2,…,n, where Δt=T/n is the time step, the approximation of the time fractional derivative at t=tm is shown in equation (29). Let θm=θmx,z be the numerical approximation to θx,z,tm and fm=fx,z,tm; then, equation (3) can be discretized as the following scheme:(40)τ−αΓ2−αa0αθm−∑k=1m−1am−k−1α−am−kαθk−am−1αθ0=∂∂xDθm∂θm∂x+∂∂zDθm∂θm∂z+∂Kθm∂z+fm.


## 3.2.2. Spatial Discretization

The governing equation (3) can be rewritten as(41)0CDtαθ−∂Kθ∂θ⋅∂θ∂z−Dθ∂2θ∂x2−∂Dθ∂θ∂θ∂x2−Dθ∂2θ∂z2−∂Dθ∂θ∂θ∂z2−f=0.From Section2.2, it can be seen that the two-dimensional fractional soil-water movement equation after adding the stable term is(42)0CDtαθ−r+h2u⋅u⋅∇r=0,where r,∇r,u are shown in Section 2.2. Then, equation (3) can be converted into(43)0CDtαθ−∂Kθ∂θ⋅∂θ∂z−Dθ∂2θ∂x2−∂Dθ∂θ∂θ∂x2−Dθ∂2θ∂z2−∂Dθ∂θ∂θ∂z2−f+h2u−∂Kθ∂θ∂2Kθ∂θ2⋅∂θ∂x⋅∂θ∂z+∂Kθ∂θ⋅∂2θ∂z∂x+3∂Dθ∂θ⋅∂θ∂x⋅∂2θ∂x2+Dθ∂3θ∂x3+∂2Dθ∂θ2∂θ∂x3+∂Dθ∂θ⋅∂θ∂x⋅∂2θ∂z2+∂2Dθ∂θ2⋅∂θ∂x⋅∂θ∂z2+∂f∂x+Dθ∂3θ∂z2∂x+2∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂z∂x−∂Kθ∂θ∂2Kθ∂θ2∂θ∂z2+∂Kθ∂θ⋅∂2θ∂z2+∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂x2+Dθ∂3θ∂x2∂z+3∂Dθ∂θ⋅∂θ∂z⋅∂2θ∂z2+Dθ∂3θ∂z3+∂2Dθ∂θ2⋅∂θ∂z⋅∂θ∂x2+2∂Dθ∂θ·∂θ∂x⋅∂2θ∂x∂z+∂2Dθ∂θ2⋅∂θ∂z3+∂f∂z=0.The approximate function atθxp,zq is(44)θxp,zq=∑φpq,ijxp,zq⋅θxi,zj,where φpq,ij represents the shape function formed by θxp,zq as the center of the support domain and θxi,zj.The approximate expressions of the unknown function, the first derivative, the second derivative, and the third derivative of the unknown function atθxp,zq can be obtained by constructing the shape function. So, the full discrete scheme of equation (3) at θxp,zq,tm is as follows:(45)τ−αΓ2−αa0αθpqm−∑k=1m−1am−k−1α−am−k−1α−am−1αθpq0=∑ij∂Kθm∂θ∂φpq,ij∂z+Dθpqm∂2φpq,ij∂x2+Dθpqm∂2φpq,ij∂z2+h2u∂Kθm∂θ⋅∂2Kθm∂θ2⋅∂2φpq,ij∂x∂z+h2u∂Kθm∂θ2∂2φpq,ij∂z∂x+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂x3+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂z2∂x+h2u∂Kθm∂θ2∂2φpq,ij∂z2+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂x2∂z+h2u∂Kθm∂θ⋅Dθm⋅∂3φpq,ij∂z3θ=θxp,zq⋅θijm+∂Dθm∂θ∑ij∂φpq,ij∂xθijm2+∂Dθm∂θ∑ij∂φpq,ij∂zθijm2+h2u∂Kθm∂θ⋅∂2Dθm∂θ2⋅∑ij∂φpq,ij∂zθijm3h2u∂Kθm∂θ3∂Dθm∂θ∑ij∂φpq,ij∂xθijm∑ij∂2φpq,ij∂x2θijm+∂2Dθm∂θ2∑ij∂φpq,ij∂xθijm3+∂Dθm∂θ∑ij∂φpq,ij∂xθijm∑ij∂2φpq,ij∂z2θijm+∂2Dθm∂θ2∑ij∂φpq,ij∂xθijm∑ij∂φpq,ij∂zθijm2+∂Dθm∂θ∑ij∂φpq,ij∂zθijm∑ij∂2φpq,ij∂x2θijm+3∂Dθm∂θ∑ij∂φpq,ij∂zθijm∑ij∂2φpq,ij∂z2θijm+∂2Dθm∂θ2∑ij∂φpq,ij∂zθijm∑ij∂φpq,ij∂xθijm2+2∂Dθm∂θ∑ij∂φpq,ij∂xθijm∑ij∂2φpq,ij∂x∂zθijm+2∂Dθm∂θ∑ij∂φpq,ij∂zθijm∑ij∂φpq,ij∂z∂xθijm+∂2Kθm∂θ2∑ij∂φpq,ij∂Zθijm2+fm+h2u∂Dθm∂θ∂fm∂x+∂fm∂zx=xp,z=zq.


## 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:(46)H1Ω=v∈L2Ω,∇v∈L2Ω,H01Ω=v∈H1Ω,v∂Ω=0,where L2(Ω) is the space of measurable functions whose square is Lebesgue integrable in Ω and the inner product form of L2(Ω) is(47)u,v=∫Ωuxvxdx,u,v∈L2Ω,u,v=∫∫Ωuxvxdxdy,u,v∈L2Ω,and the norm in L2Ω is(48)θ=θ,θ1/2=∫Ωθ2dx1/2,θ=θ,θ1/2=∫∫Ωθ2dxdy1/2.Lemma 1 (see [22]).
Letα∈0,1,alα=l+11−α−l1−α,l=0,1,2,…, then(1)
1=a0α>a1α>a2α>⋯>alα>0;alα⟶0, when l⟶∞;(2)
1−αl−α<al−1α<1−αl−1−α,l≥1.Lemma 1 gives the magnitude relationship of the coefficients in the time-discrete format.Lemma 2 (see [23]).
Let {Th}, 0 < h ≤ 1, denote a quasiuniform family of subdivisions of a polyhedral domain Ω⊂Rd′. Let K^,P,N be a reference finite element space such that P⊂Wl,pK^∩Wm,qK^ is a finite-dimensional space of functions on K^,N is a basis for P′, where 1 ≤ p ≤ ∞, 1 ≤ q ≤ ∞, and 0 ≤ m ≤ l. For K ∈ Th, let (K, PK, NK) be the affine equivalent element, and Vh = {v: v is measurable and v∣K∈PK,∀K∈Th}. Then, there exists C = C (l, p, q) such that(49)∑K∈ThvWi,pKp1/p≤Chm−l+min0,d′/p−d′/q∑K∈ThvWm,qKq1/q.for all v∈Vh.Lemma 3 (see [24]).
Let Ω be a given region with boundary∂Ω, Th be a uniform partition on Ω, h is the size of unit e in Th, and the finite element space Vh=v:v∣e∈Pke,∀e∈Th on Th, where Pk is a polynomial space whose number of times does not exceed k; then, for the function uh∈Vh in the finite element space, there exists a constant C independent of h satisfying(50)uhk,p,Ω≤Chl−k+n1/p−1/quhl,q,Ω,where 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 ∂Ω, Th be a uniform partition on Ω, h is the size of unit e in Th, and the finite element space Vh=v:v∣e∈Pke,∀e∈Th on Th, where Pk is a polynomial space whose number of times does not exceed k; then for the function uh∈Vh in the finite element space, there exists a constant C independent of h satisfying:(51)∇uh≤Ch−1uh.Lemma 4 (see [25]).
IfVh is a two-dimensional linear finite element space, then when C0=12 (C0 = 12), there is(52)∇uh0,2≤C0h−1uh0,2,where m = 0, p = 2 in um,p,Ω, and for any C < C0, there is uh ∈ Vh 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, whenVh is a linear finite element space, the constant C0 in the inverse estimation inequality is a specific number.Lemma 5 (see [26]).
GivenEt≤ρt+C∫0TEsds, thanks to the integral form of Gronwall’s inequality, it is derived that(53)Et≤ρteCT,0<t≤T.


### 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 Vh⊂H01, where(54)Vh=v∈H01Ω∩C0Ω¯∣v∣e∈Pke,∀e∈Th.Theorem 1.
Supposeθm∈Vh⊂H01Ω,m=0,1,2,…,n. θm is the numerical solution of equation (1); then(55)θm≤Cθ0+λCmax0≤l≤nfl,where C is a positive constant.Proof.
We will prove the above result by mathematical induction.
Form = 1, we have(56)a0αθ1−a0αθ0=λ∂∂zDθ1∂θ1∂z+λ∂Kθ1∂z+λf1,where λ=ταΓ2−α.
Multiplying equation (56) by θ1 and integrating on Ω, and then using integration by parts, we have(57)a0αθ1,θ1=a0θ0,θ1−λDθ1,∂θ1∂z2−λKθ1,∂θ1∂z+λf1,θ1.
Due to0<m≤Dθ≤M, Kθ≤Lθ, and K′θ≤L,(58)θ12≤a0θ0,θ1+λM,∂θ1∂z2+λf1,θ1.
From Lemma 1, we have(59)θ12≤θ0,θ1+λM,∂θ1∂z2+λf1,θ1.
Using Schwarz’s inequality and the inverse inequality in Lemmas 2 and 3, we obtain(60)θ12≤θ0θ1+λM∂θ1∂z2+λf1θ1≤θ0θ1+λMch−2θ12+λf1θ1,hence,(61)1−λMch−2θ1≤θ0+λmax0≤l≤nfl,when 1−λMch−2>0, then(62)θ1≤Cθ0+λCmax0≤l≤nfl.
Suppose that now we have proven(63)θk≤Cθ0+λCmax0≤l≤Nfl,k=1,2,3,…,m−1.
Multiplying equation (30) by θm and integrating on Ω, and then using integration by parts, we have(64)θm,θm=∑k=1m−1am−k−1α−am−kαθk,θm+am−1αθ0,θm−λDθm,∂θm∂z2−λKθm,∂θm∂z2+λfm,θm.
By a derivation process similar tom = 0, we obtain(65)1−λMch−2θm≤∑k=1m−1am−k−1α−am−kαθk+am−1αθ0+λmax0≤l≤nfl≤∑k=1m−1am−k−1α−am−kαCθ0+λCmax0≤l≤nfl+am−1αθ0+λmax0≤l≤nfl≤Ca0α−am−1α+a0αθ0+λCa0α−am−1α+a0αmax0≤l≤nfl≤Cθ0+λCmax0≤l≤nfl.
Due to1−λMch−2>0,(66)θm≤Cθ0+λCmax0≤l≤nfl,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θm∈Vh⊂H01Ω,m=0,1,2,…,n. θm is the numerical solution of equation (3); then(67)θm≤ρeCT,where C and ρ are positive constants.Proof.
Multiplying equation (41) by θm and integrating on Ω, and then using Green’s formula, we have(68)θm,θm=∑k=1m−1am−k−1α−am−kαθk,θm+anαθ0,θm−λDθm∂θm∂x,∂θm∂x−λDθm∂θm∂z,∂θm∂z+λ∂Kθm∂z,θm+λfm,θm.
By0<m≤Dθ≤M, Kθ≤Lθ, K′θ≤L, and Schwarz inequality, we obtain(69)θm2+2λM∂θm∂x∂θm∂z≤∑k=1m−1am−k−1α−am−kαθkθm+am−1αθ0θm+λM∂θm∂x2+λM∂θm∂z2+2λM∂θm∂x∂θm∂z+λL∂θm∂zθm+λfmθm.
For the derivative part on the right side of the inequality, using the inverse inequality in Lemmas 2 and 3, there is(70)λM∂θm∂x2+2λM∂θm∂x∂θm∂z+λM∂θm∂z2+λL∂θm∂zθm≤λMc1h−2θm2+λLc2h−1θm2.
Hence,(71)1−λMc1h−2−λLc2h−1θm≤∑k=1m−1am−k−1α−am−kαθk+am−1αθ0+λfm,when 1−λMc1h−2−λLc2h−1≥0, we have(72)θm≤11−λMc1h−2−λLc2h−1∑k=1m−1am−k−1α−am−kαθk+am−1α1−λMc1h−2−λLc2h−1θ0+λ1−λMc1h−2−λLc2h−1fm.
AssumeC2=max1≤k≤m−1am−k−1α−am−kα, there is(73)θm≤C21−λMc1h−2−λLc2h−1∑k=1m−1θk+11−λMc1h−2−λLc2h−1am−1αθ0+λfm=C∑k=1m−1θk+Dam−1αθ0+λfm=C∑k=1m−1θk+ρ.i.e.,(74)θm≤ρ+C∑k=1m−1θk.
Applying Lemma 5, we obtain(75)θm≤ρeCT.


## 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 Vh⊂H01, where(54)Vh=v∈H01Ω∩C0Ω¯∣v∣e∈Pke,∀e∈Th.Theorem 1.
Supposeθm∈Vh⊂H01Ω,m=0,1,2,…,n. θm is the numerical solution of equation (1); then(55)θm≤Cθ0+λCmax0≤l≤nfl,where C is a positive constant.Proof.
We will prove the above result by mathematical induction.
Form = 1, we have(56)a0αθ1−a0αθ0=λ∂∂zDθ1∂θ1∂z+λ∂Kθ1∂z+λf1,where λ=ταΓ2−α.
Multiplying equation (56) by θ1 and integrating on Ω, and then using integration by parts, we have(57)a0αθ1,θ1=a0θ0,θ1−λDθ1,∂θ1∂z2−λKθ1,∂θ1∂z+λf1,θ1.
Due to0<m≤Dθ≤M, Kθ≤Lθ, and K′θ≤L,(58)θ12≤a0θ0,θ1+λM,∂θ1∂z2+λf1,θ1.
From Lemma 1, we have(59)θ12≤θ0,θ1+λM,∂θ1∂z2+λf1,θ1.
Using Schwarz’s inequality and the inverse inequality in Lemmas 2 and 3, we obtain(60)θ12≤θ0θ1+λM∂θ1∂z2+λf1θ1≤θ0θ1+λMch−2θ12+λf1θ1,hence,(61)1−λMch−2θ1≤θ0+λmax0≤l≤nfl,when 1−λMch−2>0, then(62)θ1≤Cθ0+λCmax0≤l≤nfl.
Suppose that now we have proven(63)θk≤Cθ0+λCmax0≤l≤Nfl,k=1,2,3,…,m−1.
Multiplying equation (30) by θm and integrating on Ω, and then using integration by parts, we have(64)θm,θm=∑k=1m−1am−k−1α−am−kαθk,θm+am−1αθ0,θm−λDθm,∂θm∂z2−λKθm,∂θm∂z2+λfm,θm.
By a derivation process similar tom = 0, we obtain(65)1−λMch−2θm≤∑k=1m−1am−k−1α−am−kαθk+am−1αθ0+λmax0≤l≤nfl≤∑k=1m−1am−k−1α−am−kαCθ0+λCmax0≤l≤nfl+am−1αθ0+λmax0≤l≤nfl≤Ca0α−am−1α+a0αθ0+λCa0α−am−1α+a0αmax0≤l≤nfl≤Cθ0+λCmax0≤l≤nfl.
Due to1−λMch−2>0,(66)θm≤Cθ0+λCmax0≤l≤nfl,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θm∈Vh⊂H01Ω,m=0,1,2,…,n. θm is the numerical solution of equation (3); then(67)θm≤ρeCT,where C and ρ are positive constants.Proof.
Multiplying equation (41) by θm and integrating on Ω, and then using Green’s formula, we have(68)θm,θm=∑k=1m−1am−k−1α−am−kαθk,θm+anαθ0,θm−λDθm∂θm∂x,∂θm∂x−λDθm∂θm∂z,∂θm∂z+λ∂Kθm∂z,θm+λfm,θm.
By0<m≤Dθ≤M, Kθ≤Lθ, K′θ≤L, and Schwarz inequality, we obtain(69)θm2+2λM∂θm∂x∂θm∂z≤∑k=1m−1am−k−1α−am−kαθkθm+am−1αθ0θm+λM∂θm∂x2+λM∂θm∂z2+2λM∂θm∂x∂θm∂z+λL∂θm∂zθm+λfmθm.
For the derivative part on the right side of the inequality, using the inverse inequality in Lemmas 2 and 3, there is(70)λM∂θm∂x2+2λM∂θm∂x∂θm∂z+λM∂θm∂z2+λL∂θm∂zθm≤λMc1h−2θm2+λLc2h−1θm2.
Hence,(71)1−λMc1h−2−λLc2h−1θm≤∑k=1m−1am−k−1α−am−kαθk+am−1αθ0+λfm,when 1−λMc1h−2−λLc2h−1≥0, we have(72)θm≤11−λMc1h−2−λLc2h−1∑k=1m−1am−k−1α−am−kαθk+am−1α1−λMc1h−2−λLc2h−1θ0+λ1−λMc1h−2−λLc2h−1fm.
AssumeC2=max1≤k≤m−1am−k−1α−am−kα, there is(73)θm≤C21−λMc1h−2−λLc2h−1∑k=1m−1θk+11−λMc1h−2−λLc2h−1am−1αθ0+λfm=C∑k=1m−1θk+Dam−1αθ0+λfm=C∑k=1m−1θk+ρ.i.e.,(74)θm≤ρ+C∑k=1m−1θk.
Applying Lemma 5, we obtain(75)θm≤ρeCT.


## 5. Results and Discussion

In this section, we show the results for two examples using the method described in previous sections. To show the accuracy of the proposed method, maximum absolute error and relative error are presented as(76)ε1=max1≤i≤Nθzi−θexactzi,ε2=∑i=1Nθzi−θexactzi2∑i=1Nθzi2,where θzi and θexactzi for i=1,2,…,N are computed by the exact and numerical solutions on points zi, respectively, and N is the number of nodal points.Gaussian function has arbitrary order derivative. We choose Gaussian function as the weight function, and the basis function is quadratic basis functionpΤ=1,x,x2 in the example.Example 1.
Consider one-dimensional time fractional soil water movement equation(77)0CDtαθ=∂∂zDθ∂θ∂z+∂Kθ∂z+f,z,t∈0,1×0,T,θz,0=0,z∈0,1,θ0,t=0,t∈0,T,θ1,t=0,t∈0,T,where Kθ=θ,Dθ=0.001θ+0.001, the exact solution is θ=t2z−z2, and fx can be obtained by substituting the exact solution into the equation.
We take ∆t = 10−3, the basis function is quadratic basis function pT=1,x,x2 and scale is 2.5 in Example 1.
Figure1 shows the comparison between the FPM and the exact solution at different times when α = 0.5 and α = 0.8, respectively. It can be seen from Figure 1 that the finite point numerical solution and the exact solution can still agree well even when T is large.
Table1 shows the ε2 errors of the FPM and the FDM under different α when T = 1, ∆t = 10−3, and the number of nodes is 11. It can be seen from Table 1 that as α increases, the ε2 errors of both the FPM and the FDM increase gradually, but the ε2 errors of the FPM is significantly smaller than the FDM.
Whenα = 0.8, Table 2 gives the ε2 errors of the FPM and the FDM under different node numbers at T = 1. Table 3 gives the ε2 errors of the FPM and the FDM at different times. Table 4 verifies the convergence order of the FPM. As can be seen from Table 2, with the increase of the number of spatial nodes, the numerical solution is closer to the exact solution, and the ε2 errors of the FPM are still higher than those of the FDM. It can be seen from Table 3 that the ε2 errors of the FPM and the FDM decrease gradually with the increase of time, but the results of proposed method are still better than those of the FDM and the ε2 errors can still reach the order of e − 005 when T = 2, which shows that the proposed method can obtain relatively stable results in simulating water infiltration process. From Table 4, it can be seen that the convergence order is closer to 1.2 as the time step decreases when α = 0.8, which shows that the convergence order is closer to the theoretical convergence order.Figure 1
Comparing the proposed method with the exact solution when ∆t = 10−3 and the number of nodes is 21, where α = 0.5 (a), α = 0.8 (b).
(a)
(b)Table 1
Comparison ofε2 errors between the FPM and the FDM under different α, when T = 1, ∆t = 10−3, and the number of nodes is 11.
α
0.4
0.5
0.6
0.7
0.8
FPM
2.5156e − 06
6.9681e − 06
1.8722e − 05
4.7683e − 05
1.1539e − 04
FDM
1.4281e − 03
1.4360e − 03
1.4519e − 03
1.4823e − 03
1.5359e − 03Table 2
Comparison ofε2 errors between FPM and FDM under different node numbers, when α = 0.8.
Time step
Time
Number of nodes
FPM
FDM
0.001
1
11
1.1539e − 04
1.5339e − 03
21
1.0193e − 04
1.4928e − 03
41
9.6788e − 05
1.4326e − 03
81
9.5715e − 05
1.4027e − 03Table 3
Comparison ofε2 errors between FPM and FDM under different times, when α = 0.8.
Time step
Number of nodes
Time
FPM
FDM
0.001
21
0.1
3.1502e − 03
5.5855e − 03
0.5
3.4321e − 04
2.5670e − 03
1.0
1.0193e − 04
1.4928e − 03
2.0
2.7596e − 05
8.0918e − 04Table 4
ε1 errors and convergence order of numerical solution, when α = 0.8 and the number of nodes is 11.
Time step
ε1
Convergence rate
1/20
3.3612e − 03
1/40
1.5831e − 03
1.0862
1/80
7.4187e − 04
1.0935
1/160
3.4394e − 04
1.1090
1/320
1.5711e − 04
1.1304
1/640
7.0722e − 05
1.1515Example 2.
Consider one-dimensional time fractional soil water movement equation:(78)0CDtαθ=∂∂xDθ∂θ∂x+∂∂zDθ∂θ∂z+∂Kθ∂z+f,x,z,t∈0,1×0,1×0,T,θx,z,0=x−x2z−z2,x,z∈0,1×0,L,θ0,z,t=0,z,t∈0,1×0,T,θ1,z,t=0,z,t∈0,1×0,T,θx,0,t=0,z,t∈0,1×0,T,θx,1,t=0,z,t∈0,1×0,T,where Kθ=θ,Dθ=0.001θ+0.001, the exact solution is θ=t2x−x2z−z2, and f(x) can be obtained by substituting the exact solution into the equation.
We take ∆t is 10−3, the basis function is quadratic basis function pΤ=1,x,x2, and scale is 2 in Example 2.
WhenT = 1 and α = 0.8, Figure 2 gives the comparison between the numerical solution of FPM and the exact solution when the number of nodes is 21 × 21. Figure 3 gives the equipotentia of the exact solution and the numerical solution when the number of nodes is 31 × 31. From Figure 3, we can see that the numerical solution of FPM is basically consistent with the exact solution.
Table5 shows the ε2 errors of the FPM and the FDM under different α when T = 1, ∆t = 10−3, and the number of nodes is 21 × 21. It can be seen from Table 5 that as α increases, the ε2 errors of the FPM is significantly smaller than the FDM. when α = 0.8. Table 6 gives the ε2 errors of the FPM and the FDM under different number of nodes. From Table 6, we can see that the ε2 errors of proposed method and the FDM can reach the same order of magnitude when the number of nodes is small, but with the increase of the number of nodes, the errors of the FPM decreases faster than those of the FDM. Table 7 gives the ε2 errors of the FPM and the FDM at different times. It can be seen from Table 7 that when T = 2, the results of the FPM are still better than those of the FDM, and the ε2 errors can still reach the order of magnitude of e − 004. Table 8 verifies the convergence order of the FPM. From Table 8, it can be seen that the convergence order is closer to 1.2 as the time step decreases when α = 0.8.Figure 2
Comparing the proposed method with the exact solution whenT = 1, ∆t = 10−3, α = 0.8, and the number of nodes is 21 × 21. (a) Exact solution. (b) FPM numerical solution.
(a)
(b)Figure 3
The equipotentia of the exact solution and the numerical solution of FPM whenT = 1, ∆t = 10−3, α = 0.8, and the number of nodes is 31 × 31. (a) FPM numerical solution. (b) Exact numerical.
(a)
(b)Table 5
Comparison ofε2 errors between the FPM and the FDM under different α, when T = 1, ∆t = 10−3, and the number of nodes is 21 × 21.
α
0.4
0.5
0.6
0.7
0.8
FPM
1.4446e − 04
1.4556e − 04
1.4688e − 04
1.4846e − 04
1.5078e − 04
FDM
7.7600e + 03
3.9952e − 04
4.0953e − 04
4.2217e − 04
4.4397e − 04Table 6
Comparison ofε2 errors between FPM and FDM under different node numbers, when α = 0.8.
Time step
Time
Number of nodes
FPM
FDM
0.001
1
11 × 11
8.5446e − 04
4.4956e − 04
21 × 21
1.5078e − 04
4.4397e − 04
31 × 31
5.1900e − 05
4.4465e − 04
41 × 41
2.3658e − 05
4.4233e − 04Table 7
Comparison ofε2 errors between FPM and FDM under different times, when α = 0.8.
Time step
Number of nodes
Time
FPM
FDM
0.001
21 × 21
0.01
1.5661e − 05
6.9692e − 05
0.1
7.5940e − 05
4.0981e − 04
1.0
1.5078e − 04
4.4397e − 04
2.0
1.6702e − 04
2.9891e − 04Table 8
ε1 errors and convergence order of numerical solution, when α = 0.8 and the number of nodes is 21 × 21.
Time step
ε1
Covergence rate
1/10
1.8574e − 04
1/20
1.0418e − 04
0.8342
1/40
5.7024e − 05
0.8694
1/80
2.7958e − 05
1.0283
1/160
1.3560e − 05
1.0440


## 6. Conclusions

In this article, the finite point method (FPM) was applied to the time fractional Richards’ equation in one- and two-dimensional cases. We have discretized the time fractional derivative using a finite difference formula and obtained a time-discrete scheme which proved to be conditionally stable. One- and two-dimensional numerical examples were given. Compared with finite difference method, the accuracy of finite point method was better than that of finite difference method. Meanwhile, numerical examples showed that the proposed method is simple and effective.



---

*Source: 8358176-2019-09-29.xml*