1.
Introduction
In this paper, a nonlinear compact finite difference scheme is studied for the following the pseudo-parabolic Burgers' equation [1]
subject to the periodic boundary condition
and the initial data
where μ>0 is the coefficient of kinematic viscosity, γ and ε>0 are two parameters, φ(x) is an L-periodic function. Parameter L denotes the spatial period. Setting ε=0, Eq (1.1) reduces to a viscous Burgers' equation [2]. Equation (1.1) is derived by the degenerate pseudo-parabolic equation [3]
where α, β, κ, γ are nonnegative constants. The derivative term {uκ(uγut)x}x represents a dynamic capillary pressure relation instead of a usual static one [4]. Equation (1.4) is a model of one-dimensional unsaturated groundwater flow.
Here, u denotes the water saturation. We refer to [5] for a detailed explanation of the model. Equation (1.1) is also viewed as a simplified edition of the Benjamin-Bona-Mahony-Burgers (BBM-Burgers) equation, or a viscous regularization of the original BBM model for the long wave propagation [6]. The problem (1.1)–(1.3) has the following conservation laws
Based on (1.6), by a simple calculation, the exact solution satisfies
where c0=(1+√L2)√E(0).
Numerical and theoretical research for solving (1.1)–(1.3) have been extensively carried out. For instance, Koroche [7] employed the the upwind approach and Lax-Friedrichs to obtain the solution of In-thick Burgers' equation. Rashid et al. [8] employed the Chebyshev-Legendre pseudo-spectral method for solving coupled viscous Burgers' equations, and the leapfrog scheme was used in time direction. Qiu al. [9] constructed the fifth-order weighted essentially non-oscillatory schemes based on Hermite polynomials for solving one dimensional non-linear hyperbolic conservation law systems and presented numerical experiments for the two dimensional Burgers' equation. Lara et al. [10] proposed accelerate high order discontinuous Galerkin methods using Neural Networks. The methodology and bounds are examined for a variety of meshes, polynomial orders, and viscosity values for the 1D viscous Burgers' equation. Pavani et al. [11] used the natural transform decomposition method to obtain the analytical solution of the time fractional BBM-Burger equation. Li et al. [12] established and proved the existence of global weak solutions for a generalized BBM-Burgers equation. Wang et al. [13] introduced a linearized second-order energy-stable fully discrete scheme and a super convergence analysis for the nonlinear BBM-Burgers equation by the finite element method. Mohebbi et al. [14] investigated the solitary wave solution of nonlinear BBM-Burgers equation by a high order linear finite difference scheme.
Zhang et al. [15] developed a linearized fourth-order conservative compact scheme for the BBMB-Burgers' equation. Shi et al. [16] investigated a time two-grid algorithm to get the numerical solution of nonlinear generalized viscous Burgers' equation. Li et al. [17] used the backward Euler method and a semi-discrete approach to approximate the Burgers-type equation. Mao et al. [18] derived a fourth-order compact difference schemes for Rosenau equation by the double reduction order method and the bilinear compact operator. It offers an effective method for solving nonlinear equations. Cuesta et al. [19] analyzed the boundary value problem and long-time behavior of the pseudo-parabolic Burgers' equation. Wang et al. [20] proposed fourth-order three-point compact operator for the nonlinear convection term. They adopted the classical viscous Burgers' equation as an example and established the conservative fourth-order implicit compact difference scheme based on the reduction order method. The compact difference scheme enables higher accuracy in solving equations with fewer grid points. Therefore, using the compact operators to construct high-order schemes has received increasing attention and application [21,22,23,24,25,26,27,28,29].
Numerical solutions for the pseudo-parabolic equations have garnered widespread attention. For instance, Benabbes et al. [30] provided the theoretical analysis of an inverse problem governed by a time-fractional pseudo-parabolic equation. Moreover, Ilhan et al. [31] constructed a family of travelling wave solutions for obtaining hyperbolic function solutions. Di et al. [32] established the well-posedness of the regularized solution and gave the error estimate for the nonlinear fractional pseudo-parabolic equation. Nghia et al. [33] considered the pseudo-parabolic equation with Caputo-Fabrizio fractional derivative and gave the formula of mild solution. Abreu et al. [34] derived the error estimates for the nonlinear pseudo-parabolic equations basedon Jacobi polynomials. Jayachandran et al. [35] adopted the Faedo-Galerkin method to the pseudo-parabolic partial differential equation with logarithmic nonlinearity, and they analyzed the global existence and blowup of solutions.
To the best of our knowledge, the study of high-order difference schemes for Eq (1.1) is scarce. The main challenge is the treatment of the nonlinear term uux, as well as the error estimation of the numerical scheme. Inspired by the researchers in [15] and [20], we construct an implicit compact difference scheme based on the three–point fourth-order compact operator for the pseudo-parabolic Burgers' equation. The main contribution of this paper is summarized as follows:
● A fourth-order compact difference scheme is derived for the pseudo-parabolic Burgers' equation.
● The pointwise error estimate (L∞-estimate) of a fourth-order compact difference scheme is proved by the energy method [36,37] for the pseudo-parabolic Burgers' equation.
● Numerical stability, unique solvability, and conservation are obtained for the high-order difference scheme of the pseudo-parabolic Burgers' equation.
In particular, our numerical scheme for the special cases reduces to several other ones in this existing paper (see e.g., [38,39]).
The remainder of the paper is organized as follows. In Section 2, we introduce the necessary notations and present some useful lemmas. A compact difference scheme is derived in Section 3 using the reduction order method and the recent proposed compact difference operator. In Section 4, we establish the key results of the paper, including the conservation invariants, boundedness, uniqueness of the solution, stability, and convergence of the scheme. In Section 4.4, we present several numerical experiments to validate the theoretical findings, followed by a conclusion in Section 5.
Throughout the paper, we assume that the exact solution u(x,t) satisfies u(x,t)∈C6,3([0,L]×[0,T]).
2.
Notations and lemmas
In this section, we introduce some essential notations and lemmas. We begin by dividing the domain [0,L]×[0,T]. For two given positive integers, M and N, let h=L/M,τ=T/N. Additionally, denote xi=ih,0≤i≤M,tk=kτ,0≤k≤N; Vh={v|v={vi},vi+M=vi}. For any grid function u, v∈Vh, we introduce
Moreover, we introduce the discrete inner products and norms (semi-norm)
The following lemmas play important roles in the numerical analysis later, and we collect them here.
Lemma 1. [15,40] For any grid functions u, v∈Vh, we have
Lemma 2. [40] For any grid function v∈Vh and arbitrary ξ>0, we have
Lemma 3. [20] Let g(x)∈C5[xi−1,xi+1] and G(x)=g″(x), we have
Lemma 4. [15,18] For any grid functions u, v∈Vh and S∈Vh satisfying
we have the following results:
(I)
(II)
Proof. The result in (2.2)–(2.3) has been described in [15], and (2.5) has been proven in [18], we only need to only prove (2.4) and (2.6). Using the definition of the operator, we have
Taking the inner product of (2.1) with vk+12, we have
Therefore, the result (2.6) is obtained. □
Remark 1. [18] Denote 1=(1,1,⋯,1)T∈Vh. If S=0 in (2.1), then we further have
3.
Construction of the compact difference scheme
Let v=uxx, then the problem (1.1) is equivalent to
According to (3.2) and (3.4), it is easy to know that
Define the grid functions U={Uki|1≤i≤M,0≤k≤N} with Uki=u(xi,tk), V={Vki|1≤i≤M,0≤k≤N} with Vki=v(xi,tk). Considering (3.1) at the point (xi,tk+12) and (3.2) at the point (xi,tk), respectively, we have
Using the Taylor expansion and Lemma 3, we have
Noticing the initial-boundary value conditions (3.3)–(3.5), we have
There is a positive constant c1 such that the local truncation errors satisfy
Omitting the local truncation error terms in (3.6) and combining them with (3.7) and (3.8), the difference scheme for (3.1)–(3.5) as follows
Remark 2. As we see from the difference equations (3.9) and (3.10), only three points for each of them are utilized to generate fourth-order accuracy for the nonlinear pseudo-parabolic Burgers' equation without using additional boundary message. This is the reason we call this scheme the compact difference scheme. In addition, a fast iterative algorithm can be constructed, as shown in the numerical part in Section 4.4.
4.
Numerical analysis
4.1. Conservation and boundedness
Theorem 1. Let {uki,vki|1≤i≤M,0≤k≤N} be the solution of (3.9)–(3.12). Denote
Then, we have
Proof. Taking an inner product of (3.9) with 1, we have
By using Remark 1 in Lemma 4, the equality above deduces to
namely
□
Theorem 2. Let {uki,vki|1≤i≤M,0≤k≤N} be the solution of (3.9)–(3.12). Then it holds that
where
Proof. Taking the inner product of (3.9) with uk+12, and applying Lemma 1, we have
With the help of (2.2) and (2.4) in Lemma 4, the equality above deduces to
Replacing the superscript k with l and summing over l from 0 to k−1, we have
which implies that
□
Remark 3. Combining Lemma 1 with Theorem 2, it is easy to know that there is a positive constant c2 such that
4.2. Existence and uniqueness
Next, we recall the Browder theorem and consider the unique solvability of (3.9)–(3.12).
Lemma 5 (Browder theorem[41]). Let (H,(⋅,⋅)) be a finite dimensional inner product space, ‖⋅‖ be the associated norm, and Π:H→H be a continuous operator. Assume
Then there exists a z∗∈H satisfying ‖z∗‖≤α such that Π(z∗)=0.
Theorem 3. The difference scheme (3.9)–(3.12) has a solution at least.
Proof. Denote
It is easy to know that u0 has been determined by (3.11). From (3.10) and (3.11), we can get v0 by computing a system of linear equations as its coefficient matrix is strictly diagonally dominant. Suppose that {uk,vk} has been determined, then we may regard {uk+12,vk+12} as unknowns. Obviously,
Denote
Then the difference scheme (3.9)–(3.10) can be rewritten as
Define an operator Π on Vh:
Taking an inner product of Π(X) with X, we have
In combination of the technique from (2.2) in Lemma 4 and the Cauchy-Schwartz inequality, we have
and
Correspondingly,
Then
and
Substituting the equality above into (4.4) and according to (2.3) in Lemma 4, we have
Thus, when ‖X‖=αk, where αk=√‖uk‖2+ε22‖δxuk‖2+ε2h224‖vk‖2, then (Π(X),X)≥0. By Lemma 5, there exists a X∗∈Vh satisfying ‖X∗‖≤αk such that Π(X∗)=0. Consequently, the difference scheme (3.9)–(3.12) exists at least a solution uk+1=2X∗−uk. Observing, when (X∗1,X∗2,⋯,X∗M) is known, (Y∗1,Y∗2,⋯,Y∗M) can be determined by (4.3) uniquely. Thus, we know vk+1i=2Y∗i−vki, 1≤i≤M exists. □
Now we are going to verify the uniqueness of the solution of the difference scheme. We have the following result.
Theorem 4. When γ=0, the solution of the difference scheme (3.9)–(3.12) is uniquely solvable for any temporal step-size; When γ≠0 and τ≤min{4Lc2|γ|(L+1),2ε23c2|γ|(2L+1)}, the solution of the difference scheme (3.9)–(3.12) is uniquely solvable.
Proof. According to Theorem 3, we just need to prove that (4.2)–(4.3) has a unique solution. Suppose that both {u(1),v(1)}∈Vh and {u(2),v(2)}∈Vh are the solutions of (4.2)–(4.3), respectively. Let
Then we have
Taking an inner product of (4.5) with u, we have
With the application of Lemma 2 and (2.3) in Lemma 4, it follow from the equality above that
By the definition of ψ(⋅,⋅) and (4.1), we have
Using the Cauchy-Schwarz inequality, Lemmas 1 and 2, we have
Similarly, we have
Substituting (4.8) and (4.9) into (4.7), we can obtain
When τ≤min{4Lc2|γ|(L+1),2ε23c2|γ|(2L+1)}, we have ui=0,1≤i≤M. □
4.3. Convergence and stability
Let h0>0 and denote
and error functions
we have the following convergence results.
Theorem 5. Let {u(x,t),v(x,t)} be the solution of (3.1)–(3.5) and {uki,vki|0≤i≤M,0≤k≤N} be the solution of the difference scheme (3.9)–(3.12). When c4τ≤13 and h≤h0, we have
Proof. Subtracting (3.9)–(3.12) from (3.6)–(3.8), we can get an error system
Taking an inner product of (4.10) with ek+12, we have
Applying (2.3) in Lemma 4, we have
Similar to the derivation in (4.8) and (4.9), we have
and
Substituting (4.15)–(4.17) into (4.14), and using (2.3)–(2.4) in Lemma 4, we obtain
Then, we have
Using Cauchy-Schwartz inequality, we can rearrange the inequality above into the following form
Replacing the superscript k with l and summing over l from 0 to k, we get
For the last item on the right-hand side of (4.18), we have
Substituting (4.19) into (4.18) and using Lemma 2, we get
We can rearrange the inequality above into the following form
Denote
In combination of (2.6) in Lemma 4 with (4.12), we have
Substituting (4.21) into (4.20), when h≤h0 and c4τ≤13, (4.20) can be rewritten as
which implies that
According to the Gronwall inequality, we have
Thus, it holds that
□
Below, we consider the stability of the difference scheme (3.9)–(3.12). Suppose that {˜uki,˜vki|1≤i≤M,0≤k≤N} is the solution of
Denote
Subtracting (3.9)–(3.12) from (4.22), we obtain the perturbation equation as follows
Theorem 6. Let {ξki,ηki|1≤i≤M,0≤k≤N} be the solution of (4.23)–(4.26). When c4τ≤13 and h≤h0, we have
where c7=12√e32c4T(L2+15ε2).
Proof. Taking an inner product of (4.23) with ξk+12, we have
Similar to the analysis technique in Theorem 5, we obtain
Replacing the superscript k with l and summing over l from 0 to k, we get
Using Lemma 2 and when h≤h0, we can rearrange the inequality above into the following form
Denote
Combining (2.6) in Lemma 4 with (4.25), we have
Substituting (4.28) into (4.27), (4.27) can be rewritten as
Then, we have
According to the Gronwall inequality, when c4τ≤13, we have
where c7=12√e32c4T(L2+15ε2).
Therefore, it holds that
□
4.4. Numerical experiments
In this section, we perform numerical experiments to verify the effectiveness of the difference scheme and the accuracy of the theoretical results. Before conducting the experiments, we first introduce an algorithm for solving the nonlinear compact scheme. Denote
where 0≤k≤N. The algorithm of the compact difference scheme (3.9)–(3.12) can be described as follows:
Step 1 Solve u0 and v0 based on (3.10) and (3.11).
Step 2 Suppose uk is known, the following linear system of equations will be used to approximate the solution of the difference scheme (3.9)–(3.12), for 1≤i≤M, we have
until
Let
In the following numerical experiments, we set the tolerance error ϵ=1×10−12 for each iteration unless otherwise specified.
When the exact solution is known, we define the discrete error in the L∞-norm as follows:
where Uki and uki represent the analytical solution and the numerical solution, respectively. Additionally, the convergence orders in space and time are defined as follows:
When the exact solution is unknown, we use a posteriori error estimation to verify the convergence orders in space and time. For a sufficiently small h, we denote
Similarly, for sufficiently small τ, we denote
Example 1. We first consider the following equation
where
The initial condition is determined by the exact solution u(x,t)=etsin(πx) with the period L=2 and T=1.
The numerical results are reported in Table 1 and Figure 1.
Table 1, we progressively reduce the spatial step-size h half by half (h=1/2,1/4,1/8,1/16,1/32) while keeping the time step-size τ=1/1000. Conversely, we gradually decrease the time step-size τ half by half (τ=1/4,1/8,1/16,1/32,1/64) while maintaining the spatial step-size h=1/50.
As we can see, the spatial convergence order approaches to the four order approximately, and the temporal convergence order approaches to two orders in the maximum norm, which are consistent with our convergence results. Comparing our numerical results with those in [42] from Table 2, we find our scheme is more efficient and accurate.
By observing the first subgraph in Figure 1, the evolutionary trend surface of the numerical solution u(x,t) with τ=1/1000, h=1/50, L=2 and T=1 is illustrated. This figure successfully reflects the panorama of the exact solution. In order to verify the accuracy of the difference scheme (3.9)–(3.12), we have drawn the numerical error surface in the second subgraph in Figure 1 with τ=1/1000, h=1/50, L=2 and T=1.
We observe that the rates of the numerical error in the maximum norm approaches a fixed value, which verifies that the difference scheme (3.9)–(3.12) is convergent. It took us 2.03 seconds to compute the spatial order of accuracy and 0.37 seconds to determine the temporal order of accuracy.
Example 2. We further consider the problem of the form
where the exact solution is unavailable.
Case I ε=1:
The numerical results are reported in Table 3 and Figure 2. The two discrete conservation laws of the difference scheme (3.9)–(3.12) are reported in Table 4. In the following calculations, we set T=1. First, we fix the temporal step-size τ=1/1000 and reduce the spatial step-size h half by half (h=50/11,50/22,50/44,50/88). Second, we fix the spatial step-size h=1/2, meanwhile, reduce the temporal-step size τ half by half (τ=1/2,1/4,1/8,1/16).
As we can see, the spatial convergence order approaches to four orders, approximately, and the temporal convergence order approaches to two orders in the maximum norm, which is consistent with our convergence results. It took us 6.74 seconds to compute the spatial order of accuracy, and 0.30 seconds to determine the temporal order of accuracy.
From Table 4, we can see that the discrete conservation laws in Theorems 1 and 2 are also satisfied. In the first graph of Figure 2, we depict the evolutionary trend surface of the numerical solution u(x,t) with τ=1/1000, h=1/2, L=50 and T=1, and this figure successfully reflects the panorama of the exact solution.
When simulating a short duration of time T=1, the impact of values ε=1 and ε=0.1 on the numerical simulation is relatively small. Therefore, in the following Case Ⅱ, we take ε=0.1 and T=10 to observe the impact of ε on the numerical simulation situation.
Case II ε=0.1:
The numerical results are reported in Table 5 and Figure 2. The two discrete conservation laws of the difference scheme (3.9)–(3.12) are reported in Table 6.
First, we fix the temporal step-size τ=1/1000 and reduce the spatial step-size h half by half (h=25/3,25/6,25/12,25/24). Second, we fix the spatial step-size h=1/2 and reduce the temporal step size τ half by half (τ=1/2,1/4,1/8,1/16).
As we can see, the spatial convergence order approaches to four orders, approximately, and the temporal convergence order approaches to two orders in the maximum norm, which is consistent with our convergence results. It took us 33.76 seconds to compute the spatial order of accuracy and 0.33 seconds to determine the temporal order of accuracy.
In the second subgraph of Figure 2, we depict the evolutionary trend surface of the numerical solution u(x,t) with τ=1/1000, h=1/2, L=50 and T=10. Compared with the first subgraph of Figure 2, smaller ε amplifies sharper transitions and wave-like behavior, whereas the larger ε makes the solution smoother.
Example 3. In the last example, we consider the problem
with the Maxwell initial conditions u(x,0)=e−(x−7)2,0≤x≤30.
Figure 3 reflects the behavior of the solutions to the pseudo-parabolic Burgers' equation. During the propagation process, we observe that the pseudo-parabolic Burgers' equation exhibits characteristics of both diffusion and advection. As we can see, the peak gradually spreads out and flattens as time progresses. Additionally, the solution moves to the right, indicating propagation direction.
The numerical scheme ensures stability, convergence, and the preservation of physical properties, which can be observed from the smooth transitions over time. This phenomenon indicates that the numerical scheme preserves the physical properties, ensuring stability and convergence.
In Figure 4, we observe that the pseudo-parabolic Burgers' equation exhibits propagation characteristics coupled with gradual damping.
From Table 7, we can see that the discrete conservation law agrees well with Theorems 1 and 2. The value of Q remains almost constant throughout the simulation, which is crucial for maintaining the physical integrity of the solution. Similarly, the phenomenon is suitable for the energy E, and these results further verify the correctness and reliability of the high-order compact difference scheme.
5.
Conclusions
We propose and analyze an implicit compact difference scheme for the pseudo-parabolic Burgers' equation, achieving second-order accuracy in time and fourth-order accuracy in space. Using the energy method, we provide a rigorous numerical analysis of the scheme, proving the existence, uniqueness, uniform boundedness, convergence, and stability of its solution. Finally, the theoretical results are validated through numerical experiments. The experimental results demonstrate that the proposed scheme is highly accurate and effective, aligning with the theoretical predictions. As part of our ongoing research [42,43,44,45,46,47,48], we aim to extend these techniques and approaches to other nonlocal or nonlinear evolution equations [49,50,51,52,53,54,55].
Use of AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
Acknowledgments
The authors would like to thank the supervisor Qifeng Zhang, who provide this interesting topic and detailed guidance. They are also guilt to Baohui Xie's work when he studied in Zhejiang Sci-Tech Universtiy. The authors are grateful to the editor and the anonymous reviewers for their careful reading and many patient checking of the whole manuscript.
The work is supported by Institute-level Project of Zhejiang Provincial Architectural Design Institute Co., Ltd. (Research on Digital Monitoring Technology for Central Air Conditioning Systems, Institute document No. 18), and National Social Science Fund of China (Grant No. 23BJY006) and General Natural Science Foundation of Xinjiang Institute of Technology(Grant No. ZY202403).
Conflict of interest
The authors declare there is no conflicts of interest.