Numerical Solutions of the Nonlinear Dispersive Shallow Water Wave Equations Based on the Space–Time Coupled Generalized Finite Difference Scheme

Abstract

This study applies the space–time generalized finite difference scheme to solve nonlinear dispersive shallow water waves described by the modified Camassa–Holm equation, the modified Degasperis–Procesi equation, the Fornberg–Whitham equation, and its modified form. The proposed meshless numerical scheme combines the space–time generalized finite difference method, the two-step Newton’s method, and the time-marching method. The space–time approach treats the temporal derivative as a spatial derivative. This enables the discretization of all partial derivatives using a spatial discretization method and efficiently handles mixed derivatives with the proposed mesh-less numerical scheme. The space–time generalized finite difference method is derived from Taylor series expansion and the moving least-squares method. The numerical discretization process only involves functional data and weighting coefficients on the central and neighboring nodes. This results in a sparse matrix system of nonlinear algebraic equations that can be efficiently solved using the two-step Newton’s method. Additionally, the time-marching method is employed to advance the space–time domain along the time axis. Several numerical examples are presented to validate the effectiveness of the proposed space–time generalized finite difference scheme.

1. Introduction

Dispersive shallow water waves are wave propagation phenomena that have received extensive attention in engineering due to their numerous applications in physics, engineering, mechanics, and mathematics. In the past few decades, several partial differential equations (PDEs), including the Camassa−Holm equation [1,2], the Degasperis−Procesi equation [2,3,4], the Fornberg−Whitham equation [3,5,6,7,8,9,10,11,12], and the Schrödinger equation [13,14,15,16,17], have been presented to describe the traveling wave phenomenon and water waves of finite depths. These PDEs and their modified forms mathematically represent transient nonlinear PDEs that involve third-order time–space mixed partial derivatives or complex terms. Therefore, obtaining the numerical solutions of those PDEs is a significant issue for analyzing the traveling wave phenomenon. In previous studies, two types of methods were mainly used to obtain solutions for such equations: the approximate method, the spectral method [18,19], and numerical simulation.

Due to the complexity of these equations, finding exact solutions is generally challenging. Therefore, numerical methods play a crucial role in solving these equations. Traditional numerical approaches, relying on temporal and spatial discretization, face difficulties when dealing with high-order mixed partial derivatives. As a result, researchers have extensively employed approximate methods to derive numerical solutions, which offer convenient and effective alternatives to classical numerical treatments. For example, Wazwaz [2] used the extended tanh method to establish the new solitary wave solutions of the modified Degasperis−Procesi and Camassa−Holm equations. Abidi and Omrani [10] applied the homotopy analysis method (HAM) to acquire the HAM results of the Fornberg−Whitham equation and compared them with Adomian’s decomposition method (ADM). Ahmad et al. [5] utilized the modified variational iteration algorithm (MVIA) to obtain the approximate numerical solutions of the Fornberg−Whitham equation and its modified form. In addition to the above-mentioned methods, there are other approximate methods such as the method of phase portrait analysis [9], He’s variational iteration method (VIM) [8], and the reproducing kernel Hilbert space method (RKHSM) [7]. Recently, due to the evolution of computer technology and numerical methods, novel numerical schemes have also been applied to solve such equations. These include the Godunov method [6], the cubic B-spline quasi-interpolation method [4], the hybrid radial basis function method [3], the linearized implicit pseudo-spectral method [18], and the Gegenbauer spectral method [19].

In the field of numerical discretization methods for PDEs, two main categories can be distinguished: mesh-based methods and meshless/free methods. Mesh-based methods include well-known techniques such as finite difference (FDM), finite volume, and finite element methods. These methods have been extensively developed and modified to solve engineering problems and analyze complex physical phenomena that are challenging to investigate through practical experiments. Despite its widespread application, mesh generation remains a challenging task. To overcome this issue, the concept of using numerical discretization methods without the need for a mesh has emerged as an important research direction. Meshless methods offer a solution by avoiding the requirement for mesh generation, making it easier to handle computational domains with complex shapes. In meshless methods, numerical discretization involves constructing interpolation functions on arbitrary distribution nodes to approximate the governing equation. Over the years, meshless methods have gained advantages such as simplified numerical procedures, easy implementation, flexibility, and the ability to construct hybrid numerical schemes tailored to specific research needs. At present, there are several meshless numerical methods, for example, the radial basis functions collocation method (RBFCM) [3], the virtual boundary meshless Galerkin method [20], the moving particle semi-implicit method [21], the radial point interpolation meshless method [22], the singular boundary method [23], the localized scheme based on boundary-type method [24,25,26], the generalized finite difference method (GFDM) [27,28,29,30,31,32,33,34,35,36,37], etc.

The GFDM is a localized domain-type meshless collocation method based on the multivariate Taylor series expansion theory and the moving least-squares method. In 1980, the idea of FDM without mesh was proposed by Liszka and Orkisz [27]. After Benito et al. [28] proposed the explicit numerical process in 2001, numerical applications related to the GFDM have been proposed in the past two decades. By utilizing the GFDM, the approximate expression of the partial derivative terms can be represented by a linear combination of the weight coefficients and function values in a local support domain. This support domain only contains the goal node and its nearby nodes. Once the nodes are distributed in the computational domain, the governing equation can be discretized into a system of algebraic equations. In addition, the resultant matrix system is a sparse matrix due to the localized property of the GFDM. The above characteristics make the GFDM capable of applying to complex problems, such as theoretical analysis [29], shallow water equation [30,31], porous media flow [32], wave propagation [33], elliptic interface problems [34], extended Fisher–Kolmogorov equation [35], stream function formulation [36], and elastic wave [37]. In this research, the space–time GFDM (ST-GFDM) was applied as the foundation of the proposed meshless numerical scheme, and it is an extended meshless method from the GFDM.

The space–time (ST) coupled approach is a numerical technique that treats the temporal derivative as one of the spatial derivatives. Prior to the use of the ST-coupled approach, time-dependent problems were typically solved by combining spatial discretization methods with temporal discretization methods. However, this approach often led to challenges in defining the accuracy and stability of the numerical scheme. The ST-coupled approach transforms a time-dependent problem of dimension (n) into a steady-state problem of dimension (n + 1). As a result, all the partial derivatives in the governing equation can be discretized using numerical methods designed for spatial discretization. This simplification enables a clearer determination of the properties of the numerical scheme, such as accuracy and stability. Recently, the ST coupled approach has been widely combined with meshless methods such as the ST localized RBFCM [38], ST kernel-based method [39], ST localized method of fundamental solutions [40,41], ST Trefftz Method [42,43], ST backward substitution method [44], and ST-GFDM [45,46,47,48,49,50,51,52,53,54]. Based on the flexibility of the ST-GFDM, researchers have applied this meshless method for engineering problems in the past few years. Examples of these applications include heat conduction problems [45,54], Burgers’ equations [46], parabolic PDEs [47], unsteady double-diffusive natural convection [48], thin elastic plate bending problems [49], thermoelasticity problems [50], equal-width equation [51], Zakharov–Kuznetsov-Modified equal-width equation [52], and hyperbolic PDEs [53], etc. The aforementioned advancements demonstrate that ST-GFDM is a promising meshless method with potential for engineering applications. In this study, the ST-GFDM is utilized in combination with the two-step Newton’s method and the time-marching approach to tackle the problem of solving dispersive shallow water waves.

After providing the motivation for the research. The governing equations of the dispersive shallow water wave are represented in Section 2. In the section on proposed ST-GFDM, the numerical process is described, and four numerical examples are simulated to verify the proposed method in Section 4. Then, the conclusions are represented in the final section. All abbreviations and descriptions in this paper are listed in

Table 1. The abbreviations and descriptions in this paper.

Abbreviations	Description
1st	first
1D	one-dimensional
2D	two-dimensional
2nd	second
3rd	third
ADM	Adomian’s decomposition method
Eq	equation
FDM-RK4	finite difference method with the fourth-order Runge–Kutta method
GFDM	generalized finite difference method
HAM	homotopy analysis method
MAE	maximum absolute error
MVIA	modified variational iteration algorithm
PDE	partial differential equation
RBFCM	radial basis functions collocation method
RKHSM	reproducing kernel Hilbert space method
RMSE	root-mean-squared error
ST	space–time
ST-GFDM	space–time GFDM
VIM	He’s variational iteration method
Superscript
ext	exact solution
num	numerical solution

for quick and easy reference for the reader.

2. The Mathematical Models of the Dispersive Shallow Water Waves

The dispersive shallow water wave represents the behavior of the wave propagation, breaking wave, and dynamic system of the fluid. It is an important phenomenon in physics and engineering. In the past few decades, several researchers have proposed the mathematical model of the dispersive shallow water wave problem, and its general form of the governing equation [1,2,11] is described as,

$$\frac{\partial u}{\partial t}-\frac{{\partial }^{3}u}{\partial x^2\partial t}+\left(\alpha +\beta u^{\omega }\right)\frac{\partial u}{\partial x}=q\frac{\partial u}{\partial x}\frac{{\partial }^{2}u}{\partial x^2}+u\frac{{\partial }^{3}u}{\partial x^3}, a\le x\le b, t\ge 0,$$

(1)

where u is the wave speed; x and t are the spatial- and temporal-axes; α, β, ω, and q are the constants. The initial conditions and boundary conditions are represented as follows:

$$\{\begin{array}{c}u\left(x,0\right)=U_0\left(x\right)\\ u\left(a,t\right)=f_1\left(t\right)\\ u\left(b,t\right)=f_2\left(t\right)\end{array}.$$

(2)

By applying the constants a, b, ω, and q, three different types of the dispersive shallow water wave equation are defined. The first one is the Fornberg–Whitham equation [5,6,7,8,11] that is given with α = 1, β = 1, ω= 1, and q = 3 as,

$$\frac{\partial u}{\partial t}-\frac{{\partial }^{3}u}{\partial x^2\partial t}+\left(1+u\right)\frac{\partial u}{\partial x}=3\frac{\partial u}{\partial x}\frac{{\partial }^{2}u}{\partial x^2}+u\frac{{\partial }^{3}u}{\partial x^3},a\le x\le b, t\ge 0,$$

(3)

and its modified form [5,7,10] is given (α = 1, β = 1, k = 2, and q = 3) as,

$$\frac{\partial u}{\partial t}-\frac{{\partial }^{3}u}{\partial x^2\partial t}+\left(1+u^2\right)\frac{\partial u}{\partial x}=3\frac{\partial u}{\partial x}\frac{{\partial }^{2}u}{\partial x^2}+u\frac{{\partial }^{3}u}{\partial x^3}.$$

(4)

The remaining two are similar: the Camassa–Holm equation [1,2] and the Degasperis–Procesi equation [2,3,4], respectively. By giving α = 0, β = 3, ω = 1, and q = 2, the Camassa–Holm equation can be denoted as

$$\frac{\partial u}{\partial t}-\frac{{\partial }^{3}u}{\partial x^2\partial t}+3u\frac{\partial u}{\partial x}=2\frac{\partial u}{\partial x}\frac{{\partial }^{2}u}{\partial x^2}+u\frac{{\partial }^{3}u}{\partial x^3},$$

(5)

and the Degasperis–Procesi equation (α = 0, β = 4, ω = 0, and q = 3) is described as

$$\frac{\partial u}{\partial t}-\frac{{\partial }^{3}u}{\partial x^2\partial t}+4u\frac{\partial u}{\partial x}=3\frac{\partial u}{\partial x}\frac{{\partial }^{2}u}{\partial x^2}+u\frac{{\partial }^{3}u}{\partial x^3},$$

(6)

As with the relationship between the Fornberg–Whitham equation and its modified form, the modified forms of the Camassa–Holm equation [2] and the Degasperis–Procesi equation [2,3,4] are represented by changing the nonlinear term $u\frac{\partial u}{\partial x}$ into $u^2\frac{\partial u}{\partial x}$. Then, they can be described as follows:

$$\frac{\partial u}{\partial t}-\frac{{\partial }^{3}u}{\partial x^2\partial t}+3u^2\frac{\partial u}{\partial x}=2\frac{\partial u}{\partial x}\frac{{\partial }^{2}u}{\partial x^2}+u\frac{{\partial }^{3}u}{\partial x^3},$$

(7)

and

$$\frac{\partial u}{\partial t}-\frac{{\partial }^{3}u}{\partial x^2\partial t}+4u^2\frac{\partial u}{\partial x}=3\frac{\partial u}{\partial x}\frac{{\partial }^{2}u}{\partial x^2}+u\frac{{\partial }^{3}u}{\partial x^3},$$

(8)

Equations (1) and (3)–(7) are highly nonlinear PDEs and have 3rd-order mixed partial derivatives as $\frac{{\partial }^{3}u}{\partial x^2\partial t}$. The complexity of dispersive shallow water wave equations has posed challenges to solving them using original explicit or implicit numerical schemes. As a result, researchers have primarily focused on proposing approximate or exact solutions for these equations. Consequently, applying numerical schemes to such equations is still in its early stages. In this study, the proposed ST-GFDM scheme is employed to solve the Fornberg–Whitham equation and its modified form, the modified Camassa–Holm equation, and the modified Degasperis–Procesi equation. The next section provides an introduction to the numerical procedures involved in the proposed ST-GFDM scheme.

3. The Proposed Meshless Numerical Scheme

The numerical procedures of the proposed ST-GFDM scheme are described in this section. The proposed numerical scheme uses ST-GFDM for the numerical discretization to obtain the nonlinear algebraic equations. Then, the nonlinear algebraic system is solved by the two-step Newton’s method. A simple technique called the time-marching method is applied to save computational resources and solve numerical problems within an unpredictable end-time.

3.1. Space–Time Generalized Finite Difference Method

In this subsection, the numerical process of ST-GFDM is described. Once the time-dependent problem is determined, the distribution nodes can be set in the corresponding computational domain. For the one-dimensional (1D) time-dependent problem, the ST approach defines a two-dimensional (2D) steady-state problem in an x-t field. Since the ST approach is applied, the distribution nodes are set both in the space- and time-axes and named ST-domain [38,40,41], as shown in Figure 1a.

A supporting domain is set up by choosing the n_s nearest nodes within the central i node in Figure 1b (see the “×” symbol), and $\left(x_i,t_i\right), i=1,2,\dots ,N_T$ is defined as the coordinate of the ith node. On the other hand, $\left(x_{i,j},t_{i,j}\right), j=1,2,\dots ,n_s$ represent the coordinates of the nodes in the supporting domain. By applying the n-order Taylor expansion and the weighting function, a residual function $B_n\left(u_i\right)$ can be written as [36,51],

$$B_n\left(u_i\right)={\displaystyle \sum _{j=1}^{n_s}{\left[\left(-u_{i,j}+{\displaystyle \sum _{0\le o,m\le n, o+m\le n}^{}\frac{1}{o!m!}\cdot \frac{{\partial }^{o+m}{u}_i}{{\partial }^{o}x{\partial }^{m}t}\cdot h_{i,j}^o\cdot l_{i,j}^m+E_n}\right)w\left(h_{i,j},l_{i,j}\right)\right]}^2},$$

(9)

where $u_i$ and $u_{i,j}$ denotes the unknowns at the ith node and unknowns inside the ith supporting domain, respectively; $E_n$ is the truncation error of the n-order Taylor expansion; o and m is the natural number; $h_{i,j}=x_i-x_{i,j}$; $l_{i,j}=t_i-t_{i,j}$; and $w$ is the quartic spline function and describes as follows:

$$w\left(h_{i,j},l_{i,j}\right)=\{\begin{array}{l}1-6{\left(\frac{d_{i,j}}{d{m}_i}\right)}^2+8{\left(\frac{d_{i,j}}{d{m}_i}\right)}^3-3{\left(\frac{d_{i,j}}{d{m}_i}\right)}^4, d_{i,j}\le d{m}_i\\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}0 \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}, d_{i,j}>d{m}_i \hfill \end{array},$$

(10)

where $d_{i,j}$ is the distance between the ith node and each node inside the ith supporting domain, and $d_{i,j}=\sqrt{h_{i,j}^2+l_{i,j}^2}$; and $d{m}_i$ is the maximum distance.

In this study, the 3rd-order ST-GFDM was applied, and Equation (9) can be written as,

$$B_3\left(u_i\right)={\displaystyle \sum _{j=1}^{n_s}{\left[\left(\begin{array}{l}{u}_i-u_{i,j}+h_{i,j}{\frac{\partial u}{\partial x}|}_i+l_{i,j}{\frac{\partial u}{\partial t}|}_i+\frac{h_{i,j}^2}{2}{\frac{{\partial }^{2}u}{\partial x^2}|}_i+h_{i,j}{l}_{i,j}{\frac{{\partial }^{2}u}{\partial x\partial t}|}_i+\frac{l_{i,j}^2}{2}{\frac{{\partial }^{2}u}{\partial t^2}|}_i+\dots \\ \frac{h_{i,j}^3}{6}{\frac{{\partial }^{3}u}{\partial x^3}|}_i+{\frac{h_{i,j}^2{l}_{i,j}^{}}{2}\frac{{\partial }^{3}u}{\partial x^2\partial t}|}_i+\frac{h_{i,j}^{}{l}_{i,j}^2}{2}{\frac{{\partial }^{3}u}{\partial x\partial t^2}|}_i+\frac{l_{i,j}^3}{6}{\frac{{\partial }^{3}u}{\partial t^3}|}_i+E_3\end{array}\right)w_{i,j}\right]}^2},$$

(11)

The order of approximations in Equation (11), which accounts for the truncation errors, is determined by the Taylor series expansion up to the 3rd order. By minimalizing $B_3\left(u_i\right)$ with respect to,

$$D_{u_i}={\left[{\frac{\partial u}{\partial x}|}_i,{\frac{\partial u}{\partial t}|}_i,{\frac{{\partial }^{2}u}{\partial x^2}|}_i,{\frac{{\partial }^{2}u}{\partial x\partial t}|}_i,{\frac{{\partial }^{2}u}{\partial t^2}|}_i,{\frac{{\partial }^{3}u}{\partial x^3}|}_i,{\frac{{\partial }^{3}u}{\partial x^2\partial t}|}_i,{\frac{{\partial }^{3}u}{\partial x\partial t^2}|}_i,{\frac{{\partial }^{3}u}{\partial t^3}|}_i\right]}^{\mathrm{T}},$$

(12)

a linear system $A{D}_{u_i}=b{U}_i$ is obtained. In this linear system, the matrix A is denoted as $A=P^{\mathrm{{\rm T}}}{W}^{2}P,$ where

$$P={\left[\begin{array}{ccccccccc}{h}_{i,1}& l_{i,1}& \frac{h_{i,1}^2}{2}& h_{i,1}{l}_{i,1}& \frac{l_{i,1}^2}{2}& \frac{h_{i,1}^3}{6}& \frac{h_{i,1}^2{l}_{i,1}}{2}& \frac{h_{i,1}{l}_{i,1}^2}{2}& \frac{l_{i,1}^3}{6}\\ h_{i,2}& l_{i,2}& \frac{h_{i,2}^2}{2}& h_{i,2}{l}_{i,2}& \frac{l_{i,2}^2}{2}& \frac{h_{i,2}^3}{6}& \frac{h_{i,2}^2{l}_{i,2}}{2}& \frac{h_{i,2}{l}_{i,2}^2}{2}& \frac{l_{i,2}^3}{6}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ h_{i,j}& l_{i,j}& \frac{h_{i,j}^2}{2}& h_{i,j}{l}_{i,j}& \frac{l_{i,j}^2}{2}& \frac{h_{i,j}^3}{6}& \frac{h_{i,j}^2{l}_{i,j}}{2}& \frac{h_{i,j}{l}_{i,j}^2}{2}& \frac{l_{i,j}^3}{6}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ h_{i,n_s}& l_{i,n_s}& \frac{h_{i,n_s}^2}{2}& h_{i,n_s}{l}_{i,n_s}& \frac{l_{i,n_s}^2}{2}& \frac{h_{i,n_s}^3}{6}& \frac{h_{i,n_s}^2{l}_{i,n_s}}{2}& \frac{h_{i,n_s}{l}_{i,n_s}^2}{2}& \frac{l_{i,n_s}^3}{6}\end{array}\right]}_{n_s\times 9},$$

(13)

$W=diag\left(w_{i,1},w_{i,2},\cdots ,w_{i,n_s}\right)$; b is a coefficient matrix whose size is $9\times \left(n_s+1\right)$ and is written as $b=\left[\begin{array}{cc}{P}^{\mathrm{{\rm T}}}w& P^{\mathrm{{\rm T}}}{W}^2\end{array}\right]$; $U_i={\left[u_i,u_{i,1},u_{i,2},\dots ,u_{i,n_s}\right]}^{\mathrm{T}}$; and the superscript T is a symbol of transpose. By solving this linear system, the approximate expression of $D_{u_i}$ can be written as,

$$D_{u_i}=A^{-1}b{U}_i=E_i{U}_i,$$

(14)

where $E_i$ is the matrix of weighting coefficients within sized $9\times \left(n_s+1\right)$. Then, the derivatives are denoted as follows:

$$\{\begin{array}{c}{\frac{\partial u}{\partial x}|}_i=w{x}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}w{x}_{i,j}^{}{u}_{i,j}}\\ {\frac{\partial u}{\partial t}|}_i=w{t}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}w{t}_{i,j}^{}{u}_{i,j}}\\ {\frac{{\partial }^{2}u}{\partial x^2}|}_i=wx{x}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}wx{x}_{i,j}^{}{u}_{i,j}}\\ {\frac{{\partial }^{2}u}{\partial x\partial t}|}_i=wx{t}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}wx{t}_{i,j}^{}{u}_{i,j}}\\ {\frac{{\partial }^{2}u}{\partial t^2}|}_i=wt{t}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}wt{t}_{i,j}^{}{u}_{i,j}}\\ {\frac{{\partial }^{3}u}{\partial x^3}|}_i=wxx{x}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}wxx{x}_{i,j}^{}{u}_{i,j}}\\ {\frac{{\partial }^{3}u}{\partial x^2\partial t}|}_i=wxx{t}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}wxx{t}_{i,j}^{}{u}_{i,j}}\\ {\frac{{\partial }^{3}u}{\partial x\partial t^2}|}_i=wxt{t}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}wxt{t}_{i,j}^{}{u}_{i,j}}\\ {\frac{{\partial }^{3}u}{\partial t^3}|}_i=wtt{t}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}wtt{t}_{i,j}^{}{u}_{i,j}}\end{array},$$

(15)

where ${\left\{w{x}_{i,j}\right\}}_{j=0}^{n_s}$, ${\left\{w{t}_{i,j}\right\}}_{j=0}^{n_s}$, ${\left\{wx{x}_{i,j}\right\}}_{j=0}^{n_s}$, ${\left\{wx{t}_{i,j}\right\}}_{j=0}^{n_s}$, ${\left\{wt{t}_{i,j}\right\}}_{j=0}^{n_s}$, ${\left\{wxx{x}_{i,j}\right\}}_{j=0}^{n_s}$, ${\left\{wxx{t}_{i,j}\right\}}_{j=0}^{n_s}$,${\left\{wxt{t}_{i,j}\right\}}_{j=0}^{n_s}$, and ${\left\{wtt{t}_{i,j}\right\}}_{j=0}^{n_s}$ are the weighting coefficients, and obtained by solving Equation (14). Equation (15) is a significant part of the proposed ST-GFDM in that the partial derivatives are approximated by a linear combination of the weight coefficients and function values in a local support domain. Thus, the discretized governing equations and boundary conditions only relate to the central i node and its nearby n_s nodes in the local support domain. In addition, the numerical process of Equations (9)–(15) will apply for each node to obtain the weighting coefficients matrix $E_i$ of each node. Finally, the resultant matrix system, formed by the algebraic equations after discretization, is a sparse matrix.

The following contexts demonstrate the numerical discretization for the governing equation and boundary conditions to form the nonlinear algebraic systems. Since the simulated start, N_b1, N_b2, N_b3, and N_b4 nodes have been set up along the ${\mathrm{\Gamma }}_1$, ${\mathrm{\Gamma }}_2$, ${\mathrm{\Gamma }}_3$ and ${\mathrm{\Gamma }}_4$, respectively. Meanwhile, N_i nodes are distributed inside the computational domain $\mathrm{\Omega }$. (see Figure 1a). The nodes along ${\mathrm{\Gamma }}_1$ are satisfied with the initial conditions $U_0\left(x\right)$, and the nodes along ${\mathrm{\Gamma }}_2$ and ${\mathrm{\Gamma }}_3$ are satisfied with the boundary conditions $f_1\left(a,t\right)$ and $f_2\left(b,t\right)$, respectively. Therefore, the algebraic equations along the ${\mathrm{\Gamma }}_1$, ${\mathrm{\Gamma }}_2$ and ${\mathrm{\Gamma }}_3$ are described as,

$$\{\begin{array}{l}{F}_i^{}=u_i-U_0\left(x_i\right), i=1,2,\dots ,N_{b1},\\ F_i^{}=u_i-f_2\left(a,t_i\right), i=N_{b1}+1,N_{b1}+2,\dots ,N_{b1}+N_{b2},\\ F_i^{}=u_i-f_1\left(b,t_i\right), i=N_{b1}+N_{b2}+1,N_{b1}+N_{b2}+2,\dots ,N_{b1}+N_{b2}+N_{b3},\end{array}$$

(16)

The remaining domains, ${\mathrm{\Gamma }}_4$ and $\mathrm{\Omega }$, are both satisfied by the governing equation, Equation (1), and they can be discretized by Equation (15) as follows:

$$\begin{array}{l}{F}_i^{}=\left\{\begin{array}{l}\left(w{t}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}w{t}_{i,j}^{}{u}_{i,j}}\right)-\left(wxx{t}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}wxx{t}_{i,j}^{}{u}_{i,j}}\right)+\dots \\ \left(\alpha +\beta u_i^{\omega }\right)\left(w{x}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}w{x}_{i,j}^{}{u}_{i,j}}\right)-\dots \\ q\left(w{x}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}w{x}_{i,j}^{}{u}_{i,j}}\right)\left(wx{x}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}wx{x}_{i,j}^{}{u}_{i,j}}\right)-\dots \\ u_i\left(wxx{x}_{i,0}^{}{u}_i+{\displaystyle \sum _{j=1}^{n_s}wxx{x}_{i,j}^{}{u}_{i,j}}\right)\end{array}\right\}, \\ i=N_{b1}+N_{b2}+N_{b3}+1,N_{b1}+N_{b2}+N_{b3}+2,\dots ,N_B,N_B+1,N_B+2,\dots ,N_T,\end{array}$$

(17)

where $N_B$ is a number of nodes set up on the whole boundary, and defined as $N_B=N_{b1}+N_{b2}+N_{b3}+N_{b4}$; $N_T$ is the number of total nodes in the computational domain, and defined as $N_T=N_B+N_i$. In Equations (16) and (17), $N_T$ algebraic equations are acquired and denoted as $F={\left[F_1^{},F_2^{},\dots ,F_i^{}\right]}^{\mathrm{T}}, i=1,2,\dots ,N_T$. Meanwhile, $N_T$ unknown values are defined as $U={\left[u_1^{},u_2^{},\dots ,u_i^{}\right]}^{\mathrm{T}}, i=1,2,\dots ,N_T$. Then, the two-step Newton’s method is applied to solve this nonlinear algebraic system, and the specific numerical process is described in the next subsection.

3.2. Two-Step Newton’s Method

In this research, the famous nonlinear solver, Newton’s method, was applied for iterating the numerical solution of the proposed ST-GFDM. Newton’s method is a classic solver for nonlinear systems based on the first-order Taylor expansion formula. Newton’s method possesses characteristics that make it easy to program and ensure fast convergence. The iteration formula for Newton’s method is written as follows:

$$U^{k+1}=U^k-{\left(J^{-1}\right)}^k{F}^k,k=1,2,3,\dots ,$$

(18)

where J is the Jacobian matrix, and the elements are expressed as $J_{i,j}=\partial F_i/\partial U_j$; k is the number of iterations; $U^{k+1}$ and $U^k$ are the vectors of unknowns at the (k +1)-th and k-th iterations, respectively. Although Newton’s method is simple to apply, the inverse matrix of J is a troublesome issue during the simulation. To avoid the computation of the inverse matrix, the two-step iteration of Newton’s method is given [46,48,51] as follows:

$$J^k\mathrm{\Delta }{U}^k=-F^k,$$

(19)

$$U^{k+1}=U^k+\mathrm{\Delta }{U}^k,$$

(20)

where $\mathrm{\Delta }{U}^k$ is the numerical increment at the kth-step iteration.

The two-step process described in Equations (19) and (20) is employed as a means to avoid computing the inverse of the Jacobian matrix, and the $\mathrm{\Delta }{U}^k$ can be solved as a linear system. Thus, the nonlinear algebraic system obtained by the ST-GFDM can be solved efficiently. Furthermore, the iterative system will converge once the convergence condition $\delta =Max\left|\mathrm{\Delta }{U}^k\right|\le {10}^{-9}$ is reached and the numerical solutions of the proposed ST-GFDM are obtained.

3.3. Time-Marching Approach

In the previous subsections, the numerical process of the proposed ST-GFDM scheme has actually been introduced, but those processes only work in a single ST domain, $\mathrm{\Omega }=\left[a,b\right]\times \left[0,\mathrm{dt}\right]$. To deal with numerical problems with unpredictable end-times or long temporal scales, a simple technique, the time-marching approach, is used in the proposed numerical scheme. The schematics of the time-marching approach are demonstrated in Figure 2. In Figure 2a, a small-scale ST domain, $\mathrm{D}1=\left[a,b\right]\times \left[0,\mathrm{dt}\right]$, is formed when the simulation starts and the distribution of nodes is set up like in Figure 1a. In this ST domain, the initial conditions $U_0\left(x\right)$ are satisfied along ${\mathrm{\Gamma }}_1^1$, ${\mathrm{\Gamma }}_2^1$, and ${\mathrm{\Gamma }}_3^1$, which satisfy the boundary conditions $f_1\left(a,t\right)$ and $f_2\left(b,t\right)$, respectively. The governing equation (Equation (1)) is satisfied on ${\mathrm{\Gamma }}_4^1$ and $\mathrm{\Omega }$. Then, the numerical processes in Section 3.1 and Section 3.3 are applied to obtain the numerical solutions in the D1 domain. Once the numerical procedures within D1 are completed, the D1 will be shifted by distance dt along the t-axis, and the 2nd small-scale ST domain $\mathrm{D}2=\left[a,b\right]\times \left[\mathrm{dt},2\mathrm{dt}\right]$ is formed as shown in Figure 2b.

The numerical processes within D2 are the same as the D1, but the boundary conditions along ${\mathrm{\Gamma }}_1^2$ are different. In ${\mathrm{\Gamma }}_1^2$, the numerical solutions at ${\mathrm{\Gamma }}_4^1$ in D1 are treated as the initial conditions for D2. Thus, the numerical solutions in D2 can be solved by the proposed ST-GFDM and the two-step Newton’s method. It should be noted that the coordinates between D1 and D2 are relative locations. Thus, the weighting coefficients of the proposed ST-GFDM in D1 are the same as those of each ST domain generated by the time-marching approach. After moving the ST domain within n times, the n-th ST domain, $\mathrm{Dn}=\left[a,b\right]\times \left[\left(\mathrm{n}-1\right)\mathrm{dt},\mathrm{n}\times \mathrm{dt}\right]$, is generated (see Figure 2c), and the numerical solutions are solved. The simulation will finish when n × dt reaches the specified end-time determined by the numerical example. To summarize the above numerical process, the flowchart is shown in Figure 3. In the next section, four numerical examples are provided to verify the proposed ST-GFDM scheme.

4. Numerical Examples

In this section, the proposed ST-GFDM is validated by testing four numerical examples. These examples include the Fornberg–Whitham equation and its modified form, the modified Camassa–Holm equation, and the modified Degasperis–Procesi equation. The following symbols are defined: T is the end time of the numerical problem, N_D is the number of times that the ST domain has been moved by the time-marching method, and is written as $N_D=T/\mathrm{dt}$. The distribution of nodes for each example is displayed in Figure 1a. To analyze the accuracy of the proposed ST-GFDM, the numerical errors are defined in this study as follows:

$$\{\begin{array}{l}{L}_{\infty ,\mathrm{n}}={\underset{1\le i\le N_T}{\max }\left|u_i^{ext}-u_i^{num}\right||}^{\mathrm{n}},\\ L_{\infty }=\max \underset{1\le \mathrm{n}\le N_D}{\left\{L_{\infty ,\mathrm{n}}\right\},}\end{array}$$

(21)

$$\{\begin{array}{l}{\mathrm{RMSE}}_{\mathrm{n}}={\sqrt{{\displaystyle \sum _{i=1}^{N_T}{\left(u_i^{ext}-u_i^{num}\right)}^2}/N_T}|}^{\mathrm{n}},\\ \mathrm{RMSE}=\max \underset{1\le \mathrm{n}\le N_D}{\left\{{\mathrm{RMSE}}_{\mathrm{n}}\right\},}\end{array}$$

(22)

where superscript ext is the exact solution; num denotes numerical solutions. In Equation (21), $L_{\infty ,\mathrm{n}}^{}$ is the maximum absolute errors (MAEs) in Dn, and $L_{\infty }^{}$ is the maximum MAEs in the entire simulation. Meanwhile, the root-mean-squared error (RMSE) in each ST domain is defined as RMSE_n in Equation (22), and RMSE is the maximum RMSE in the whole simulation. In addition, the initial guesses are the initial conditions of each numerical problem for iteration in D1, and the numerical solutions in Dn − 1 are initial guesses of numerical processes in Dn.

4.1. Example 1

The Fornberg–Whitham equation [9,11] is solved in this subsection. This mathematical model is presented to describe the behavior of the breaking wave and has solitary or traveling wave solutions. The exact solution of the first numerical example is described as,

$$u^{ext}\left(x,t\right)=e^{\left(\frac{x}{2}-\frac{2t}{3}\right)},$$

(23)

and the initial condition and boundary conditions are represented as,

$$U_0\left(x\right)=e^{\frac{x}{2}}, f_1\left(t\right)=e^{\left(\frac{a}{2}-\frac{2t}{3}\right)}, f_2\left(t\right)=e^{\left(\frac{b}{2}-\frac{2t}{3}\right)}.$$

(24)

The following parameters are applied in this numerical example: a = −5, b = 5, dt = 0.5, T = 5, and N_D = 10. In Figure 4, the numerical solutions are obtained by applying N_T = 8421 and n_s = 23. The numerical solutions at different time levels are compared with the exact solutions and demonstrated in Figure 4a. This figure shows that the numerical solutions agree well with the exact solutions. Figure 4b is a 3D surface plot of numerical solutions to observe the numerical behavior of the Fornberg–Whitham equation through time.

Table 2. Numerical errors of Example 1 with n_s = 20 using different N_T at different t, convergence rates, and CPU-time.

	t = 1 (n = 2)		t = 3 (n = 6)		t = 5 (n = 10)		Convergence Rates (${\mathit{L}}_{\mathit{\infty },10}$ Are Used)	CPU-Time
NT	${\mathit{L}}_{\mathit{\infty },2}$	RMSE2	${\mathit{L}}_{\mathit{\infty },6}$	RMSE6	${\mathit{L}}_{\mathit{\infty },10}$	RMSE10
2211	1.725 × 10−4	2.424 × 10−5	2.251 × 10−5	4.081 × 10−6	3.250 × 10−5	5.891 × 10−6	-	0.7829 s
4215	7.653 × 10−5	8.978 × 10−6	1.206 × 10−5	1.780 × 10−6	1.703 × 10−5	2.693 × 10−6	1.918	1.5629 s
8421	3.267 × 10−5	3.071 × 10−6	6.677 × 10−6	8.855 × 10−7	9.153 × 10−6	1.214 × 10−6	1.742	3.5874 s
12,025	2.147 × 10−5	1.773 × 10−6	5.152 × 10−6	6.536 × 10−7	6.891 × 10−6	8.334 × 10−7	1.543	5.9342 s
18,631	1.313 × 10−5	9.162 × 10−7	3.856 × 10−6	4.747 × 10−7	4.927 × 10−6	5.949 × 10−7	1.528	10.6508 s

,

Table 3. MAEs and RMSEs of Example 1 with N_T = 8421 by using different n_s at different t.

	t = 1 (n = 2)		t = 3 (n = 6)		t = 5 (n = 10)
ns	${\mathit{L}}_{\mathit{\infty },2}$	RMSE2	${\mathit{L}}_{\mathit{\infty },6}$	RMSE6	${\mathit{L}}_{\mathit{\infty },10}$	RMSE10
20	3.267 × 10−5	3.071 × 10−6	6.677 × 10−6	8.855 × 10−7	9.153 × 10−6	1.214 × 10−6
23	5.158 × 10−5	5.737 × 10−6	1.961 × 10−5	2.011 × 10−6	1.301 × 10−5	1.563 × 10−6
26	5.093 × 10−5	5.713 × 10−6	3.067 × 10−5	3.064 × 10−6	1.264 × 10−5	1.566 × 10−6

and

Table 4. Comparisons of the numerical results of Example 1 with exact solutions at various locations at different time levels (N_T = 8421 and n_s = 23).

		t = 1 (n = 2)		t = 3 (n = 6)		t = 5 (n = 10)
		uext	ST-GFDM	uext	ST-GFDM	uext	ST-GFDM
x	−4	0.06948345	0.06948238	0.01831564	0.01831622	0.00482795	0.00483149
	−2	0.18887560	0.18887409	0.04978707	0.04978139	0.01312373	0.01312198
	0	0.51341712	0.51342701	0.13533528	0.13532764	0.03567399	0.03566356
	2	1.39561243	1.39565330	0.36787944	0.36788606	0.09697197	0.09696040
	4	3.79366789	3.79371117	1.00000000	1.00001868	0.26359714	0.26359633
$L_{\infty ,n}$		5.158 × 10−5		1.961 × 10−5		1.301 × 10−5

are provided to verify the accuracy and parameter sensitivity of the proposed method. The MAEs and RMSEs are calculated using N_T and n_s at different time levels. In Table 2, the results for different time levels are displayed, showing the MAEs and RMSEs obtained for various N_T. These error values demonstrate that increasing N_T leads to more accurate numerical solutions. Meanwhile, the convergence rate and CPU times are shown in Table 2. Since the proposed numerical scheme is a 3rd-order ST-GFDM combined with Newton’s method, according to the theoretical convergence rates of these two methods, the theoretical convergence rates of the proposed numerical scheme should be 2nd-order convergence. From Table 2, the numerical performance is between 1.5–2. On the other hand, the CPU times are provided in Table 2 as well. Table 2 reveals that simulations are completed within 4 s when N_T is below 9000, and the errors table indicates that satisfactory numerical solutions can already be achieved at this N_T. With an end time of 5 for the first example, the proposed numerical scheme can successfully complete the simulation within the problem’s specified end time using approximately 9000 nodes. Despite the increase in CPU time as N_T increases, the CPU time remains below 10 s, demonstrating the efficiency of the proposed ST-GFDM in solving problems with long-time scales.

In Table 3, MAEs and RMSEs are also presented using different n_s. The errors indicate that the parameter n_s exhibits low sensitivity. Additionally, both tables demonstrate no accumulation of errors throughout the numerical simulation. This observation supports the effectiveness of the ST approach in mitigating the accumulation of numerical errors, a concern that hybrid numerical schemes often encounter when applied to time-dependent problems.

In order to further validate the capability of the proposed numerical scheme, the numerical solutions located at specific positions (x,t) are presented in Table 4 and compared with the exact solutions. In this table, the function values are taken to eight decimal places to show the precision accuracy of the proposed ST-GFDM. Additionally, the proposed ST-GFDM is compared with the HAM [10], ADM [7], RKHSM [7], MVIA-I [5], and the 2nd-order finite difference method combined with the fourth-order Runge–Kutta method (FDM-RK4). Those comparisons are shown in

Table 5. Comparisons of the numerical results of Example 1 with the exact solution and the other approximation methods at various locations at t = 5 (n = 10, N_T = 8421, and n_s = 23).

x	uext	HAM [10]	ADM [7]	RKHSM [7]	MVIA-I [5]	FDM-RK4	ST-GFDM	MAEs
−4	0.00482795	0.00487526	0.00317192	0.00482795	0.00482795	0.0048810	0.00482958	1.626 × 10−6
−2	0.01312373	0.01325233	0.00862217	0.01312373	0.01312373	0.0127930	0.01312294	7.893 × 10−7
0	0.03567399	0.03602358	0.02343750	0.03567400	0.03567399	0.0349080	0.03566866	5.337 × 10−6
2	0.09697197	0.09792224	0.06370972	0.09697195	0.09697197	0.1008894	0.09696298	8.992 × 10−6
4	0.26359714	0.26618027	0.17318095	0.26359729	0.26359714	0.3172236	0.26359021	6.925 × 10−6

. In these comparative data, it can be seen that the proposed ST-GFDM performs better than the HAM and ADM. Although the accuracy is weaker than the RKHSM and MVIA-I, solutions with numerical errors below 10⁻⁶ or even 10⁻⁷ can be obtained, which is already an acceptable numerical result for engineering applications.

To compare with the original numerical scheme for time-dependent problems, the FDM-RK4 is applied to obtain the numerical results within ∆x= 0.5 and ∆t = 0.0001. The numerical solutions of the FDM-RK4 show that the results have a maximum absolute error of the order of 10⁻¹. In our test, the smaller ∆x needs a smaller ∆t for FDM-RK4, and the simulation will diverge during the simulation because the end time could not reach t = 5. In Figure 5, it can be clearly seen that the numerical error of ST-GFDM is maintained at the same magnitude during the simulation, and the numerical results are more accurate than those of FDM-RK4. This numerical performance shows that using the space–time approach, the numerical errors will not accumulate during simulation for each ST domain. The above discussions show that the proposed ST-GFDM can accurately and stably solve the Fornberg–Whitham equation. Since the capability of the proposed ST-GFDM has been verified in this numerical example, the following numerical tests are solved using specific parameters, and the obtained numerical results are compared with other approximation algorithms.

4.2. Example 2

The second example is a numerical test for the modified Fornberg–Whitham equation [3,5,6,7,9]. The exact solution is given as follows:

$$u^{ext}\left(x,t\right)=\frac{3}{4}\left(\sqrt{15}-5\right){\mathrm{sech}}^2\left(c\left(x-\left(5-\sqrt{15}\right)t\right)\right),$$

(25)

where $c=\frac{1}{20}\sqrt{10\left(5-\sqrt{15}\right)}$. The initial condition and boundary conditions are denoted as follows:

$$\{\begin{array}{l}{U}_0\left(x\right)=\frac{3}{4}\left(\sqrt{15}-5\right){\mathrm{sech}}^2\left(cx\right), \\ f_1\left(t\right)=\frac{3}{4}\left(\sqrt{15}-5\right){\mathrm{sech}}^2\left(c\left(a-\left(5-\sqrt{15}\right)t\right)\right), \\ f_2\left(t\right)=\frac{3}{4}\left(\sqrt{15}-5\right){\mathrm{sech}}^2\left(c\left(b-\left(5-\sqrt{15}\right)t\right)\right).\end{array}$$

(26)

The following parameters are used in this numerical example: a = −10, b = 10, dt = 0.2, T = 1, N_D = 10, N_T = 40,421, and n_s = 20. The numerical solutions are displayed in Figure 6. In Figure 6a, the comparison of the exact solutions and the numerical solutions at specific time levels (t = 0.2, 0.4, 0.8, and 1) is shown. These numerical results are in good agreement with the exact solutions. In addition, the MAEs at specific points (x,t) are listed in

Table 6. MAEs of Example 2 with the exact solution and the other numerical methods.

	x = 2.5			x = 5
t	MVIA [5]	VIM [8]	ST-GFDM	MVIA [5]	VIM [8]	ST-GFDM
0.02	3.780 × 10−5	1.180 × 10−4	1.883 × 10−7	9.966 × 10−6	2.124 × 10−5	4.272 × 10−8
0.04	7.240 × 10−5	2.363 × 10−4	1.949 × 10−7	1.778 × 10−5	4.797 × 10−5	4.077 × 10−8
0.06	1.036 × 10−4	3.547 × 10−4	1.996 × 10−7	2.333 × 10−5	8.029 × 10−5	3.867 × 10−8
0.08	1.313 × 10−4	4.731 × 10−4	2.042 × 10−7	2.653 × 10−5	1.183 × 10−4	3.647 × 10−8
0.1	1.552 × 10−4	5.914 × 10−4	2.087 × 10−7	2.727 × 10−5	1.622 × 10−4	3.419 × 10−8

with the errors of the MVIA [5] and VIM [8]. Obviously, the error performance of the proposed numerical scheme is better than that of the other two approximation algorithms. Also, it verifies that the proposed ST-GFDM is fully capable of handling such complex PDEs.

4.3. Example 3

The 3rd example is the modified Camassa–Holm equations [1,2]. The exact solution is given as follows:

$$u^{ext}\left(x,t\right)=-2{\mathrm{sech}}^2\left(\frac{x}{2}-t\right),$$

(27)

The following are applied in this numerical example: a = −10, b = 10, dt = 0.5, T = 2, N_D = 4, N_T = 37,231, and n_s = 20. The initial condition and boundary conditions are obtained by introducing the above parameters into Equation (27). The numerical solutions are illustrated in Figure 7. In these plots, the obtained numerical results are well compared with the exact solution in Figure 7a. Moreover, Figure 7b demonstrates the physical behavior of the wave transmission. Additionally,

Table 7. Comparisons of the numerical results of Example 3 with the exact solution and the other numerical methods at various locations at t = 0.05 and t = 0.1.

x	t	uext	ST-GFDM	VIM [2]	MAEs
6	0.05	−0.021795977	−0.021798486	−0.019790189	2.508360 × 10−6
8	0.05	−0.002963750	−0.002960395	−0.002682964	3.355673 × 10−6
6	0.1	−0.024074444	−0.024070611	−0.019848303	3.832791 × 10−6
8	0.1	−0.003275195	−0.003265414	−0.002684058	9.781412 × 10−6

presents the numerical results at various locations, comparing them with the numerical results of the VIM [2]. The table shows that the proposed ST-GFDM successfully solves the Camassa–Holm equation and provides accurate numerical solutions.

4.4. Example 4

The final example is the modified Degasperis−Procesi equation [2,3]. The exact solution is given as follows:

$$u^{ext}\left(x,t\right)=-\frac{15}{8}{\mathrm{sech}}^2\left(\frac{x}{2}-\frac{5t}{4}\right),$$

(28)

The boundary condition and initial condition can be obtained from Equation (28). The following parameters are used: a = -10, b = 10, dt = 0.5, T = 2, N_D = 4, N_T = 37,231, and n_s = 20. The numerical solutions and data are given in Figure 8 and

Table 8. Comparisons of the numerical results of Example 4 with the exact solution and the other numerical methods at various locations at t = 0.05 and t = 0.1.

x	t	uext	ST-GFDM	VIM [2]	MAEs
6	0.05	−0.020951831	−0.020948113	−0.018566923	3.717910 × 10−6
8	0.05	−0.002851900	−0.002848801	−0.002682964	3.099397 × 10−6
6	0.1	−0.023714184	−0.023719627	−0.018635025	5.443036 × 10−6
8	0.1	−0.003229560	−0.003227788	−0.002516809	1.771923 × 10−6

, respectively. Due to the similarity between the modified Degasperis−Procesi equation and the modified Camassa–Holm equation, the numerical solutions’ behavior is also similar. These comparisons show that the proposed ST-GFDM consistently yields accurate numerical results, as demonstrated in the previous numerical examples. As described earlier, the proposed ST-GFDM effectively solves the shallow water wave problem and demonstrates stability in simulating problems over long time scales.

5. Conclusions

This study uses the proposed ST-GFDM to numerically solve nonlinear shallow water wave equations, specifically the modified Camassa–Holm equation, the modified Degasperis–Procesi equation, the Fornberg–Whitham equation, and its modified form. The ST-GFDM discretizes the governing equations in the ST domain, eliminating the need for mesh generation and simplifying the numerical procedures. The resulting nonlinear algebraic system is solved using the two-step Newton’s method. Adopting ST-GFDM makes the computational process more straightforward and avoids the complexities associated with mesh generation. One significant innovation of the ST-GFDM is treating the temporal derivative as a spatial derivative. The ST-GFDM enables the discretization of all partial derivatives and the efficient handling of mixed derivatives. Due to the Taylor series expansion and the moving least-squares method, the proposed meshless numerical scheme forms a nonlinear algebraic system, and this system has a sparse Jacobian matrix. Then, the two-step Newton–Raphson method is applied to efficiently solve this nonlinear algebraic system. To deal with long-time-scale numerical problems, the time-marching method is applied for moving the ST domain and saving computational resource, thereby addressing the shortcomings of the ST approach.

Four numerical examples are presented to verify the proposed ST-GFDM scheme. In the first numerical case, by varying different parameters, it is observed that increasing the N_T leads to more accurate numerical solutions. Additionally, adjusting the parameter n_s in a supporting domain demonstrates that n_s has low sensitivity. The numerical solutions have good agreements with the exact solutions and the solutions of other approximation methods in four examples. Furthermore, the accumulation of numerical errors during the time progression is insignificant. The proposed ST-GFDM scheme offers several advantages, including accurate treatment of transient problems, computational efficiency, avoidance of mesh generation, and flexibility in handling complex PDEs. This observation highlights the scheme’s capability to effectively solve transient problems while maintaining numerical accuracy. Future research can explore the application of the ST-GFDM scheme in computational fluid dynamics and tackle more challenging problem domains.