Shallow Water Waves with Surface Tension by Laplace–Adomian Decomposition

1 Introduction

The theory of shallow water waves is a long-standing topic of interest in the fluid dynamics community ever since the solitary waves were discovered during a horse-drawn carriage in the canal of Edinburgh. Several mathematical models have been proposed, and they have all been analyzed over time. The first model that was extensively studied in this context is the Korteweg-de Vries (KdV) equation, which was followed by the modified KdV equation. Subsequently, many models, to address shallow water wave dynamics, emerged with time. They are the Boussinesq equation, Kawahara equation, Camassa-Holm equation and several others.

Shallow water wave dynamics, though traditionally rooted in geophysical and coastal engineering studies, hold significant relevance in materials science and fluid-structure interaction (FSI). These wave models capture essential features of dispersive and nonlinear wave propagation, making them valuable tools for understanding wave-induced stresses and responses in flexible or deformable structures. In materials science, such insights are crucial when evaluating the mechanical behavior of thin films, membranes, and soft composites subjected to dynamic fluid environments.

In the context of FSI, shallow water wave models help simulate the interaction between wave fronts and material surfaces, revealing critical information about fatigue, resonance, and structural deformation under transient loads. For instance, micro- and nano-scale materials used in sensors or biomedical devices may interact with fluid interfaces in ways that mirror shallow wave behavior, particularly when capillarity and surface tension play a role. Additionally, the incorporation of surface tension effects into shallow water wave equations enhances the applicability of these models to scenarios involving liquid–solid boundaries, such as in microfluidic channels or soft robotics.

The effects of surface tension have been incorporated into just a small number of the most recent models that were previously mentioned. This aspect was first taken up during 2002 and 2006 [1,2], respectively. This was for the Boussinesq equation. Later the model was extensively studied analytically with the effect of surface tension included. It must be noted that the surface tension effect leads to the inclusion of sixth-order dispersion that extends the familiar Boussinesq equation that comes with fourth-order dispersion. The solitary wave and shock wave solution to this sixth-order Boussinesq equation was recovered by the aid of undetermined coefficients during 2023 [3]. The conservation laws were also identified using the multipliers approach. Subsequently, the perturbed version of the sixth-order Boussinesq equation was addressed and its solitary waves and shock waves were retrieved and reported during 2025 [4].

Moving further along, after the analytical studies, it is now time to take a look at this model from a numerical perspective. The current paper does exactly that. This sixth-order Boussinesq equation is studied numerically in the current paper with the aid of the Laplace-Adomian decomposition method (LADM). This integration scheme is an extension of the pre-existing versions of the Adomian decomposition, namely the improved Adomian decomposition method. This improved version of the Adomian decomposition scheme was successfully implemented in a couple of models from Quantum Optics [5,6].

In the current paper, the LADM integration scheme retrieves solitary waves and shock waves. The error analysis of this scheme is also displayed, for solitary waves as well as shock waves, and the observation is that it is impressively low, of the order of 10−8. The details are chalked out in the rest of the paper after a quick revisitation of the model and a succinct recapitulation of the scheme.

2 The Model of the Sixth-Order Boussinesq Equation

The Boussinesq equation characterizing the motion of long waves traveling bidirectionally in shallow water under the influence of gravity is expressed as qtt=qxx+aqxxxx+(q2)xx, (1) where a is a non-zero constant. This equation is well-posed when a = −1, but it is ill-posed when a = 1 [7,8,9]. An enhanced Boussinesq equation that is well-posed is represented by

A number of variants of Eqs. (1) and (2) have been examined in the literature. A generalization implies replacing the nonlinear term q2 with a nonlinear differential function f(q). A further generalization involved accounting for high dispersion components and incorporating surface tension effects. These generalizations result in the equation

The independent variables are the spatial variable x and the temporal variable t, whereas the dependent variable, representing the wave structure, is represented as q(x, t). The wave operator is represented by the first two terms, where k corresponds to the wave number. The parameter m is the general power law parameter, and the nonlinearity coefficient is denoted by c. This parameter provides the model with a generalized nature. In addition, the fourth-order dispersion parameters, which are caused by the effects of ai, are connected in the original model from the Navier–Stokes equation. On the other hand, the coefficients of bi come from the surface tension effect.

Eq. (3) is referred to as the Sixth-Order Boussinesq Equation (6BE) and was formulated based on the basic principles in the year 2006 in [1]. Later, during 2023, several analytical studies with the model have been conducted such as the retrieval of the analytical form of solitary and shock waves, recovery of the conservation laws and such [3]. The current paper addresses the study of shallow water waves, with the effect of surface tension included, numerically, for the first time.

SBE is used in various mathematical models for real-world situations, like studying vibrations in nonlinear atomic chains [10], solving some microstructure problems [11], and examining the movement of inviscid flow in shallow layers of fluid [2]. An interesting study on waves and solitons can be found at [12]. These applications highlight the versatility of the SBE, as it provides insights into complex phenomena across different scientific fields. Furthermore, ongoing research continues to explore its potential in new areas, enhancing our understanding of both theoretical and practical aspects of fluid dynamics.

Eq. (3) has been extensively studied, analyzed, and resolved using a variety of methodologies in the specific case of m = 1 [13,14,15,16,17,18]. The results suggest that the fundamental ideas are more profoundly understood as a consequence of the consistent behavior observed across various approaches. Additionally, these findings open the door to the investigation of complex scenarios in which m ≠ 1. In the current paper, we introduce a novel method for solving (3) when the parameter value is m ≠ 1. The proposed method is the Adomian decomposition method improved with the combination of the Laplace transform.

3 Waves Type Solutions

This section will illustrate the general form of both solitary waves and shock waves, which are solutions to the 6BE. The solution structure of such waves inherently generates the parametric constraints. The subsequent information is also listed and presented.

3.1 Solitary Wave Type Solutions

The authors in [3,4] have only recently found the solitary wave type solutions for Eq. (3) using the undetermined coefficients technique, which is expressed as q(x,t)=A1sechβ[B1(x−vt)], (4) where A1 denotes the wave amplitude, B1 denotes the inverse width, and v is the corresponding wave velocity, the parameter β is determined to be β=22m−1,(5) from which m ≠ 1/2.

The amplitude A1 of the wave is A1=(2m+1)(b1+b2k2)2cb22m−1,(6) from which the restriction naturally arises:

The velocity of the wave is given by v=−b1b2,(8) with the restriction

This demonstrates that the 6BE will have solitary waves as long as the two sixth-order dispersion terms have opposing signatures.

The inverse width B1 of the wave is B1=|2m−1|2b1+b2k2a1b2−a2b1,(10) with the restriction (b1+b2k2)(a1b2−a2b1)>0,(11) another quick observation leads to m > 1/2.

Consequently, the nonlinearity of the power law indicates the presence of solitary waves, provided that the constraint on m is supported.

3.2 Shock Wave Type Solutions

The authors in [3,4] have only recently found the shock wave type solutions for Eq. (3) using the undetermined coefficients technique, which is expressed as q(x,t)=A2tanhβ[B2(x−vt)],(12) where A2 denotes the wave amplitude, B2 denotes the inverse width, and v is the corresponding wave velocity, the parameter β takes the value of β = 1, from which m = 3/2.

The amplitude A2 of the wave is A2=−(b1+b2k2)c,(13) from which the restriction naturally arises:

The velocity of the wave is given by v=−b1b2,(15) with the restriction

This confirms that the 6BE will have solitary waves as long as the two sixth-order dispersion terms have opposing signatures.

The inverse width B2 of the wave is B2=b1+b2k22(a2b1−a1b2),(17) with the restriction

As a consequence, the nonlinearity of the power law suggests the existence of shock waves, provided that the constraint m = 3/2 is satisfied.

4 The Laplace-Adomian Decomposition Method

This section shows the general procedure for a numerical treatment of Eq. (3), which is based on the initial conditions provided for solitary waves and shock waves. The present scenario employs the Laplace–Adomian decomposition method (LADM), originally established by Refs. [19,20]. This method helps to break down complicated equations into simpler parts, making it easier to find accurate solutions for both solitary and shock wave situations.

Before establishing the steps of the method, we can write Eq. (3) as: qtt=k2qxx−c(q2m)xx−a1qxxxx−a2qxxtt−b1qxxxxxx−b2qxxxxtt,(19) which, in terms of differential operators, can be expressed as Dttq(x,t)=(L+N)q(x,t),(20) where D2t represents the second time derivative and L is the linear differential operator acting on q, defined by Lq(x,t)=(k2D2x−a1D4x−a2D2x2t−b1D6x−b2D4x2t)q(x,t),(21) where the symbolism is D i x j t = ∂ i + j ∂ i x ∂ j t .

In addition, in (20), N symbolizes a nonlinear operator acting on q by

The following are the steps that are involved in the procedure of the method:

Step 1:

Apply the temporal Laplace transform L to both sides of Eq. (19) to obtain

With the assumption that L { q ( x , t ) } = q ( x , S ) and the known properties associated with Laplace transform differentiation, we derive

The values of q(x, 0), qxx(x, 0), q4x(x, 0), etc., are calculated from the initial conditions associated with the waveforms listed in Section 3.

Step 2:

The function that is unknown q may be expressed as an infinite set of summands by using the typical Adomian decomposition method: q(x,t)=∑n=0∞qn(x,t),(25) where each component qn(x, y) must be determined iteratively.

Furthermore, the nonlinear term N is expressed as an infinite series of Adomian polynomials in a single variable, Nq(x,t)=∑n=0∞Qn(q0,…,qn),(26) where Qn denotes all Adomian polynomials [21].

The Adomian polynomials Qn in a single variable may be derived for all types of nonlinearity, computed using the following formulae [21]:

These polynomials provide a systematic way to handle nonlinear terms in differential equations, enabling the construction of a series solution that converges to the desired function. By applying these polynomials, one can effectively separate the linear and nonlinear components, simplifying the analysis and solution of complex problems.

Substituting Eqs. (27) and (26) into Eq. (25) gives

Utilizing the linearity of the Laplace transform in Eq. (28) and equating both sides of the equation, we get the subsequent recursive formula: qm(x,S)=1S2L{(k2D2x−a1D4x−b1D6x)qm−1(x,t)+Qm−1(q0,…,qn)},m≥1,(29) where q0 is given by

Step 3:

We will consider the inverse Laplace transform L − 1 and the formula L − 1 1 S n = t n − 1 ( n − 1 ) ! . Subsequently, applying L − 1 to Eq. (30) yields the following:

Similarly, applying the inverse Laplace transform in Eq. (29) to obtain:

Eqs. (31) and (32) allow us to determine each term of the function sequence { q n } n = 0 ∞ .

The solution derived from the suggested approach is expressed as an infinite series, represented by:

The deduced series solution may converge to an exact solution, required that such a solution exists. Alternatively, the series solution may be employed for numerical applications. For further information concerning the convergence of the suggested technique, consult Refs. [22,23].

As a result, a truncated series solution with N components is obtained as

The series solution (34) works well for numerical applications.

The subsequent section presents numerical examples to demonstrate the algorithm’s high accuracy and efficiency as offered by the suggested technique.

5 Numerical Simulation of Wave-Type Solutions

This section will simulate wave solutions for Eq. (3) utilizing the algorithm explained in Section 4. We will derive the initial conditions for each simulation from the general waveforms presented in Section 3. These initial conditions will provide the necessary parameters to accurately model the wave behavior as described in the equations. By referencing the general waveforms, we ensure that our simulations are grounded in established theoretical frameworks. The Mathematica 14.2.1 software is employed to generate the numerical simulation results.

5.1 Simulation of Solitary Wave Solutions and Graphical Representations

We will consider the three nonlinear cases in Table 1 when simulating solitary waves using LADM. It is crucial to keep in mind that m > 1/2 is one of the conditions for the existence of solitary waves in Section 3.1. This analysis is essential for validating the theoretical predictions and ensuring the reliability of the simulation results.

To implement the algorithm offered by LADM via Eqs. (31) and (32), we must take into account the general solitary waveforms illustrated in Section 3.1; that is, we must consider:

Moreover, using Formula (27), the first Adomian polynomials to decompose the nonlinear term q3/2 in the case m = 3/4 are given by:

These polynomials serve as the foundation for constructing a convergent series solution to the nonlinear problem at hand.

Also, using Formula (27), the first Adomian polynomials to decompose the nonlinear term q3 in the case m = 3/2 are given by:

Finally, the first Adomian polynomials to decompose the nonlinear term q5 in the case m = 5/2 are given by:

The graphical results of the solitary wave simulation using the algorithm developed by LADM are illustrated in Fig. 1, Fig. 2 and Fig. 3, which display the three cases enumerated in Table 1. The simulations illustrate the behavior of solitary waves in a selection of conditions, emphasizing the significant differences in propagation speed and amplitude. Additionally, we have included a comparison of the error generated by the Adomian decomposition method (ADM) in each figure. This method produces a slightly higher level of error than the LADM. The data used in the simulations conducted by the ADM is identical to that used in the LADM simulations. These visual representations offer valuable insights into the dynamics of wave actions and their implications for future research.

Figure 1: (a) 3D solitary wave simulated by LADM, (b) 2D density depiction of the corresponding wave, and (c) 2D time evolution for different time points. (d) error in the simulation by LADM, and (e) error in the simulation by ADM. We consider a scenario where m = 3/4 and N = 15 steps are used.

Figure 2: (a) 3D solitary wave simulated by LADM, (b) 2D density depiction of the corresponding wave, and (c) 2D time evolution for different time points. (d) error in the simulation by LADM, and (e) error in the simulation by ADM. We consider a scenario where m = 3/2 and N = 15 steps are used.

Figure 3: (a) 3D solitary wave simulated by LADM, (b) 2D density depiction of the corresponding wave, and (c) 2D time evolution for different time points. (d) error in the simulation by LADM, and (e) error in the simulation by ADM. We consider a scenario where m = 5/2 and N = 15 steps are used.

5.2 Simulation of Shock Wave Solutions and Graphical Representations

We will consider the three nonlinear cases in Table 2 when simulating shock waves using LADM. It is necessary to consider that the conditions established in Section 3.2 require the existence of shock waves for 6BE only if β = 1. Consequently, the sole supported nonlinearity is q3, which corresponds to m = 3/2. This analysis is essential for validating the theoretical predictions and ensuring the reliability of the simulation results.

Also, using Formula (27), the first Adomian polynomials to decompose the nonlinear term q3 in the case m = 3/2 are the same as those obtained in the second simulation example in Table 1.

The graphical results of the solitary wave simulation using the algorithm developed by LADM are illustrated in Fig. 4, Fig. 5 and Fig. 6, which display the three cases enumerated in Table 2. The simulations illustrate the behavior of solitary waves in a selection of conditions, emphasizing the significant differences in propagation speed and amplitude. Additionally, we have included a comparison of the error generated by the Adomian decomposition method (ADM) in each figure. This method produces a slightly higher level of error than the LADM. The data used in the simulations conducted by the ADM is identical to that used in the LADM simulations. These visualisations show wave action dynamics and its implications for future study.

Figure 4: (a) 3D shock wave simulated by LADM, (b) 2D density depiction of the corresponding wave, and (c) 2D time evolution for different time points. (d) error in the simulation by LADM, and (e) error in the simulation by ADM. For the case (i) of Table 2 with N = 15 steps.

Figure 5: (a) 3D shock wave simulated by LADM, (b) 2D density depiction of the corresponding wave, and (c) 2D time evolution for different time points. (d) error in the simulation by LADM, and (e) error in the simulation by ADM. For the case (ii) of Table 2 with N = 15 steps.

Figure 6: (a) 3D shock wave simulated by LADM, (b) 2D density depiction of the corresponding wave, and (c) 2D time evolution for different time points. (d) error in the simulation by LADM, and (e) error in the simulation by ADM. For the case (iii) of Table 2 with N = 15 steps.

6 Conclusions

The current paper retrieves solitary waves and shock waves that emerged from the sixth-order Boussinesq equation, which models shallow water wave dynamics in presence of surface tension. The Laplace-Adomian decomposition scheme has made this recovery possible with a stunningly low error measure. The solitary waves and shock wave solutions are numerically addressed using this integration algorithm. The sound results of the work pave the way for future research in this direction with the model.

Later, the shallow water wave dynamics with multi-layer fluids will be addressed. The corresponding vector coupled model of the sixth-order Boussinesq equation will be formulated and subsequently analyzed analytically as well as numerically. The results of such research activities will be disseminated over time.