Exact and numerical approaches for solitary and periodic waves in a (2+1)-dimensional breaking soliton system with adaptive moving mesh

1.   Introduction

Research on nonlinear evolution equations (NLEEs) plays a significant role in establishing many fields of science and technology. These studies encompass areas such as nonlinear dynamics ^[1,2], fiber optics ^[3,4], plasma physics ^[5], high-frequency waves ^[6,7], and ocean engineering ^[8]. Solitons, a unique wave phenomenon, preserve their structure and energy over extensive distances, which is crucial for comprehending the dynamics of intricate nonlinear systems. Solitons are pivotal in analyzing nonlinear systems due to their stability and capacity to encapsulate fundamental insights into system behavior ^[9]. The identification and examination of solitonic solutions not only deepen our theoretical comprehension but also promote practical advancements in innovations derived from nonlinear phenomena ^[10]. Due to the critical significance of solitonic solutions in comprehending advanced nonlinear frameworks across many applications, substantial efforts have been dedicated to devising reliable methods for their compilation and validation. Over recent decades, researchers have devised numerous techniques to analyze and formulate solitonic solutions for NLEEs. The methods encompass the Hirota bilinear method ^[11,12], direct algebraic technique ^[13], hyperbolic sine and hyperbolic cosine technique ^[14], modification and extension of classical Bäcklund transformation ^[15,16], Darboux transformation ^[17], breather limit method ^[18], basic (G'/G)-expansion method ^[19], generalized exponential rational function method ^[20,21], exponential rational function method, improved modified extended hyperbolic tangent function method ^[22], among others.

Our main objective of this work is to obtain the (2+1)-dimensional generalized breaking soliton system (GBSS). This system is a reliable theoretical framework for simulating the dispersion, nonlinearity, and cross-wave coupling of media such as shallow water waves, optical fiber systems, and magnetized plasma. Researchers have focused on specific instances of GBSS, including ^[23,24,25], whereas we seek to formulate the equation in a more comprehensive and generalized manner. We aim to reveal dynamic behaviors and solutions that may have been overlooked in prior case-specific research by expanding its scope. The generalized form provides a cohesive framework that facilitates the comprehension of various nonlinear phenomena encompassed by the GBSS, hence enhancing theoretical understanding. We enhance previous methodologies and contribute to the expanding literature on more compact and efficient methods for analyzing nonlinear evolution equations. This work is fundamentally based on the GBSS, which is articulated by the following equations:

The parameters $\alpha, \beta, \gamma, \delta$, and $\lambda$ are real, whereas the functions $W = W(x, y, t)$ and $V = V(x, y, t)$ represent the values provided across the two structures, respectively. $W_t$ dictates the temporal evolution of the wave. The third derivative with respect to $x$ yields $\alpha W_{x x x}$, representing $x$-axial dispersion. Dispersion is a mechanism that distributes the wave, so mitigating the steepening caused by nonlinear phenomena. The mixed derivative $\beta W_{x x y}$ characterizes the cross-dispersion in the spatial domains. Consequently, it illustrates the intricate interplay of waves in both directions. The nonlinear term $\gamma W W_x$ is the nonlinear self-steepening. The wave profile steepens as its amplitude increases. This elucidates why waves of greater amplitude propagate more rapidly and the wave profile steepens. Furthermore, the variables $\delta W V_x$ and $\lambda W_x V$ represent bidirectional interactions between the two wave fields, incorporating the motion of the waves along both axes. It is clarified how forward-traveling waves affect the motion of backward-traveling waves. The second equation of the system, $W_y = V_x$, defines the relationship between spatial gradients in $x$ and $y$, which reduces the degrees of freedom.

In the context of fluid dynamics, GBSS concerns the interaction of surface waves on the water surface and internal waves lying beneath the surface, which are induced by mechanical forces in the various density layers. While surface waves caused by the wind and gravitational effects have a short period and travel fast, internal surface waves have longer wavelengths and move at much slower speeds. The GBSS records these interactions by equilibrating nonlinear effects from wave steepening (terms at the $\gamma W W_x$ level) with dispersive effects that induce wave dispersion (the $\alpha W_{x x x}$ term). The cross-interaction terms $\delta W V_x$ and $\lambda W_x V$ represent the energy transfer between surface and interior waves during contact. The situation may involve internal waves modulated by surface waves leading to soliton-type nonlinear interactions that change the shape of the wave, possibly inducing phenomena such as soliton fission ^[26]. In nonlinear optics, the (2+1) dimensional GBSS described above is motivated by and gives the interaction of light pulses, or more specifically, solitons propagating through optical fibers, where the nonlinear refractive index and its dispersion are adjusted to allow the generation of solitons. As these light pulses traverse the fiber, they are affected by a nonlinear medium, resulting in pulse steepening and dispersion, which causes pulse broadening. The GBSS incorporates these dynamics via terms such as $\alpha W_{x x x}$ for dispersion and $\gamma W W_x$, which denotes nonlinear effects. The cross-interaction terms $\delta W V_x$ and $\lambda W_x V$ denote the coupling between two light pulses, which may propagate in distinct polarisation modes or adjacent fibers, facilitating energy transfer and intricate interactions. One such interaction is the steepening or fragmentation into sub-pulses (the breaking soliton phenomenon) of a single pulse. This modeling is essential for optical communications: A comprehensive comprehension of soliton activity and interaction can enhance data transmission quality and extend transmission distances ^[27]. Beyond the interest from their framework of the GBSS, there also exists a practical interest in a number of different fields. Consider, for example, the application of GBSS in hydrodynamics; it can model the movement of surface and internal waves in oceans and lakes, which is critical for understanding energy transfer and wave interactions in these domains. In plasma physics, the system can capture the behavior of plasma waves in a controlled fusion device, where the interaction of solitons plays an important role in sustaining the stability of the plasma.

The generalized direct algebraic method is effective for partial nonlinear evolution and wave equations. This method derives solutions by positing trial functions as exponential or rational functions and determining their coefficients through balancing principles. It extends the type of solitary wave solutions with analytic functions of the real exponential solutions to the linearized equation ^[13]. The method has been used for various possible equations, adding KdV, Boussinesque, and Schrödinger equations ^[28,29]. The phase space analysis ^[30] has been carried out to elucidate and, to some extent, prove the method. There have been new extensions of the technique, like the new extended direct algebraic method in which ^[31] applied this technique to Tzitzica-type evolution equations and other nonlinear models. The construction of exact solutions for various nonlinear equations confirms the method's effectiveness and accuracy ^[29]. Other relevant ones are the Newton polygon method and its extensions, which are proposed to solve generalized algebraic equations ^[32] as well as general ones for solving linear algebraic problems ^[33].

Analytical methods, including the generalized direct algebraic method, yield exact solutions; however, numerical methods have also played an essential role in obtaining approximate solutions for nonlinear PDEs for which there are no or very complex closed forms. Among the numerous numerical methods, adaptive moving mesh methods have been developed as a robust approach for problems with sharp gradients, singularities, and multi-scale features. Using these methods, the computational mesh is dynamically adjusted to regions of interest to improve accuracy and efficiency and reduce the cost of computation. The adaptive moving mesh methods designed and investigated in recent papers show the ability to solve high-order parabolic PDEs and capture intricate solution features. In addition, other works have attended to the same issue using finite difference schemes and showed unconditional stability over a wide range of parameters, thereby increasing the trust of the numerical solutions ^[34]. It should be emphasized that spectral methods provide an alternative possibility, which is most of the time very precise, in solving high-order nonlinear evolution equations. Recent attempts covering wave phenomena and dispersion have proven the efficiency of spectral collocation and similar approaches^[35,36]. However, these spectral methods are complemented by our adaptive moving mesh framework, which actively moves grid points to capture sharp features, minimizing the chance of complex features being under-resolved. This complementarity ensures that our approach can be used in conjunction or in contrast with spectral methods for a wider scope of challenging wave problems.

In recent decades, Moving Mesh Partial Differential Equations (MMPDEs) have been raised as one of the economical methods for solving time-dependent PDEs. A significant foundation for adaptive moving mesh techniques was set by Huang and Russell ^[37] when they introduced the MMPDEs. The paradigm of MMPDE is based on a monitor function that can quantify specific features of the solution, such as curvature, arc length, or error estimates for adaptive mesh point redistribution. Concentrating mesh points in critical regions enhances computational efficiency and solution accuracy. The equidistribution principle, which is at the core of this approach, ensures that the integral of the monitor function is constant over each sub-domain, thus balancing the resolution to make the domain's borders less problematic. Budd and other contributors supplemented these guidelines by dealing with difficulties such as mesh tangling and blow-up problems that may arise when solving higher-order PDEs. Moreover, the exact and numerical solutions of the system (1.1) are discussed. Section 2 is devoted to finding the exact solutions of the system using the generalized direct algebraic method. In Section 3, we deal with obtaining numerical solutions using the adaptive moving mesh method based on the finite difference methods.

2.   Exact solutions of GBSS

Considering the development equation with physical fields $W(x, y, t)$ and $V(x, y, t)$ in the variables $x,$ $y,$ and $t,$ we express it in the following form:

To investigate travelling-wave solutions of the system (1.1), we assume:

where $k$ and $c$ are arbitrary constants, while $w$ denotes the wave speed. Under this assumption, the non-linear evolution (2.1) reduces to an ordinary differential equation (ODE):

Here, $P_2$ is a polynomial in $q(\zeta), \; p(\zeta),$ and the derivatives of $p(\zeta)$ and $q(\zeta)$ concerning $\zeta$. In this step, we utilize the traveling wave solution to merge the original PDE system into an ODE and simplify the mathematics. This method substitutes the partial derivatives in $x$ and $y$ through a single derivative with respect to $\zeta$, which accounts for the propagation of the wave along a spatially chosen combination of directions. Such change is often helpful in finding solitonic or periodic wave solutions that correspond to some known analytical forms or can be used for further physical interpretation ^[38].

2.1.   Generalized direct algebraic method and ODE derivation

According to the generalized direct algebraic method, the solutions in relation to the function $H(\zeta)$, as solutions of (2.3), are sought in the form:

where $\mu_k$ and $\sigma_k$ are constants to be determined, and $N$ is an integer fixed by balancing the highest degree of the nonlinear terms against the highest order of the derivatives. The function $H(\zeta)$ satisfies

where $\varepsilon$ is a user-specified constant, usually $\pm1$, and the parameters $\tau_i$ for $i = 0, 1, 2, 3, 4$ are provided in Table 1 ^[39]. The table presents particular relationships between the values of $\tau_i$ that lead to elementary functions, such as hyperbolic sine and cosine, trigonometric functions, or rational forms depending on the sign and magnitude of the coefficients.

At this point, we apply the above method to derive analytical solutions to the GBSS system (1.1). The traveling-wave substitutions (2.2) simplify the PDEs into a set of ordinary differential equations:

By integrating the second equation, Eq (2.6), yields

Substituting (2.7) back into the first equation, Eq (2.6), and integrating once with respect to $\zeta$ leads to:

Based on the generalized algebraic technique, the solutions of (2.8) are expressed as follows:

such that $H(\zeta)$ satisfies

2.2.   Analytical solutions: Case studies

For each of the scenarios outlined in Table 1, we employed symbolic computation (e.g., the Wolfram Language) to determine the specific constants $\mu_0, \mu_1, \mu_2, \sigma_{1}, \sigma_{2}, c, k$, and others. As a result, we obtain the following families of solutions:

First case. $\tau_0 = 0$, $\tau_1 = 0$, $\tau_2 > 0$, $\tau_3 = 0$, $\tau_4 < 0$,

Second case. $\tau_0 = \tfrac{\tau_2^2}{4 \tau_4}$, $\tau_1 = 0$, $\tau_2 < 0$, $\tau_3 = 0$, $\tau_4 > 0$,

Third case. $\tau_0 = 0$, $\tau_1 = 0$, $\tau_2 < 0$, $\tau_3 = 0$, $\tau_4 > 0$,

Fourth case. $\tau_0 = \frac{\tau_2^2}{4 \tau_4}$, $\tau_1 = 0$, $\tau_2 > 0$, $\tau_3 = 0$, $\tau_4 > 0$,

Fifth case. $\tau_0 = 0$, $\tau_1 = 0$, $\tau_2 > 0$, $\tau_3 \neq 0$, $\tau_4 = 0$,

The resulting solutions show a variety of waveforms that some are localized solitons (notably involving $\operatorname{sech}^2$, $\tanh^2$, or related hyperbolic functions), while others manifest as periodic waves (often involving $\sec^2$ or $\csc^2$). The choice of signs for $\tau_2$ and $\tau_4$ plays a key role in shaping the wave profile. For instance, positive $\tau_2$ together with negative $\tau_4$ can give rise to localized solitonic structures that decay rapidly, whereas certain trigonometric ratios (like $\tan$ or $\sec$) dominate when $\tau_2$ and $\tau_4$ are both positive, leading to repeating, periodic solutions. By altering these parameter constraints, one effectively tunes the system from sharply localized waveforms toward smoother, periodic oscillations, revealing the flexibility of the model in describing diverse physical phenomena.

2.3.   Graphical analysis and discussion

In this subsection, we summarize the major findings from the graphical analysis while mentioning the analytical solutions achieved with $\tau_2 = 2.5$ (or $-2.5$), $\alpha = -0.2$, $\beta = 3.1$, $\gamma = -0.5$, $\delta = 3.7$, $\lambda = 0.5$, $k = 0.4$, and $c = 0.2$ which capture the subtle interplay of nonlinearity and dispersion. In the last case with the refined values of $\sigma_1 = 0.5$ and $\tau_3 = 0.8$, the parametrizations visually show the so-called sensitivity of solution behavior to these parameters, which are critical aspects of the model's stability as well as its dynamic responses.

Figures resulting from the generalized direct algebraic method showcase different solution families to the nonlinear system, each of which is related to a distinct value of the parameter set. In these figures, parameters $k$ and $c$ represent the wave numbers or velocities in the $x$ and $y$ axes, and the system's dispersive as well as nonlinear features are controlled by $\alpha, \beta, \gamma, \delta$, and $\lambda$ while the auxiliary parameters $\tau_0, \tau_1, \tau_2, \tau_3$, and $\tau_4$ become apparent in the course of the solution, determining the form of the solution a priori.

In Figure 1, the solutions represented as hyperbolic secant functions are obtained when $\tau_0 = 0$, $\tau_1 = 0, \tau_2 >0, \tau_3 = 0, \tau_4 < 0$ is satisfied. The $\operatorname{sech}$ function represents a hyperbolic secant, which has bright solitons as it decays very quickly from its amplitude maxima. In addition, this corresponds to the waveforms of $W_1$ and $W_2$ (and associated $V$ components) being localized non-dispersive pulses that travel with nearly constant shape, which implies a balance of nonlinearity and dispersion in the system.

Figure 2 corresponds to the region $\tau_0 = {\tau_2^2}/{4\tau_4}$, $\tau_1 = 0$, $\tau_2 < 0$, $\tau_3 = 0$, and $\tau_4>0$. Here, the solutions contain the functions hyperbolic tangent and hyperbolic cosecant. The $\tanh$ function gives rise to kink-type structures that join various asymptotic states, and the $\operatorname{csch}$ function generates pulses that decay very rapidly. The plots for $W_{4, 5, 6, 7}$, (and their corresponding $V$ -components) exhibit steep wavefronts and sharp transitions, which is the characteristic of these hyperbolic functions and the strong localized nature of the solutions.

In Figure 3, the parameters are set to $\tau_0 = 0$, $\tau_1 = 0$, $\tau_2 < 0$, $\tau_3 = 0$, $\tau_4>0$, which gives rise to solutions with a secant profile. As for the previous cases, the secant function is periodic and has vertical asymptotes. This creates repeating minima and maxima in the solutions $W_{10}$ and $W_{11}$ (with their respective $V$ -components), which appear mathematically as a periodic function. Physically, while this suggests a coherent wave structure, the singularities need, however, to be treated with caution in real-word applications.

The parameter set $\tau_0 = {\tau_2^2}/{4\tau_4}$, $\tau_1 = 0$, $\tau_2>0$, $\tau_3 = 0$, and $\tau_4>0$ is considered in Figure 4. The solutions are given as the tangent and cosecant functions. As these functions are both periodic and have singular behavior for certain arguments, the plots corresponding to $W_{13}$ through $W_{17}$ (and also the $V$ -components) appear to oscillate repeatedly with steep changes, which, depending on the physical system, can be viewed as resonant or unstable wave modes.

Figure 5, on the other hand, shows the case where $\tau_0 = 0$, $\tau_1 = 0$, $\tau_2 > 0$, $\tau_3\neq 0$, and $\tau_4 = 0$. In this case, the solutions have the form of a shifted $\operatorname{sech}^2$ form. This describes a bright soliton that has a smooth bell-shaped curve that decays rapidly away from the peak. The plots for $W_{19}$ through $W_{22}$ (and their corresponding $V$ -components) exhibit the typical localized behavior of such solitons with the additional parameter $\tau_3$ controlling the pulse's height and width. Collectively, these figures depict the wave structures that the system is capable of supporting. It is evident that with proper modifications of the parameters, one can attain transitioned solitary waves that are either periodic or even singular, which is a clear manifestation of nonlinearity and dispersion interplay in the model.

3.   Numerical solution

Moving mesh methods within r-adaptive spatial discretizations are useful for time-dependent problems through the use of coordinate transforms. Such techniques utilize mesh density functions in order to modify the mesh according to different features in the solution. This time, we utilized the two-dimensional MMPDEs on the specified physical ($\Omega_p$) and computational ($\Omega_c$) domains, over which we were able to extend our previous works, thus enabling the mesh to adapt based on a specified mesh density function. Before discussing the formulation of the equations, it is worth noting one of the benefits of the moving mesh approach is the capability of concentrating mesh points in areas of steep gradients or variation in the solution ^[40]. This significantly reduces computational workload when the solution dynamics are adjustable and ensures accurate capture of solution dynamic changes. Note that the method's effectiveness relies on specifying adequate mesh monitor functions and robust mesh-movement PDEs. In the core equations, MMPDE5 is commonly expressed in two-dimensional form as follows:

where $\boldsymbol{\eta} = (\eta, \xi)$ is the computational coordinates, $\boldsymbol{x} = (x, y)$ is the mapping of the physical domain, and $\Psi(\boldsymbol{x}, t)$ is the monitor mesh function, which depends on the solutions' variations. $\nu$ is the scalar parameter that scales the temporal speed responsive mesh moves. The right-hand side of Eq (3.1) is given by

Thus, with such an addition, we obtain two PDEs: One for $x(\eta, \xi, t)$ and another for $y(\eta, \xi, t)$. The governing structure of the equations should be the same, while the dependent variable distinguishes the equations.

In this case, $\Psi(\boldsymbol{x}, t)$ is dependent on a scalar $\nu$ that modifies how closely the mesh tracks the changes in the monitor function. MMPDEs are effective for adaptive mesh generation; however, they carry extra computational costs and may result in challenges like mesh tangling if not carefully handled. As addressing the adaptive mesh through PDEs in addition to the PDEs in question becomes extraneous, computational costs escalate. Even so, coordinate transformations that are MMPDE-based for 2D problems have their own limitations in this regard, which may cause mesh tangling and loss of mesh regularity. Among such methods, optimal transportation theory, for instance, those based on Monge-Ampère and Parabolic Monge-Ampère (PMA) equations, is considered to have more advantages in that they require fewer equations or lead to meshes that are more robust and regular ^[41]. Judging from previous remarks, the selection of a suitable monitor function $\Psi$ to be used to mesh is of great significance. They may also be selected for the set of error estimates: For example, polynomial interpolation or truncation error, or on the set of the solution features of the PDEs, such as arc length or curvature:

where $n = 2$ or $4$ and $\alpha$ and $\beta$ are adaptivity parameters that can change in space to accommodate sudden changes in the characteristics of the solution. In the process of designing and selecting monitoring functions, accuracy must be balanced against the extra effort required to measure it. Adaptive or highly oscillatory solutions may require more refined and specific monitor functions. In contrast, uncomplicated functions better serve more straightforward solutions. Thus, adjusting the monitoring function becomes a standard approach in MMPDE integration to obtain a smoother computational domain and reduce MMPDE embedding ^[40]. This confirms the theoretical and practical progress that are mutually integrated into the efficient use of adaptive moving mesh methods for simulation purposes.

3.1.   Semi-discretization scheme on adaptive grids

To transition from the computational domain $\Omega_c,$ represented by the coordinate $\boldsymbol{\eta} = (\eta, \xi)$, to the physical domain $\Omega_p$ defined by the coordinate $\boldsymbol{x}$, the adaptive moving mesh method employs a coordinate transformation.

where $\Omega_c$ is given as:

and $\Omega_p$ associated with the solution $W(x(\boldsymbol{\eta}, t), t)$ and $V(x(\boldsymbol{\eta}, t), t)$ as follows:

where the endpoints of $\Omega_p$ are:

In this context, $\Delta \eta = \frac{b-a}{L_x}$ and $\Delta \xi = \frac{d-c}{L_y}$ represent the uniform step sizes within $\Omega_c$, where $L_x$ and $L_y$ are positive integers indicating the number of mesh points in the $x$ and $y$ directions, respectively. The next step is to approximate the spatial derivatives using standard second-order centered differences. The term $\frac{\partial}{\partial \eta}\Bigl(\Psi\, \frac{\partial x}{\partial \eta}\Bigr)$ is handled by first calculating $\frac{\partial x}{\partial \eta}$ at $\Psi$ midpoints between grid nodes and then multiplying it by an appropriately averaged $\Psi$. A typical discretization in $\eta$ for the right hand side of Eq (3.1) is provided as

and similarly in $\xi$ direction. In practice, to have better stability, $\Psi_{j+1 / 2, k}$ is defined at midpoints by an average as follows:

Thus, the semi-discretization in $\boldsymbol{\eta}$ with continuing temporal time is provided as

where FD stands for finite difference. In order to account for the boundary conditions, one usually locks $\boldsymbol{x}$ at the domain edges if the physical boundary is not moving, thus preserving the outermost grid mesh at the required boundary positions.

As a result, inserting these finite-difference approximations into the MMPDE5 equations gives a set of ordinary differential equations (ODEs) in time concerning $x_{j, k}(t)$ and $y_{j, k}(t)$. These ODEs can be solved using classical integration methods for ordinary differential equations such as forward Euler, Runge-Kutta, and implicit schemes. The choice of a particular method relies upon the requirements for stability and accuracy of the results. The outcomes provide a discrete evolution of the mesh coordinates in time, which is adjusted according to the directed solution variation, controlled by the monitored function $\Psi$.

To implement the use of adaptive meshes for the discretization of the undertaken PDE, it is necessary to express the derivatives in the new coordinate system. By relying on the reality of $W_y = V_x$, mentioned above, system (1.1) can be rewritten using the Lagrangian form as follows:

where

and $J = x_{\eta} y_\xi - x_\xi y_{\eta}$.

In practice, the calculation of every partial derivative in combination with new coordinates requires careful handling of the terms such as $x_\eta, x_\xi, y_\eta,$ and $y_\xi$. Besides, the transformations based on Jacobians guarantee that the mesh can be smoothly adjusted to capture critical solution details, enabling localized refinement in important areas. Even though these transformations increase the difficulty level, the accuracy and speed that can be achieved make them worth the effort. Each of the terms in Eqs (3.5) and (3.6) is discretized using a computational finite-difference scheme in $\boldsymbol{\eta}$ as:

where

It is possible to mention here that $y_\xi = J \eta_x, \, y_{\eta} = -J \xi_x, \, x_\xi = -J \eta_y, \, x_{\eta} = J \xi_y.$ These connections are meaningful because they ensure that the physical terms capture the mesh and are consistent with the undertaken system. If these transformations are precisely done, higher-order derivatives with respect to $\boldsymbol{x}$ could be computed efficiently while still using the simpler $\boldsymbol{\eta}$ coordinate system.

Substituting (3.6)–(3.8) into (3.5), yield

In Figure 6, the evolution of the numerical solution of the (2+1)-dimensional GBSS problem, along with the adaptive temporal evolution of the mesh, is shown. In Figure 6(a), the solution $W_{22}(x, \, y = 2, \, t)$ for fixed $y = 2$ is demonstrated for $t = 0$ to $t = 10s$. The exact (analytical) solution is shown with a dashed black line, while the numerical solution obtained from the non-uniform mesh is shown in solid light blue (see the insets within in Figure 6(a)). Their proximity is consistent with the second-order accuracy discussed in the error analysis. In addition, Figure 6(b) shows the non-uniform mesh for $x \in [-20, \, 20]$ with the regions of steep gradients well captured. The figures demonstrate that adaptive methods give a good prediction of wave phenomena, aligning well with the principles and methods discussed in the following sections.

3.2.   Stability of the discrete MMPDE5

To see if the discrete version of MMPDE5 is stable, we can use a Von Neumann (Fourier) stability analysis. For simplicity, assume the monitor function $\Psi$ is a constant. Then, each component (say $x$) of the mesh mapping follows an equation like a 2D diffusion process,

where, $\eta$ and $\xi$ are the computational coordinates. We reformulate it as

Subsequently, we discretize spatially and employ an explicit time-stepping technique. This results in a conventional finite-difference formulation:

To study the stability, we consider $x_{j, k}^n = \kappa^n e^{\, i(j\, \theta + k\, \phi)}$. After substituting into the finite-difference equation, we solve for $\kappa$, which is the amplification factor. For stability, we require that $|\kappa |\le 1$. In two dimensions, the most restrictive case occurs when $\theta$ and $\phi$ take on values of $\pi$, resulting in $\cos(\theta) = \cos(\phi) = -1$. This leads to a condition on the time step:

The scheme is stable if $\Delta t$ fulfills this requirement. In practice, if $\Psi$ is not a constant, one often chooses as the maximum value of $\Psi$ or a sufficiently small $\Delta t$ to guarantee, to some extent, stability.

3.3.   Stability of the numerical scheme

The Von Neumann analysis, often known as Fourier analysis, is utilized to assess the stability of the scheme (3.5). This analysis is utilized only within the framework of linear systems. We first reduce the system (1.1) into a singular equation by utilizing the relationship $W_y = V_x$. Therefore, (1.1) is expressed as

Thus, we propose that the linear version is represented by

where the constant quantities $s_0 = \gamma W$, $s_1 = \delta W$, and $s_2 = \lambda V$ are given by,

As the Von Neumann stability analysis depends on a fixed spatial discretization, we substitute

We believe that our stability analysis is conservative; that is, if the scheme is stable for the largest and most restrictive step size, it will remain stable for all finer resolutions. Hence, we assume a solution form:

Consequently, we may write

Substituting (3.12) into (3.11) and doing some operations, yield

Hence,

where

Thus,

Equation (3.15) indicates the absence of growth. Thus, the stability condition of the Von Neumann analysis is satisfied. As a result, the scheme is unconditionally stable.

3.4.   Error analysis

Within the adaptive moving mesh framework, a non-uniform grid with local variable spacing is used for spatial discretization. Let the computational coordinates be $\eta_j$ and $\xi_k$ and local spacings be defined by

that vary with the indices $j$ and $k$ as well as with time $t$. Therefore, the local truncation error needs to be assessed based on these variable spacings.

Let us define the exact solution $W(\eta, \xi, t)$ of the transformed system and, for the adaptive grid point $(\eta_j, \xi_k)$ at time $t_n$, let us denote its numerical approximation as $W_{j, k}^{n}$. We describe the local error as follows:

The adaptive finite-difference scheme, which is comparable to (3.9), can be expressed in a general format as:

where $\mathcal{L}$ and $\mathcal{N}$ correspond to the discretized linear and, if available, nonlinear operators, respectively. Most importantly, these operators now involve finite differences calculated with the local spacings $\Delta \eta_j$ and $\Delta \xi_k$. For example, the spatial derivative with respect to $\eta$ can be expressed by a non-uniform formula:

with equivalent formulations for higher-order derivatives and for derivatives in $\xi$.

By substituting the exact solution into (3.17) and comparing it with the continuous equation, we define the local truncation error at $(\eta_j, \xi_k, t_{n+1})$ as

At this point, we use Taylor series expansions at $(\eta_j, \xi_k, t_{n+1})$, while carefully adjusting for the non-uniform mesh sizes. Thus, we derive

provided that the adaptive grid is smooth. In other words, if the variations of $\Delta \eta_{j}$ and $\Delta_{\xi, k}$ are both limited and smooth such that $\Delta_{\eta, j+1}-\Delta \eta_{j} = O(\Delta^2)$ (and likewise for $\Delta_{\xi, k}$), then the spatial discretization keeps its local second-order accuracy. Specifically, we can describe the truncation error as

where the constants $C_1$, $C_2$, and $C_3$ depend on higher order derivatives of $W$ together with the smoothness of the coordinate change in the computational and physical domains.

With respect to Table 2, the results show that the adaptive finite difference scheme has an estimated second order of spatial convergence. As the number of points increases at $t = 10$ with the relative $L_2$ error, the CPU time increases because of the additional computational effort required. The data in Figure 7 is represented in log-log form, with the left figure shows CPU time in relation to $\mathrm{Effective}\, L_x$. In contrast, the right shows the relative error as it pertains to $\mathrm{Effective}\, L_x$. The slope of the error curve validates the data shown in the right figure---meaning the scheme has achieved a second-order error in spatial accuracy of about $2$. Finer grids visually showcase the computational trade-off in the left panel. All these conclusions conform to the previous theoretical analysis, showcasing that the adaptive strategy achieves optimization through accuracy and the mesh with significant gradients.

4.   Conclusions

Using a generalized direct algebraic approach, we established exact solutions to the (2+1) dimensional generalized breaking soliton system (GBSS). According to the results of our analysis, the system accommodates a wide range of waveforms, including solitary and periodic forms. These results indicate that when non-linearity is balanced with dispersion, the GBSS serves as an effective model for intricate wave processes in fluid dynamics, plasma physics, and non-linear optics, which greatly depend on the wave phenomena. Additionally, we outlined a numerical scheme using the adaptive moving mesh techniques with the MMPDE5 framework. The adaptive approach effectively focuses grid points at steep gradient locations, accurately resolving solution features with minimized computations. Stability and error analyses confirm that the numerical method converges in space with second-order accuracy. Taken together, the analytical and the numerical approaches yield helpful information regarding the GBSS dynamics and provide a basis for further study of non-linear wave processes.

Author contributions

Amer Ahmed: Conceptualization, methodology, formal analysis, validation, writing original draft; A. R. Alharbi (arharbi@taibahu.edu.sa): Conceptualisation, investigation, data curation, visualisation; Haza S. Alayachi (hsshareef@taibahu.edu.sa): Investigation, funding acquisition; Ishak Hashim (ishak_h@ukm.edu.my): Supervision, writing review and editing. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The author, Amer Ahmed, is grateful to Taibah University for its generous support in the form of a fully paid scholarship, which significantly helped in successfully completing this work.

Conflict of interest

The authors declare no conflict of interest.