Shuyun Yang¹, Yaxing Xiao¹, Yifan Cao¹, Yonghui Qin^1,2*
1Guangxi Colleges and Universities Key Laboratory of Data Analysis and Computation, Guangxi Key Laboratory of Automatic Detecting Technology and Instruments, College of Mathematics and Computing Science, Guilin University of Electronic Technology, Guilin, China.
²Center for Applied Mathematics of Guangxi, Guilin University of Electronic Technology, Guilin, China.
DOI: 10.4236/jamp.2026.142033   PDF    HTML   XML   104 Downloads   399 Views   Citations

Abstract

Darcy’s law is the fundamental equation describing the flow of a fluid through a porous medium. Combined with the principle of mass conservation, it leads to the diffusion Equation (e.g., the groundwater flow equation). In this paper, the Legendre-Petrov-Galerkin method is developed for solving the parabolic problem with Dirichlet boundary conditions based on Darcy’s law-preserving. This problem is transformed into an equivalent first-order system by introducing a flux based on Darcy’s law. Our scheme is based on the Legendre Galerkin method, and the right hand side term is processed using the Legendre/Chebyshev-Gauss-Lobatto points. The time direction is approximated by the Crank-Nicolson method. The algebraic system with a sparse coefficient matrix is obtained by selecting the appropriate basis function. Error estimate of the semi-discrete scheme is given by using Gronwall’s inequality (integral form) and Darcy’s law. Numerical examples show that our scheme has the high-order spectral accuracy and it preserves Darcy’s law.

Keywords

Parabolic Problem, Legendre-Petrov-Galerkin, Crank-Nicolson, Darcy’s Law, Legendre/Chebyshev Gauss-Lobatto Points

Share and Cite:

1. Introduction

Darcy’s law serves as the fundamental mathematical model for the flow of an incompressible fluid through a porous medium. This equation is used to describe subsurface water movement. It was formulated by Henry Darcy in 1856 based on experiments with water flowing through sand beds [1]. In this famous literature, let $P$ be the Volumetric flow rate and $K$ denote the Hydraulic conductivity, respectively. $A$ is denoted by cross-sectional area perpendicular to the flow. Let $\text{Δ}h$ be difference in hydraulic head between two points (head represents fluid potential energy) and its distance be denoted by $L$ . Thus, the most common form of the equation is:

$P=-KA\frac{\text{Δ}h}{L},$ (1.1)

where the term $\text{Δ}h/L$ is called the hydraulic gradient. Thus, as in [2] [3], (1.1) is simplified the equation to

$P=-{\kappa }^{\frac{1}{2}}\nabla U$ .(1.2)

The model combines Darcy’s law, which relates the fluid velocity (or flux) to the pressure gradient with the mass conservation equation [1] [4] [5]. It establishes the scientific foundation for the concept of fluid permeability in geosciences and hydrogeology. Although originally established experimentally as an expression of momentum conservation, Darcy’s law has since been derived from the Navier-Stokes equations through homogenization.

Nowadays, the numerical solution of Darcy’s model has been extensively studied in [6]-[9] and the references therein. Crucially, the appropriate functional setting for porous media flow takes the velocity field in $H\left(\text{div},\Omega \right)$ and the pressure in $L^2\left(\Omega \right)$ . This leads to a saddle-point problem, which is well-posed due to the inf-sup conditions satisfied at the continuous level, and allows one to derive stability estimates for both the pressure and the divergence of the velocity as [10]. Therefore, the research on high-order and stable methods for the problem satisfying Darcy’s flow still poses certain challenges, and the current research results are relatively limited. In this paper, we proposed the spectral-Petrov-Galerkin method for the parabolic equation based on Darcy’s law.

In this paper, we consider the spectral-Petrov-Galerkin method for the following model satisfying Darcy’s law as

$\{\begin{array}{ll}{U}_t-\nabla \left(\kappa \nabla U\right)=f\left(x,y,t\right),\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right],\hfill \\ U\left(x,y,t\right)=0,\hfill & \left(x,y,t\right)\in \partial \Omega \times \left(0,T\right],\hfill \\ n\cdot \nabla U=0,\hfill & \left(x,y,t\right)\in \partial \Omega \times \left(0,T\right],\hfill \\ U\left(x,y,0\right)=U_0\left(x,y\right),\hfill & \left(x,y\right)\in \Omega .\hfill \end{array}$ (1.3)

where $\Omega =\left(a,b\right)\times \left(c,d\right)$ , $\partial \Omega$ is a Lipschitz continuous boundary, and $n$ is the unit outward normal vector to $\partial \Omega$ . Herein, the scalar-valued coefficient $\kappa$ is denoted by the symmetric positive definite matrix-valued permeability field and $f\left(x,y,t\right)$ is the given smooth function. Now, we focus on satisfying Darcy’s law (1.2). Thus, an equivalent first order system of (1.3) can be given as

$\{\begin{array}{lll}{U}_t+{\kappa }^{\frac{1}{2}}\nabla \cdot P=f\left(x,y,t\right),\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right],\hfill & \left(\text{mass conservation}\right),\hfill \\ P+{\kappa }^{\frac{1}{2}}\nabla U=0,\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right],\hfill & \left(\text{Darcy}’\text{s law}\right)\hfill \end{array}$ (1.4)

with the same initial value

$U\left(x,y,0\right)=U_0\left(x,y\right),P\left(x,y,0\right)=-{\kappa }^{\frac{1}{2}}\nabla U_0\left(x,y\right),$

and the boundary conditions as (1.3).

Recently, spectral and spectral element methods have been widely applied to solve numerical solutions of partial differential equations arising in many scientific research and engineering technology, such as, the electromagnetic scattering [11]-[14], the eigenvalue problems [15]-[18], the multiple solutions problem [19] [20], and the Volterra integral equations [21] [22]. More recently, the integration of spectral concepts with modern machine learning techniques has emerged, as seen in the application of Chebyshev polynomials as activation functions in deep neural networks for solving neuro-cognitive models [23], and the use of the second Chebyshev wavelet method for inverse nodal problems [24].

The Petrov-Galerkin method is an efficient algorithm, which has been widely used in many problems [25]-[27]. The test function processing method of Petrov-Galerkin spectral method is similar to that of tau method, that is, the test function of this method does not have to satisfy the boundary condition. In [28], this method is applied to a class of nonlocal convection-dominated diffusion problems. Discontinuous Petrov-Galerkin method with perfectly matched layers for time-harmonic problems posed on unbounded domains is given in [29]. In [30], the Legendre Petrov-Galerkin and collocation method for the generalized Korteweg-de Vries equation is derived; the results show that the method has good stability.

In this paper, the Legendre-Petrov-Galerkin spectral method is proposed for solving the parabolic problem (1.3) based on the first order system. The fully-discrete scheme is given by using the Crank-Nicolson method in time discretization. By using Darcy’s law that is used for the flow of a viscous fluid in a permeable medium as in [2], the proposed approach solves the solution and its flux simultaneously, which is analogous to the famous mixed finite element method or discontinuous Galerkin method [31]-[34]. It is well-known that the mixed finite element methods or the discontinuous Galerkin method have been widely applied to the numerical solutions of the partial differential equation arising in the field of physical and engineering problems [31]-[33] [35]. As point out in [36] [37], the numerical scheme based on this framework is proved to obey mass conservation and unconditional stability. It is worth emphasizing that preserving Darcy’s law at the discrete level is crucial for physical simulations, as it ensures a physically consistent velocity field and guarantees local mass conservation, which are fundamental for obtaining reliable results in applications such as groundwater flow and reservoir modeling. The stability analysis of the semi-discrete scheme is derived based on Gronwall’s inequality (integral form) and Darcy’s law, the convergence analysis is also obtained. Our numerical examples are given to test the spectral accuracy of our scheme.

This paper is organized as follows. In Section 2, we introduce some notations and present the Legendre-Petrov-Galerkin spectral scheme for the parabolic equation based on the first-order reformulation. Section 3 is devoted to the implementation details of the proposed method. The stability and convergence analysis of the semi-discrete scheme are provided in Section 4. Section 5 describes the fully discrete scheme and its implementation aspects. Section 6 presents numerical examples to illustrate the spectral accuracy and conservation properties of our method. Finally, concluding remarks are given in Section 7.

2. Notations and Numerical Schemes

2.1. Notations

In this paper, standard notations are used for Sobolev spaces and corresponding norms. Thus, let $H^{\sigma }\left(\Lambda \right)$ be the classical weighted Sobolev space with a integer $\sigma \ge 0$ , where its norm and seminorm are denoted by ${\Vert \cdot \Vert }_{L^{\sigma }\left(\Lambda \right)}$ and ${\left|\cdot \right|}_{L^{\sigma }\left(\Lambda \right)}$ , respectively. We set $L^2\left(\Lambda \right)=H^0\left(\Lambda \right)$ if $\sigma =0$ , where the inner product and the norm are denoted by ${\left(\cdot ,\cdot \right)}_{\Lambda }$ and ${\Vert \cdot \Vert }_{\Lambda }$ , respectively. Further, we set $\Omega :={\Lambda }_x\times {\Lambda }_y$ , where ${\Lambda }_x:=\left(a,b\right)$ and ${\Lambda }_y:=\left(c,d\right)$ . Now, let us introduce the Sobolev space that are tensor products of Sobolev space on ${\Lambda }_{\xi }$ , where $\xi =x$ or $y$ . For simplicity, let ${\Omega }_d$ stand for the open hypercube $\Omega$ , where $d=1,2,3$ . The generic point in ${\Omega }_d$ is denoted by $x=\left(x_1,x_2,\cdots ,x_d\right)$ . For any $v\in L^2\left(\Omega \right)$ , we associate the $d$ -functions $v_i$ defined $v_i\left(x_i\right)\left(x_1,\cdots ,x_{i-1},x_{i+1},\cdots ,x_d\right)=v\left(x_1,x_2,\cdots ,x_d\right)$ , $1\le i\le d$ . For any nonnegative real numbers $r$ and $s$ , the Sobolev space is defined as

$H^r\left({\Lambda }_i;H^s\left({\Omega }_{d-1}\right)\right):=\left\{v\in L^2\left(\Omega \right):v_i\in H^r\left({\Lambda }_i;H^s\left({\Omega }_{d-1}\right)\right)\right\}.$

Thus, we derive

$L^2\left({\Omega }_d\right)=L^2\left({\Lambda }_i;L^2\left({\Omega }_{d-1}\right)\right),H^r:=H^j\left({\Lambda }_i;H^{r-j}\left({\Omega }_{d-1}\right)\right).$

Next, we fixed $d=2$ for the convenience of narration. Thus, for any $v\in L^2\left(\Omega \right)$ , let

$H_0^1\left(\Omega \right)=\left\{v\in H^1\left(\Omega \right):{v|}_{\partial \Omega }=0\right\},$

$H^{1,0}\left(\Omega \right)=\left\{v\in H^1\left({\Lambda }_x;L^2\left({\Lambda }_y\right)\right):v\left(x,c,t\right)=v\left(x,d,t\right)=0\right\},$

$H^{0,1}\left(\Omega \right)=\left\{v\in L^2\left({\Lambda }_x;H^1\left({\Lambda }_y\right)\right):v\left(a,y,t\right)=v\left(b,y,t\right)=0\right\},$

and

$H\left(\text{div};\Omega \right):=\left\{v\in L^2{\left(\Omega \right)}^2:\nabla \cdot v\in L^2\left(\Omega \right)\right\},$

$H_0\left(\text{div};\Omega \right):=\left\{v\in H\left(\text{div};\Omega \right):v^{\left(1\right)}\in H^{1,0}\left(\Omega \right),v^{\left(2\right)}\in H^{0,1}\left(\Omega \right)\right\},$

where ${\Vert v\Vert }_{\text{div}}:={\Vert v\Vert }_{H\left(\text{div};\Omega \right)}={\left({\Vert v\Vert }^2+{\Vert \nabla \cdot v\Vert }^2\right)}^{\frac{1}{2}}$ .

2.2. Numerical Scheme

In this subsection, we present our scheme for (1.4). To this end, the weak formula of (1.4) is to find $\left(U,P\right)\in H_0^1\left(\Omega \right)\times H_0\left(\text{div};\Omega \right)$ such that

$\{\begin{array}{ll}\left(U_t,\phi \right)+{\kappa }^{\frac{1}{2}}\left(\nabla \cdot P,\phi \right)=\left(f,\phi \right),\hfill & \phi \in H^1\left(\Omega \right),t\in \left(0,T\right],\hfill \\ \left(P,\psi \right)-{\kappa }^{\frac{1}{2}}\left(U,\nabla \cdot \psi \right)=0,\hfill & \psi \in H_0\left(\text{div};\Omega \right),t\in \left(0,T\right].\hfill \end{array}$ (2.1)

Now, we turn to consider the discrete schemes of the weak form (2.1). Let ${\mathbb{P}}_N\left({\Lambda }_{\xi }\right)$ denote the set of all the polynomials of the degree at most $N$ on $I_{\xi }$ with the variable $\xi :=x$ or $y$ . Thus, we set

${\mathbb{Q}}_N\left(\Omega \right):={\mathbb{P}}_N\left({\Lambda }_x\right)\times {\mathbb{P}}_N\left({\Lambda }_y\right),$

$V_N:=H_0^1\left(\Omega \right)\cap {\mathbb{Q}}_N\left(\Omega \right),$

$W_N:=H_0\left(\text{div},\Omega \right)\cap {\left[{\mathbb{Q}}_N\left(\Omega \right)\right]}^2.$

For simplicity, we also introduce the LGL or CGL nodes on $\overline{\Omega }$ , which are well-suited for spectral approximations due to their high-order accuracy. They are defined by:

$x_j:=\left(x_{j_1},y_{j_2}\right),j:=\left(j_1,j_2\right),0\le j_1,j_2\le N.$

As in [38] [39], the LGL interpolation operator or CGL interpolation operator on $\Omega$ is defined by $ℐℐ:C\left(\overline{\Omega }\right){\to }_N\left(\Omega \right)$ such that

$\left[ℐ{ℐ}_{N}v\right]\left(x_j\right)=v\left(x_j\right),j\in {\mathbb{Z}}^2\cap {\left[0,N\right]}^2,$

where $ℐ{ℐ}_{N}v:={ℐ}_{N,x}\circ {ℐ}_{N,y}$ . Therefore, the semi-discrete scheme of (2.1) is to find $\left(u_N\left(t\right),p_N\left(t\right)\right)\in V_N\times W_N$ such that

$\{\begin{array}{ll}\left(u_{Nt},\phi \right)+{\kappa }^{\frac{1}{2}}\left(\nabla \cdot p_N,\phi \right)=\left(ℐ{ℐ}_{N}f,\phi \right),\hfill & \phi \in V_N,t\in \left(0,T\right],\hfill \\ \left(p_N,\psi \right)-{\kappa }^{\frac{1}{2}}\left(u_N,\nabla \cdot \psi \right)=0,\hfill & \psi \in W_N,t\in \left(0,T\right],\hfill \end{array}$ (2.2)

where the initial value $u_N\left(x,y,0\right)=ℐ{ℐ}_N{U}_0$ and $p_N\left(x,y,0\right)=-ℐ{ℐ}_N\left[{\kappa }^{1/2}\nabla U_0\right]$ .

To obtain the fully-discrete scheme, we show our method by using the Crank-Nicolson discretization in time for (2.2). Now, let us introduce a general nonuniform in time variable $t$ as

$0=t_0<t_1<t_2<\cdots <t_{n_T}=T,n_T=\lfloor T/\tau \rfloor ,$

where $\tau :={\tau }_k=t_k-t_{k-1},0\le k\le n_T$ . Thus, we set

$S_{\tau }:=\left\{k\tau ,k=0,1,\cdots ,n_T,n_T\tau =T\right\}.$

For simplicity, we denote $v^k\left(x\right):=v\left(x,y,t_k\right)$ by $v^k$ and define

$v_t^k=\frac{1}{\tau }\left(v^{k+1}-v^k\right),{\overline{v}}^k=\frac{1}{2}\left(v^{k+1}+v^k\right).$(2.3)

Thus, the fully-discrete scheme of (2.1) is to find $\left(u_N^k,p_N^k\right)\in V_N\times W_N$ such that

$\{\begin{array}{ll}\left(u_{Nt}^k,\phi \right)+{\kappa }^{\frac{1}{2}}\left(\nabla \cdot {\overline{p}}_N^k,\phi \right)=\left(ℐ{ℐ}_N{\overline{f}}^k,\phi \right),\hfill & \phi \in V_N,t\in S_{\tau }\hfill \\ \left({\overline{p}}_N^k,\tau \psi \right)-{\kappa }^{\frac{1}{2}}\left({\overline{u}}_N^k,\tau \nabla \cdot \psi \right)=0,\hfill & \psi \in W_N,t\in S_{\tau },\hfill \end{array}$(2.4)

where the initial value $u_N^0\left(x,y,0\right)=ℐ{ℐ}_N{U}_0$ and $p_N^0\left(x,y,0\right)=-ℐ{ℐ}_N\left[{\kappa }^{\frac{1}{2}}\nabla U_0\right]$ .

3. Implementation

In this subsection, a simple description of the implementation of our scheme (2.4) is presented. Let ${\left\{L_i\left(\xi \right)\right\}}_{i=0}^N$ be the set of all Legendre polynomials with degree less than $N$ , which satisfy the orthogonality relation as

${\displaystyle {\int }_{-1}^1{L}_i\left(\xi \right)L_j\left(\xi \right)\text{d}x}=\{\begin{array}{ll}\frac{2}{2i+1},\hfill & i=j,\hfill \\ 0,\hfill & i\ne j.\hfill \end{array}$

As in [40], we recommend the following basis functions

${\varphi }_0\left(\xi \right)=\frac{1}{2}\left(1-\xi \right),{\varphi }_1\left(\xi \right)=\frac{1}{2}\left(1+\xi \right)$ (3.1)

for the boundary conditions and

${\varphi }_i\left(\xi \right)=L_i\left(\xi \right)-L_{i-2}\left(\xi \right),i\ge 2$ (3.2)

for the interior basis functions of ${\Lambda }_{\xi }$ , where $\xi =x,y$ . Thus, let us introduce the mass and stiffness matrices as

$M={\left(M_{ij}\right)}_{i,j=0}^N,S={\left(S_{ij}\right)}_{i,j=0}^N,K={\left(K_{ij}\right)}_{i,j=0}^N,$

where

$M_{ij}=\left({\varphi }_i\left(\xi \right),{\varphi }_j\left(\xi \right)\right),S_{ij}=\left({{\varphi }^{\prime }}_i\left(\xi \right),{{\varphi }^{\prime }}_j\left(\xi \right)\right),K_{ij}=\left({\varphi }_i\left(\xi \right),{{\varphi }^{\prime }}_j\left(\xi \right)\right).$

Now, we present the algebra system of our schemes (2.4), respectively. For any $\left(u_N,p_N\right)\in V_N\times W_N$ in (2.4), by using the tensor-type basis function of (3.1) and (3.2), $\left(u_N,p_N\right)$ are expanded into

$u_N={\displaystyle \sum }_{k,l=2}^N{\widehat{u}}_{kl}\left(t\right){\varphi }_k\left(x\right){\varphi }_l\left(y\right),$(3.3)

$p_N=\left(\begin{array}{c}{p}_N^{\left(1\right)}\\ p_N^{\left(2\right)}\end{array}\right)={\displaystyle \sum }_{k=0}^N{\displaystyle \sum }_{l=2}^N{\widehat{p}}_{k,l}^{\left(1\right)}\left(t\right){\Phi }_{kl}^{\left(1\right)}+{\displaystyle \sum }_{k=2}^N{\displaystyle \sum }_{l=0}^N{\widehat{p}}_{k,l}^{\left(2\right)}\left(t\right){\Phi }_{kl}^{\left(2\right)},$(3.4)

where ${\widehat{U}}^k,{\widehat{P}}_1^k$ and ${\widehat{P}}_2^k$ are the expansion coefficient of the approximate solution in the basis function,

${\Phi }_{kl}^{\left(1\right)}:=\left(\begin{array}{c}{\varphi }_k\left(x\right){\varphi }_l\left(y\right)\\ 0\end{array}\right),{\Phi }_{kl}^{\left(2\right)}:=\left(\begin{array}{c}0\\ {\varphi }_k\left(x\right){\varphi }_l\left(y\right)\end{array}\right).$

By applying the Crank-Nicolson method to transform (2.3) into (2.4), we obtain the following algebraic system as

$A{x}^{k+1}=B{x}^k+f,$ (3.5)

where

$x^k=\left(\begin{array}{c}{\widehat{u}}^k\\ {\widehat{p}}_1^k\\ {\widehat{p}}_2^k\end{array}\right),\widehat{u}:=vec\left(\widehat{U}\right),{\widehat{p}}_1:=vec\left({\widehat{P}}_1\right),{\widehat{p}}_2:=vec\left({\widehat{P}}_2\right),$

$\begin{array}{l}A=\left(\begin{array}{ccc}{M}^{\text{T}}\otimes M& \overline{\kappa }{K}^{\text{T}}\otimes M& \overline{\kappa }{M}^{\text{T}}\otimes K^{\text{T}}\\ \overline{\kappa }{K}^{\text{T}}\otimes M& \overline{\tau }{M}^{\text{T}}\otimes M& 0\\ \overline{\kappa }{M}^{\text{T}}\otimes K^{\text{T}}& 0& \overline{\tau }{M}^{\text{T}}\otimes M\end{array}\right),\\ B=\left(\begin{array}{ccc}{M}^{\text{T}}\otimes M& -\overline{\kappa }{K}^{\text{T}}\otimes M& -\overline{\kappa }{M}^{\text{T}}\otimes K^{\text{T}}\\ -\overline{\kappa }{K}^{\text{T}}\otimes M& -\overline{\tau }{M}^{\text{T}}\otimes M& 0\\ -\overline{\kappa }{M}^{\text{T}}\otimes K^{\text{T}}& 0& -\overline{\tau }{M}^{\text{T}}\otimes M\end{array}\right).\end{array}$

where $\overline{\kappa }={\kappa }^{\frac{1}{2}}$ , $\otimes$ denotes the tensor product of matrices. Finally, some suitable solvers can be used to (3.5) by depending boundary conditions of (1.3).

4. Preliminaries

Herein, let ${\pi }_N^{\xi }:L^2\left(I\right)\to {\mathbb{P}}_N\left({\Lambda }_{\xi }\right)$ be the orthogonal projection operator on ${\Lambda }_{\xi }$ as in [38] [39] [41] [42], where $\xi =x$ or $y$ , namely, for any $w\in L^2\left({\Lambda }_{\xi }\right)$ , we have

${\left({\pi }_N^{\xi }w-w,v\right)}_{L^2\left({\Lambda }_{\xi }\right)}=0,v{\in }_N\left({\Lambda }_{\xi }\right).$ (4.1)

Further, we define ${\Pi }_N:={\pi }_N^x\circ {\pi }_N^y$ and $\Omega :={\Lambda }_x\otimes {\Lambda }_y$ . Thus, as in ([38], (5.8.13) and (5.8.14)), ${\Pi }_N:L^2\left(\Omega \right){\to }_N\left(\Omega \right)$ is the $L^2$ -orthogonal projection operator on $\Omega$ , for any $w\in L^2\left(\Omega \right)$ ,

$\left({\Pi }_{N}w-w,v\right)=0,\forall v\in {\mathbb{Q}}_N^0\left(\Omega \right).$ (4.2)

For any $w\in H^r\left(\Omega \right)\cap H_0^1\left(\Omega \right)$ and $l=0,1$ , the following results are given in ([41], Theorem 7.2 and Theorem 14.2 and [39]),

${\Vert {\Pi }_{N}w-w\Vert }_{H^l\left(\Omega \right)}\le C{N}^{l-r}{\left|w\right|}_{H^r\left(\Omega \right)},r\ge 1,$ (4.3)

${\Vert ℐ{ℐ}_{N}w-w\Vert }_{H^l\left(\Omega \right)}\le C{N}^{l-r}{\left|w\right|}_{H^r\left(\Omega \right)},r\ge d.$ (4.4)

5. Stability and Convergence for Our Semi-Discrete Scheme

In this section, we give a simple proof of our semi-discrete Scheme (2.2). Here, let ${\tilde{u}}_N$ and ${\tilde{p}}_N$ be the errors of $u_N$ and $p_N$ , respectively. Assume that $\tilde{f}$ is the error of $f$ . Thus, from (2.1) and (2.2) we have the error equation as $\left({\tilde{u}}_N\left(t\right),{\tilde{p}}_N\left(t\right)\right)\in V_N\times W_N$ such that, for any $t\in \left(0,T\right]$ ,

$\{\begin{array}{ll}\left({\tilde{u}}_N,\phi \right)+{\kappa }^{\frac{1}{2}}\left(\nabla \cdot {\tilde{p}}_N,\phi \right)=\left(\tilde{f},\phi \right),\hfill & \phi \in V_N,\hfill \\ \left({\tilde{p}}_N,\psi \right)-{\kappa }^{\frac{1}{2}}\left({\tilde{u}}_N,\nabla \cdot \psi \right)=0,\hfill & \psi \in W_N.\hfill \end{array}$ (5.1)

If taking $\phi ={\tilde{u}}_N$ and $\psi ={\tilde{p}}_N$ in (5.1), we have

$\{\begin{array}{l}{\Vert {\tilde{p}}_N\Vert }^2={\kappa }^{\frac{1}{2}}\left({\tilde{u}}_N,\nabla \cdot {\tilde{p}}_N\right),\\ \left({\tilde{u}}_N,{\tilde{u}}_N\right)+{\kappa }^{\frac{1}{2}}\left(\nabla \cdot {\tilde{p}}_N,{\tilde{u}}_N\right)=\left(ℐ{ℐ}_{N}f,{\tilde{u}}_N\right).\end{array}$ (5.2)

Therefore, we derive

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert {\tilde{u}}_N\Vert }^2+{\Vert {\tilde{p}}_N\Vert }^2=\left(\tilde{f},{\tilde{u}}_N\right).$

By the Cauchy–Schwarz and Young inequalities

$\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert {\tilde{u}}_N\Vert }^2+{\Vert {\tilde{p}}_N\Vert }^2\le \frac{1}{2}{\Vert \tilde{f}\Vert }^2+\frac{1}{2}{\Vert {\tilde{u}}_N\Vert }^2.$(5.3)

Integrating in time of (5.4), we obtain

${\Vert {\tilde{u}}_N\left(t\right)\Vert }^2+2{\displaystyle {\int }_0^t{\Vert {\tilde{p}}_N\left(s\right)\Vert }^2\text{d}s}\le {\Vert {\tilde{u}}_N\left(0\right)\Vert }^2+{\displaystyle {\int }_0^t{\Vert \tilde{f}\left(s\right)\Vert }^2\text{d}s}+{\displaystyle {\int }_0^t{\Vert {\tilde{u}}_N\left(s\right)\Vert }^2\text{d}s}.$(5.4)

Let $E\left(t\right):={\Vert {\tilde{u}}_N\left(t\right)\Vert }^2$ and $R\left(t\right):={\Vert {\tilde{u}}_N\left(0\right)\Vert }^2+{\displaystyle {\int }_0^t{\Vert \tilde{f}\left(s\right)\Vert }^2\text{d}s}$. Thus, we have

$E\left(t\right)\le R\left(t\right)+{\displaystyle {\int }_0^{t}E\left(s\right) \text{d}s}.$

By Gronwall’s inequality (integral form), we obtain

$E\left(t\right)\le R\left(t\right)+{\displaystyle {\int }_0^{t}R\left(s\right){\text{e}}^{t-s}\text{d}s}.$

Noting that $R\left(t\right)$ is non-decreasing for any $t\in \left(0,T\right]$ and $R\left(s\right)\le R\left(t\right)$ for $s\le t$ , we have

$E\left(t\right)\le R\left(t\right)+R\left(t\right){\displaystyle {\int }_0^t{\text{e}}^{t-s}\text{d}s}=R\left(t\right){\text{e}}^t.$

Therefore, from above inequality and (5.6), we derive that for any $t\in \left(0,T\right]$ ,

$\begin{array}{c}{\Vert {\tilde{u}}_N\left(t\right)\Vert }^2+{\displaystyle {\int }_0^t{\Vert {\tilde{p}}_N\left(s\right)\Vert }^2\text{d}s}\le {\Vert {\tilde{u}}_N\left(t\right)\Vert }^2+2{\displaystyle {\int }_0^t{\Vert {\tilde{p}}_N\left(s\right)\Vert }^2\text{d}s}\\ \le {\text{e}}^t\left({\Vert u_N\left(0\right)\Vert }^2+{\displaystyle {\int }_0^t{\Vert \tilde{f}\left(s\right)\Vert }^2\text{d}s}\right).\end{array}$(5.5)

For a fixed time interval $\left[0,T\right]$ , we take $C=e^T$ to get

${\Vert {\tilde{u}}_N\left(t\right)\Vert }^2+{\displaystyle {\int }_0^T{\Vert {\tilde{p}}_N\left(t\right)\Vert }^2\text{d}t}\le C\left({\Vert {\tilde{u}}_N\left(0\right)\Vert }^2+{\displaystyle {\int }_0^T{\Vert \tilde{f}\left(t\right)\Vert }^2\text{d}t}\right).$(5.6)

Next, let $U^{*}={\Pi }_{N}U$ and $P^{\text{*}}={\Pi }_{\text{div},N}P$ . Thus, we set

$e_u:=u_N-U^{\text{*}},{\eta }_u:=U^{\text{*}}-U,e_p:=p_N-P^{\text{*}},{\eta }_p:=P^{\text{*}}-P.$

By utilizing the properties of orthogonal projection (4.1), we have

$\{\begin{array}{ll}\left(U^{\text{*}},\phi \right)+{\kappa }^{\frac{1}{2}}\left(\nabla \cdot P^{\text{*}},\phi \right)=\left(f,\phi \right)+{\kappa }^{\frac{1}{2}}\left(\nabla \cdot {\eta }_p,\phi \right),\hfill & \phi \in V_N,\hfill \\ \left(P^{\text{*}},\psi \right)-{\kappa }^{\frac{1}{2}}\left(U^{\text{*}},\nabla \cdot \psi \right)=0,\hfill & \psi \in W_N.\hfill \end{array}$ (5.7)

From (2.2) and (5.7), for any $t\in \left(0,T\right]$ the error equation is obtained as

$\{\begin{array}{ll}\left(e_u,\phi \right)+{\kappa }^{\frac{1}{2}}\left(\nabla \cdot e_p,\phi \right)=\left(ℐ{ℐ}_{N}f-f,\phi \right)+{\kappa }^{\frac{1}{2}}\left({\eta }_p,\nabla \phi \right),\hfill & \phi \in V_N,\hfill \\ \left(e_p,\psi \right)-{\kappa }^{\frac{1}{2}}\left(e_u,\nabla \cdot \psi \right)=0,\hfill & \psi \in W_N.\hfill \end{array}$ (5.8)

Taking $\phi =e_u$ and $\psi =e_p$ for $\nabla e_u=e_p$ in (5.8), we have

$\begin{array}{l}\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert e_u\Vert }^2+{\Vert e_p\Vert }^2\le \frac{1}{2}{\Vert {ℐ}_{N}f-f\Vert }^2+\frac{1}{2}{\Vert e_u\Vert }^2+\frac{\kappa }{2}{\Vert {\eta }_p\Vert }^2+\frac{1}{2}{\Vert \nabla e_u\Vert }^2\\ \stackrel{\left(4.3\right),\left(4.4\right)}{\le }C{N}^{-2s}{\Vert f\Vert }_{H^s\left(\Omega \right)}^2+\frac{1}{2}{\Vert e_u\Vert }^2+C{N}^{-2r}{\Vert P\Vert }_{H^r\left(\Omega \right)}^2+\frac{1}{2}{\Vert e_p\Vert }^2.\end{array}$

As in (5.6), we have

$\begin{array}{l}{\Vert e_u\left(t\right)\Vert }^2+{\displaystyle {\int }_0^T{\Vert e_p\left(t\right)\Vert }^2\text{d}t}\\ \le C{\Vert e_u\left(0\right)\Vert }^2+{\displaystyle {\int }_0^T\left(C{N}^{-2s}{\Vert f\Vert }_{H^s\left(\Omega \right)}^2+C{N}^{-2r}{\Vert P\Vert }_{H^r\left(\Omega \right)}^2\text{d}t\right)}.\end{array}$

Now, we consider the error estimate $e\left(t\right)$ with $t=0$ by using (4.3) and (4.4),

$\begin{array}{c}\Vert e\left(0\right)\Vert =\Vert \left(ℐ{ℐ}_N-{\Pi }_N\right)U_0\Vert \\ \le \Vert U_0-ℐ{ℐ}_N{U}_0\Vert +\Vert U_0-{\Pi }_N{U}_0\Vert \\ \le C{N}^{-r}{\Vert U_0\Vert }_{H^r\left(\Omega \right)}.\end{array}$

Theorem 5.1. Assume that $r,s>1$ and

$U\in H^1\left(0,T;H_0^1\left(I\right)\cap H^r\left(I\right)\right)\cap C\left(0,T;C^1\left[0,T\right]\right),$

$P\in L^2\left(0,T;H_0\left(\text{div};\Omega \right)\cap H^r\left(\Omega \right)\right),$

$g\in L^2\left(0,T;H^s\left(I\right)\right).$

Then there exists a positive constant $C$ such that

$\Vert u\left(t\right)-U\left(t\right)\Vert \le C{N}^{-r}{\Vert P\Vert }_{H^r\left(\Omega \right)}+C{N}^{-s}{\Vert f\Vert }_{H^s\left(\Omega \right)},t\in \left(0,T\right].$

6. Numerical Examples

In this section, some numerical examples are given to verify the effectiveness and high accuracy of our method. Darcy’s law (1.2) with approximate solutions $\left(p_N,u_N\right)$ is also verified via numerical examples. For simplicity, we define $E_u:=u-u_N$ , $E_p:=p-p_N$ . Herein, we apply the Legendre-Petrov-Galerkin method (2.4) with the interpolation on CGL nodes to solve problems (1.3), and compare the results with other spectral methods.

Example 6.1. Consider (1.3) with $\kappa =1$ , $\left(x,t\right)\in \left(0,\pi \right)\times \left(0,1\right]$ , and the initial value condition $U\left(x,0\right)=\sin \left(12x\right)-0.5\sin \left(8x\right)$ . The Dirichlet boundary condition is given as: $u\left(0,t\right)=u\left(\pi ,t\right)=0$ . Assume that the exact solution and Darcy’s equation are given by

$U\left(x,t\right)=\sin \left(12x\right){\text{e}}^t-0.5\sin \left(8x\right){\text{e}}^{-2t},$

$P\left(x,t\right)=-12\cos \left(12x\right){\text{e}}^t+4\cos \left(8x\right){\text{e}}^{-2t}.$

In this example, we confirm the spectral accuracy of the proposed method (2.4) and Darcy’s law with respect to the corresponding numerical solution.

Figure 1 shows images of the exact solution $U$ and the approximate solution $u_N$ at the final time $T=0.1$ , and its corresponding error function $u_N-U$ is also given, where $N=64$ and $\tau ={10}^{-5}$ . The maximum errors and $L^2$ -errors of $E_u$ and $E_p$ at $T=1$ are given in Table 1, where the time steps are taken $\tau ={10}^{-3}$ and $\tau ={10}^{-4}$ , the degree of polynomial $N$ is increasing from 28 to 64. Note that our scheme has high order accuracy and is an effective method from the numerical results given in Figure 1 and Table 1. From Figure 2, we can see that all the error ratios in $L^2$ -norm of $p_N+\nabla u_N$ are less than or equal to 10⁻⁹. Thus, the discrete form of Darcy’s law (1.2) is satisfied.

(a)

(b)

Figure 1. The plot of solutions and their corresponding error function at $T=0.1$ in Example 6.1, where $N=64$ and $\tau ={10}^{-5}$ , (a) The image of the exact solution $U\left(T\right)$ and the approximate solution $u_N\left(T\right)$ , For the exact solution $U$ and its numerical solution $u_N$ ; (b) The image of the error function $E_u$ , For the error function $E_u$ .

Figure 2. The $L^2$ -error evolution for Darcy’s flow with our scheme taking time step-size as $\tau \in \left\{0.1,{10}^{-2},{10}^{-3}\right\}$ for Example 6.1, where $N=128$ . (a) for $\tau =0.1$ , (b) for $\tau ={10}^{-2}$ , (c) for $\tau ={10}^{-3}$ .

Table 1. The maximum errors and $L^2$ -errors of Scheme (2.4) with $\tau ={10}^{-3}$ and $\tau ={10}^{-4}$ at $T=1$ in Example 6.1.

$\tau$	$N$	$\Vert E_u\Vert$	$\Vert E_p\Vert$	${\Vert E_u\Vert }_{\infty }$	${\Vert E_p\Vert }_{\infty }$
10−3	28	1.9495e−04	6.9862e−02	1.4481e−04	4.7403e−02
	32	2.7855e−06	1.1517e−03	2.0081e−06	7.8518e−04
	64	2.2190e−09	5.1843e−08	2.1977e−09	2.4567e−08
10−4	28	1.9495e−04	6.9862e−02	1.4481e−04	4.7403e−02
	32	2.7855e−06	1.1517e−03	2.0100e−06	7.8520e−04
	64	2.2230e−11	5.2136e−10	2.2033e−11	2.4843e−10

Example 6.2. Consider the parabolic problem (1.3) with some different constant coefficients $\kappa$ and $\left(x,t\right)\in \left(-1,1\right)\times \left(0,1\right]$ , where the exact solution and Darcy’s equation are given by

$U\left(x,t\right)=\sin \left(\pi x\right)\cos \left(\pi t\right)-\sin \left(2\pi x\right){\text{e}}^{\frac{1}{2}t},$

$P\left(x,t\right)=-{\kappa }^{\frac{1}{2}}\left(\pi \cos \left(\pi x\right)\cos \left(\pi t\right)-2\pi \cos \left(2\pi x\right){\text{e}}^{\frac{1}{2}t}\right).$

In this example, we consider the effectiveness of our scheme (2.4) for the parabolic equation with different values of $\kappa$ .

For $\kappa =0.1$ , $N=28$ and $\tau ={10}^{-5}$ , the images of the exact solution $U$ and its approximate solution $u_N$ at $T=10$ are presented in Figure 3, and the corresponding error function is also shown. In Table 2, we take $N$ increase from 14 to 22 for $\tau ={10}^{-5}$ and take the time step decrease from $\tau ={10}^{-1}$ to $\tau ={10}^{-5}$ for $N=24$ , the $L^2$ -errors for $E_u$ and $E_p$ at $T=1$ for $\kappa =5$ and $\kappa =12$ are listed. Here, we also show that Darcy’s law is satisfied from Figure 4.

Table 2. Errors with the $L^2$ -norm of Scheme (2.4) with $\tau \in \left\{{10}^{-1},{10}^{-3},{10}^{-5}\right\}$ and $T=1$ in Example 6.2 for $\kappa =5$ and $\kappa =12$ .

$\tau$	$N$	$\kappa =5$		$\kappa =12$
		$\Vert E_u\Vert$	$\Vert E_p\Vert$	$\Vert E_u\Vert$	$\Vert E_p\Vert$
10−5	14	1.8479e−05	7.1165e−03	1.8481e−05	1.1025e−02
	16	7.9116e−07	3.5902e−04	7.9124e−07	5.5621e−04
	18	2.5723e−08	1.3639e−05	2.5725e−08	2.1129e−05
	20	6.6188e−10	4.0565e−07	6.6192e−10	6.2844e−07
	22	1.3873e−11	9.7240e−09	1.3871e−11	1.5063e−08
10−1	24	3.2197e−05	5.1087e−04	5.8322e−06	1.4703e−04
10−3		3.1639e−09	5.0238e−08	5.5801e−10	1.4190e−08
10−5		3.2477e−13	1.9095e−10	2.5707e−13	2.9152e−10

(a)

(b)

Figure 3. The plot of solutions and the corresponding error function at $T=10$ in Example 6.2 for $\kappa =0.1$ , where $N=28$ and $\tau ={10}^{-5}$ , (a) The image of the exact solution $U\left(T\right)$ and the approximate solution $u_N\left(T\right)$ , For the exact solution $U$ and its numerical solution $u_N$ ; (b) The image of the error function $E_u$ , For the error function $E_u$ .

Figure 4. $L^2$ -error evolution for Darcy’s flow with taking time step-size as $\tau \in \left\{0.1,{10}^{-2},{10}^{-3}\right\}$ in Example 6.2, where $N=32$ . (a) for $\tau =0.1$ , (b) for $\tau ={10}^{-2}$ , (c) for $\tau ={10}^{-3}$ .

Example 6.3. Consider the parabolic problem (1.3) with $\kappa =1$ and $\left(x,t\right)\in \left(-1,1\right)\times \left(0,1\right].$ Assume that the following exact solution and Darcy’s equation are given as

$U\left(x,t\right)=\exp \left(\frac{1}{t+\alpha }\right){\left(\text{ch}x-\text{ch}1\right)}^2,$

$P\left(x,t\right)=2\exp \left(\frac{1}{t+\alpha }\right)\left(\cosh x-\cosh 1\right)\sinh x.$

Here, we consider the effectiveness of the proposed scheme for solving the parabolic problem (1.3) with the exact solution including a parameter $\alpha$ .

In Figure 5, we draw the plots of the exact solutions $U$ and the approximate solutions $u_N$ , and the error function at $T=1$ , where $N=24$ , $\tau ={10}^{-4}$ , and $\alpha =5$ . In Table 3, for $\tau ={10}^{-4}$ and $\alpha =1.13$ , we show the maximum errors and the $L^2$ -errors of our scheme (2.4) with $N$ increasing from 6 to 14 at $T=1$ .

(a)

(b)

Figure 5. The plot of solutions and its corresponding error function at $T=1$ in Example 6.3, where $N=24$ , $\tau ={10}^{-4}$ and $\alpha =5$ , (a) The image of the exact solution $U\left(T\right)$ and the approximate solution $u_N\left(T\right)$ , For the exact solution $U$ and its numerical solution $u_N$ ; (b) The image of the error function $E_u$ , For the error function $E_u$ .

Table 3. The maximum errors and $L^2$ -errors of Scheme (2.4) with $\tau ={10}^{-4}$ at $T=1$ in Example 6.3 for $\alpha =1.13$ .

$N$	$\tau$	$\Vert E_u\Vert$	$\Vert E_p\Vert$	${\Vert E_u\Vert }_{\infty }$	${\Vert E_p\Vert }_{\infty }$
6	10−4	2.1935e−04	9.3115e−03	3.2534e−04	6.0122e−03
8		2.3022e−06	1.3123e−04	3.0395e−06	8.3223e−05
10		1.7046e−08	1.1689e−06	2.2786e−08	7.2355e−07
12		6.4845e−10	8.7346e−09	8.1825e−10	5.5337e−09
14		6.4608e−10	2.4025e−09	7.7765e−10	1.2105e−09

In Table 4, we give the $L^2$ -errors of Scheme (2.4) for $\alpha =-10$ and $\alpha =10$ , where $\tau$ takes values from 10⁻¹ to 10⁻⁵ for $N=18$ and $N$ is increasing from 8 to 20 for $\tau ={10}^{-4}$ . The numerical results show that our method has high accuracy. Similarly, Figure 6 also shows that our scheme preserves Darcy’s law.

Table 4. The $L^2$ -errors of Scheme (2.4) with $\tau ={10}^{-1}$ , $\tau ={10}^{-3}$ , $\tau ={10}^{-4}$ and $\tau ={10}^{-5}$ at $T=1$ in Example 6.3 for $\alpha =-10,10$ .

$\tau$	$N$	$\alpha =-10$		$\alpha =10$
		$\Vert E_u\Vert$	$\Vert E_p\Vert$	$\Vert E_u\Vert$	$\Vert E_p\Vert$
10−1	18	1.7927e−07	6.5189e−07	1.5496e−07	5.6902e−07
10−3		1.7909e−11	6.5112e−11	1.5496e−11	5.6895e−11
10−5		1.5958e−12	6.8045e−12	3.4015e−13	3.0327e−12
10−4	8	1.2836e−06	1.9702e−05	1.5709e−06	2.3851e−05
	12	5.1595e−11	8.8627e−10	6.3144e−11	1.0721e−09
	16	2.8067e−13	1.1194e−12	1.0870e−13	5.1159e−13
	18	2.7178e−13	1.1923e−12	3.3849e−13	1.4844e−12
	20	2.7511e−13	1.2555e−12	7.0803e−14	2.3958e−13

Figure 6. $L^2$ -error evolution for Darcy’s flow with varying time steps $\tau$ in Example 6.3, where $N=24$ . (a) for $\tau =0.1$ , (b) for $\tau ={10}^{-2}$ , (c) for $\tau ={10}^{-3}$ .

Example 6.4. Consider the parabolic problem (1.3) with $\kappa =1$ and $\left(x,t\right)\in \left(-1,1\right)\times \left(0,1\right]$ . Assume that the exact solution and Darcy’s equation are given as follows

$U\left(x,t\right)=\sin \left(\pi x\right){\text{e}}^{-\frac{3}{2}{\pi }^{2}t},$

$P\left(x,t\right)=\pi \cos \left(\pi x\right){\text{e}}^{-\frac{3}{2}{\pi }^{2}t}.$

In this example, we apply our scheme (2.4) to solve the numerical solution of Example 6.4 and compare our scheme with the Legendre Galerkin Chebyshev collocation least squares method (LGCC-LS) proposed in [43].

In this example, we apply our scheme (2.4) to solve the numerical solution of Example 6.4 and compare our scheme with the Legendre Galerkin Chebyshev collocation least squares method (LGCC-LS) proposed in [43]. The images of the exact solution $U$ and its approximation solution $u_N$ are shown in Figure 7, where $N=24$ and $\tau ={10}^{-6}$ , and its corresponding error function at $T=1$ is also given. Similarly, Figure 8 also shows that our scheme preserves Darcy’s law. The $L^2$ -errors derived by our scheme (2.4) at $T=1$ are shown in Table 5, where $N$ is increasing from 8 to 20 for $\tau ={10}^{-3}$ and $\tau ={10}^{-5}$ . The $L^2$ -errors obtained from the LGCC-LS method with $T=1$ are also listed in Table 5.

Example 6.5. Consider the two-dimensional parabolic Equation (1.3) with $\kappa =1$ and $\left(x,y,t\right)\in {\left(-1,1\right)}^2\times \left(0,T\right]$ . Assume that the exact solution $U$ and Darcy’s equation are given as

$U\left(x,y,t\right)=\text{exp}\left(\text{sin}\left(t\right)\right)\text{cos}\left(\frac{\pi }{2}x\right)\text{cos}\left(\frac{\pi }{2}y\right),$

$P\left(x,y,t\right)=\left(\begin{array}{l}{P}^{\left(1\right)}\left(x,y,t\right)\hfill \\ P^{\left(2\right)}\left(x,y,t\right)\hfill \end{array}\right)=\left(\begin{array}{l}-\frac{\pi }{2}\exp \left(\sin \left(t\right)\right)\sin \left(\frac{\pi }{2}x\right)\cos \left(\frac{\pi }{2}y\right)\hfill \\ -\frac{\pi }{2}\exp \left(\sin \left(t\right)\right)\cos \left(\frac{\pi }{2}x\right)\sin \left(\frac{\pi }{2}y\right)\hfill \end{array}\right).$

Table 5. The $L^2$ -errors of Scheme (2.4) and the LGCC-LS method with $\tau ={10}^{-3}$ and $\tau ={10}^{-5}$ at $T=1$ in Example 6.4.

$\tau$	$N$	Scheme (2.4)		LGCC-LS [43]
		$\Vert E_u\Vert$	$\Vert E_p\Vert$	$\Vert E_u\Vert$	$\Vert E_p\Vert$
10−3	8	2.6520e−09	1.2201e−02	2.7450e−09	3.9857e−03
	12	2.5863e−09	7.2884e−06	3.0428e−09	3.2813e−06
	16	2.6220e−09	1.8153e−08	3.5135e−09	1.4404e−08
	20	2.6443e−09	1.7952e−08	3.9282e−09	1.6058e−08
10−5	8	1.3944e−10	1.2201e−02	2.7219e−10	3.9857e−03
	12	2.5926e−13	7.2785e−06	3.0719e−13	3.2842e−06
	16	2.6213e−13	1.1010e−09	3.5132e−13	6.2322e−10
	20	2.6434e−13	1.9048e−12	3.9270e−13	1.0826e−11

(a)

(b)

Figure 7. The plot of solutions and the corresponding error function with $T=1$ in Example 6.4, where $N=24$ and $\tau ={10}^{-6}$ , (a) The image of the exact solution $U\left(T\right)$ and the approximate solution $u_N\left(T\right)$ , For the exact solution $U$ and its numerical solution $u_N$ ; (b) The image of the error function $E_u$ , For the error function $E_u$ .

Figure 8. $L^2$ -error evolution for Darcy’s flow with varying time steps $\tau$ in Example 6.4, where $N=24$ . (a) for $\tau =0.1$ , (b) for $\tau ={10}^{-2}$ , (c) for $\tau ={10}^{-3}$ .

Here, $U_0$ and $f$ are obtained from above exact solution, respectively. In this example, we test the effectiveness of our scheme (2.4) for solving the two-dimensional parabolic Equation (1.3).

Figure 9 shows the images of the approximate solution $u_N$ and the error between the exact solution $U$ and its numerical solution $u_N$ , where $N=14$ and $\tau ={10}^{-5}$ . Table 6 lists errors in the $L^2$ -norm of our scheme with increasing $N$ from 8 to 14, where $\tau ={10}^{-5}$ . These numerical results verify that our scheme is effective for solving two-dimensional parabolic equation.

(a)

(b)

Figure 9. The images of the exact solution $U$ and its numerical solution $u_N$ , and the corresponding error function $E_u$ for Example 6.5, where $T=1$ , $N=14$ , and $\tau ={10}^{-5}$ . (a) For the exact solution $U$ and it approximate solution $u_N$ , (b) For the error function $E_u$ .

Table 6. Errors in the $L^2$ -norm at $T=1$ of our scheme with $\tau ={10}^{-5}$ for Example 6.5.

$N$	$\tau$	$\Vert E_u\Vert$	$\Vert E_{p_1}\Vert$	$\Vert E_{p_2}\Vert$
8	10−5	1.6656e−04	0.0020	0.0020
10		4.4773e−06	4.4773e−05	4.4773e−05
12		4.4801e−06	6.0700e−05	6.0700e−05
14		4.6018e−06	6.2208e−05	6.2208e−05

In order to verify Darcy’s flow equation, we consider the error graphs of each component of $p_N$ when $T=10$ and $N=32$ respectively. Among them, Figure 10 corresponds to the time step $\tau =0.1$ , Figure 11 corresponds to $\tau ={10}^{-2}$ , and Figure 12 corresponds to $\tau ={10}^{-3}$ . These numerical results show our scheme preserves Darcy’s equation.

Figure 10. Images of error functions for each component of Darcy’s law, where $N=32$ , $T=10$ , and $\tau =0.1$ . (a) For the image of the error function $E\left(p_N^{\left(1\right)}\right)$ , (b) For the image of the error function $E\left(p_N^{\left(2\right)}\right)$ .

Figure 11. Graphs of error functions for each component of Darcy’s law equation, where $N=32$ , $T=10$ , and $\tau ={10}^{-2}$ . (a) For the image of the error function $E\left(p_N^{\left(1\right)}\right)$ , (b) For the image of the error function $E\left(p_N^{\left(2\right)}\right)$ .

Figure 12. Images of error functions for each component of Darcy’s law, where $N=32$ , $T=10$ , and $\tau ={10}^{-3}$ . (a) For the image of the error function $E\left(p_N^{\left(1\right)}\right)$ , (b) For the image of the error function $E\left(p_N^{\left(2\right)}\right)$ .

7. Conclusion

In this paper, the Legendre-Petrov-Galerkin method has been developed for solving the parabolic equation satisfying Darcy’s law. Consequently, the proposed scheme can derive the approximate solutions of the original solution $U$ and its flux variable $P$ simultaneously. The stability analysis of the semi-discrete scheme has been given by using Gronwall’s inequality (integral form) and Darcy’s equation and the corresponding error estimate has been derived. Some numerical results of one-dimensional and two-dimensional problems are presented to test the high order precision and the effectiveness of our scheme. In future work, we shall give an extension of our scheme for solving the nonlinear model that satisfies Darcy’s equation and its optimal error estimate.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (No.12561067), the Science and Technology Project of Guangxi (No. GuikeAD25069086), Innovation Project of GUET Graduate Education (No. 2025YCXS127), and the Guangxi Key Laboratory of Automatic Detecting Technology and Instruments (YQ22106).