Analytical and numerical investigation of exact wave patterns and chaotic dynamics in the extended improved Boussinesq equation

2.1 Modified Kudryashov method

The modified Kudryashov method is an enhanced technique used to solve nonlinear evolution equations, building upon the traditional Kudryashov method [33]. The central improvement in the modified version involves replacing the natural exponential base e with an arbitrary base A [32,34]. Such a modification enhances the ability of problem solvers to deal with some nonlinear equations as more flexibility is provided by the freedom to choose an arbitrary base to suit the structure of the equation. Within this framework, a solution of a nonlinear evolution equation, for example, of Eq. (7), is constructed as a polynomial of the form:

where r i are the constant values, ϒ ( ζ ) is a necessary function that remains to be solved for, and m is the polynomial degree of the polynomial. The function ϒ ( ζ ) satisfies the following differential equation:

where A is a constant that must satisfy two conditions, A > 0 and A ≠ 1 . This step provides a generalization of a common function used in the process of nonlinear differential equations, expanding the resolution’s reach outside the conventional principles. The solution to this differential Eq. (10) is sought

where v is an arbitrary constant. This form allows the nonlinear evolution equation to take a more manageable algebraic form. Using the modified Kudryashov method with m = 2 , the solution can be written as follows:

When inserting this solution into the initial equation and using the characteristics of ϒ ( ζ ) , the problem is reduced to a problem of a system of algebraic equations in terms of r 0 , r 1 , and r 2 . The solutions of these equations are the exact solutions of the nonlinear problem, thus confirming the efficiency of the modified Kudryashov method.

This technique from a physicist perspective is also offering much flexibility that can be of much importance in systems where the behavior or the dynamics do not conform to traditional exponential perspective. When modeling physical systems, it is often seen that wave propagation and reaction–diffusion processes contain scaling laws, which make arbitrary base A more intuitive. This extra parameter can be useful to incorporate elaborate solutions in modeling nonlinear complex phenomena, and hence, the technique is of great use in several applications. For a certain class of nonlinear systems, characterized by reacting amplitude-dependent logarithmic terms, we derive precise traveling wave solutions.

• 1) he analysis begins with the assumption a 0 = 0 thus putting aside possible constant contributions to the system dynamics and focusing attention on the nonlinear behavior responsible by the values a 1 and a 2 . The constants a 1 and a 2 are expressed as follows:

  a 1 = − 12 ln 2 ( A ) ln 2 ( A ) − 1 , a 2 = 12 ln 2 ( A ) ln 2 ( A ) − 1 ,

  which capture the effect of the amplitude A on the system’s oscillatory modes. These parameters demonstrate a nonlinear dependence on A , indicating the system’s logarithmic responsiveness to changes in amplitude. The traveling wave s speed is determined by the following relationship:

  s = − 1 1 − ln 2 ( A ) ,

  showing that the amplitude and wave speed have a nontrivial dependency relation. The limiting speed of the wave increases without bound as ln 2 ( A ) approaches 1; this is typically understood as a point of a critical parameter in the system where change in wave behavior is likely to occur. Such a critical behavior points toward the possibility of phase transitions or bifurcations in the response of the system. The corresponding exact traveling wave solution is given by

  (13) Θ 1 ( x , t ) = − 12 v ln 2 ( A ) A s t + x ( ln 2 ( A ) − 1 ) ( v A s t + x + 1 ) 2 .

  The solution achieved here features strongly nonlinear traveling waves, whereby the amplitude and factors associated with its propagation are logarithmic functions of A . From the form of the solution, however, it can be expected that such waves may be either amplified or damped depending on the coupling of the nonlinear terms. In particular, the squared denominator implies that some parts of the space-time could be dominated by strong damping of the oscillations, while other parts like the poles could promote further the oscillations, provided that the conditions at the poles are right.

• This derived solution has its importance in explaining wave propagation phenomena within the nonlinear system and especially the critical or singular behavior of that system as the amplitudes approach predetermined thresholds. Such solutions have applications in various branches of physics such as hydrodynamics, plasmoids, nonlinear optics, among many others, where the relationship between amplitude speed and the nonlinearity is indeed a key factor in determining the responses of the system to various measures. It accommodates a stable and iterative wave form over certain bounded ranges where it is dominant the behavior of the system is related to certain parameters. For example, where ln 2 ( A ) is less than the critical value of 1, it implies that the interpreted wave speed is very well defined; hence, the wave propagates in a consistent or rather orderly manner. This is a specific kind of waveform, in which the patterns are consistent and motion is stable. As the value of ln 2 ( A ) approaches 1, the wave is associated with chaotic dynamics. This is the point where the critical wave speed becomes undefined, thus pointing to transition points of the wave in the system. From the viewpoint of “chaotic dynamics,” it may be viewed as a point that marks a transitory phase between order and chaos. Chaos appears when the wave starts becoming oscillatory with no fixed manner – in which the wave may or may not subsequently get “turned” or pass through the peak or nadir in reference to the time and the space relationship to the amplitude.

• 2) For the case where a 0 = − 2 ln 2 ( A ) ln 2 ( A ) + 1 , a 1 = 0 , and a 2 = 0 , the traveling wave equation simplifies significantly. The wave speed s is

  s = − 1 ln 2 ( A ) + 1 ,

  signifying that there is a nonlinear relationship between the wave speed and the amplitude A would repeat the sign. Given the aforementioned characteristics, the form provides an accurate and easy description for the motion of the wave within a nonlinear medium.

• The corresponding solution for the traveling wave is

  (14) Θ 2 ( x , t ) = − 2 ln 2 ( A ) ln 2 ( A ) + 1 .

  This solution describes a wave and explains that its amplitude is a logarithmic function of A , which does not change with time or space. Since the solution does not vary with either x or t , it indicates that the wave configuration is uniform (stationary) and resolution of space and time variables is not required. Here, a stable wave configuration is emphasized and is determined only by amplitude A . It depicts a stable wave where the wave is solely a function of time or space. This precise wave pattern represents a stationary wave, which behaves consistently. Despite the stability presented in the solution, chaotic dynamics might develop when the amplitude A is at a critical value or when small disturbances occur. Chaotic dynamics of nonlinear systems are often state dependent.

• 3) For the case where

  a 0 = − 2 ln 2 ( A ) ln 2 ( A ) + 1 , a 1 = 12 ln 2 ( A ) ln 2 ( A ) + 1 , a 2 = − 12 ln 2 ( A ) ln 2 ( A ) + 1 , and s = − 1 ln 2 ( A ) + 1 .

  The corresponding solution for the traveling wave is

  (15) Θ 3 ( x , t ) = − 2 ln 2 ( A ) ( v 2 A 2 ( − s t + x ) − 4 v A − s t + x + 1 ) ( ln 2 ( A ) + 1 ) ( v A − s t + x + 1 ) 2 .

  This solution captures the intricacy of wave interactions between A , which includes the amplitude, space x , and time t . The characteristic of the wave and its velocity of propagation is affected by even small variations in the amplitude A , which also exposes the degree of sensitivity of the system to any change in A . Under specific conditions, these dynamics can produce stable and predictable wave patterns known as “exact wave patterns.” However, any small variations in the initial conditions or amplitude A can lead to chaotic behavior, where the wave’s motion becomes irregular and unpredictable.

• 4) In the case where a 0 = − 2 ln 2 ( A ) ln 2 ( A ) + 1 , a 1 = 12 ln 2 ( A ) ln 2 ( A ) + 1 , and a 2 = − 12 ln 2 ( A ) ln 2 ( A ) + 1 , the exact solution for the traveling wave is derived as follows. The wave speed is given by

  s = 1 ln 2 ( A ) + 1 ,

  which demonstrates the nonlinear dependence of the wave speed on the amplitude A .

• The corresponding solution for the traveling wave is

  (16) Θ 4 ( x , t ) = − 2 ln 2 ( A ) ( v 2 A 2 x − 4 v A s t + x + A 2 s t ) ( ln 2 ( A ) + 1 ) ( A s t + v A x ) 2 .

  This solution describes complex wave dynamics with nonlinear interactions between the amplitude A , spatial position x , and time t . The exponential terms and the multiple factors involving A highlight that the wave’s shape and propagation are highly sensitive to changes in A , x , and t , leading to intricate wave behavior depending on the specific values of these parameters.

Figure 1 represents the analytical solutions to a certain equation, obtained using the modified Kudryashov method. These 2D plots offer insights into the behavior of the studied system with varying time and spatial variables. Each colored curve in both plots represents the evolution of an analytical wave solution over time and space. The horizontal axis represents the position x , while the vertical axis represents the analytical solution of the wave Θ 3 ( x , t ) at specific time values t . In Figure 1 (left), the parameters are given as A = 0.7 , v = 20 , and x 0 = − 4 , with time ranging from t = 0 to t = 15 , the number of points N = 1,000 , and space ranging from x = − 20 to x = 20 . Physically, this figure reflects an exact wave pattern, where the wave’s movement through space and time is regular and stable. All the curves show that the system exhibits predictable behavior, where the wave gradually rises and reaches a steady value. This indicates a stable, predictable wave behavior known as exact wave patterns in nonlinear systems. In Figure 1 (right), the parameters are changed to A = 0.5 and v = 200 . It can be noted that the right curves have a larger dispersion and a greater degree of variability than the left one, which is indicative the higher velocity v causes more irregular and intricate oscillations. This may be indicative of the beginning of chaos. The increase of the velocity v leads to greater nonlinear coupling of the amplitude and spatial modulation, which manifests in chaotic dynamics in which the waves undergo irregular oscillations of varying periods.

Figure 1

Modified Kudryashov method was used to obtain 2D plots of the analytical solutions to Eq. (4): (left) The parameters are given by A = 0.7 , v = 20 , x 0 = − 4 with t = 0 → 15 , N = 1,000 , and x = − 20 → 20 , and (right) the right plot shows the behavior when A = 0.5 , v = 200 .

2.2 IGREM method

In this section, we present the IGREM approach, which is quite effective in solving nonlinear differential equations [36–38]. The motivation for this approach is to develop a consistent and practical technique that is useful in addressing complex equations, which are common in many applied mathematics areas. We shall investigate the structure as well as the usage of the IGREM technique in details. The IGREM method is used on the equation to rewrite Eq. (7) as

where v i are the arbitrary constants. The function Φ ( ζ ) satisfies the following nonlinear differential equation:

From a general perspective, it is a scaled Riccati equation. In the case when ρ = 0 and γ = 0 , Eq. (17) is in the simplest form of a Riccati or Bernoulli equation, which is a standard equation in many mathematical models.

2.2.2 Traveling wave solutions

Taking into account that the homogeneous balance is m = 2 , the final form of the traveling wave solution can be expressed as

(18) Θ ˆ ( ζ ) = v 0 + v 1 Φ ( ζ ) + v 2 Φ ( ζ ) 2 .

Substituting Eq. (18) into Eq. (4) and simplifying, we derive a system of algebraic equations. Solving this system yields the corresponding soliton solutions. The exact traveling wave solutions are obtained for different families, depending on the discriminant μ = γ 2 − 4 ρ σ and associated conditions:

1. When μ = γ 2 − 4 ρ σ > 0 and ρ ≠ 0 , the solutions are divided into two cases:

Case 1:

(19) Θ 5 ( x , t ) = − 3 μ sech 2 1 2 μ ( x 0 − s t + x ) μ − 1 ,

Case 2:

(20) Θ 6 ( x , t ) = 6 i μ ( μ − 1 ) ( sinh ( μ ( x 0 − s t + x ) ) − i ) .

Family 2: When μ = γ 2 − 4 ρ σ < 0 , ρ σ ≠ 0 , the solutions are divided into two cases:

Case 3:

(21) Θ 7 ( x , t ) = − 3 μ sec 2 1 2 − μ ( x 0 − s t + x ) μ − 1 .

Case 4:

(22) Θ 8 ( x , t ) = − 3 μ csc 2 1 2 − μ ( x 0 − s t + x ) μ − 1 .

The IGREM method provides an effective framework for deriving exact traveling wave solutions to nonlinear differential equations. By examining different cases based on the discriminant μ , we obtain explicit analytical solutions that can be applied to various physical models exhibiting soliton behavior.

Figures 3(a) and 4(a) display the 3D and 2D plots for Θ 5 , respectively, demonstrating how the solutions evolve over the time interval from 0 to 20.

Figure 2

IGREM method was used to find the analytical solution for Θ 5 ( x , t ) . The study parameters were as follows: ρ = 2.5 , σ = 1.3 , γ = 3.6417 , x 0 = − 8 with t = 0 → 20 , m x = 8,000 and x = − 40 → 20 .

Figure 3

IGREM method was used to find the analytical solution for Θ 5 ( x , t ) and Ω 5 ( x , t ) , as shown in (a) and (c) as 3D surface plots. (b) and (d) Adaptive moving mesh numerical solution results. The study parameters were as follows: ρ = 2.5 , σ = 1.3 , γ = 3.6417 , x 0 = − 8 with t = 0 → 20 , m x = 8,000 , and x = − 40 → 20 .

Figure 2 depicts the analytical solution Θ 5 ( x , t ) obtained using the IGREM method, showcasing exact wave patterns characteristic of solitary waves (solitons) in nonlinear systems. In the 3D plot, the wave moves along the x -axis, preserving its localized structure over time, indicating energy conservation and stability. The 2D plot displays the wave profile at different time steps, further emphasizing the wave’s stability and lack of deformation. These solutions highlight the precise nature of solitons, where the interplay between nonlinearity and dispersion allows the wave to maintain its shape during propagation. However, in nonlinear systems, chaotic dynamics can emerge under specific conditions, such as perturbations or variations in system parameters, leading to irregular and unpredictable wave behavior. This transition from exact solutions to a simultaneous understanding of many phenomena is vital to the theory of nonlinear waves in a shallow water area.

Figure 4

IGREM method was used to find the analytical solution for Θ 5 ( x , t ) and Ω 5 ( x , t ) , as shown in (a) and (c) as 2D plots. (b) and (d) Adaptive moving mesh numerical solution results. The study parameters were as follows: ρ = 2.5 , σ = 1.3 , γ = 3.6417 , x 0 = − 8 with t = 0 → 20 , m x = 8,000 , and x = − 40 → 20 .

3.1 Numerical solutions using a uniform mesh

In this section, a uniform mesh method is applied to a physical domain of size [ 0 , L ] to obtain numerical solutions for the system in Eq. (24). The domain is divided into m x subintervals, each of uniform width Δ x = L m x , with the points denoted as [ x i , x i + 1 ] , where x i = ( i − 1 ) Δ x for all x i ∈ [ 0 , L ] and i = 1 , 2 , 3 , … , m x + 1 . The discretization of system (24) uses finite difference operators with a uniform subinterval width Δ x , and is expressed as follows:

the associated boundary conditions for system (24) are

The initial conditions are given by

It is important to note that this numerical system is solved using the FORTRAN ODE solver known as the DASPK solver [35].

3.2 Adaptive mesh numerical solutions

An adaptive method of a mesh is described to be employed in this section to address the numerical problem of mapping from the computational domain [ 0,1 ] to the physical domain [ 0 , L ] . The mapping is written in the form of x = x ( η , t ) , where η : [ 0,1 ] → [ 0 , L ] for t > 0 . In this formulation, the following coordinates introduced a physical coordinate, x , and a computational coordinate, η . The variables are defined as

where the association x = x ( η , t ) applies. The adaptive mesh refinement splits the geometrical area from x 0 = 0 to x m = L into m + 1 sub-intervals. The mesh x i will be given by x i = x ( η i , t ) , where η i = i m and i = 0,1 , … , m . Controlled differentiation allowed us to obtain the following:

where Θ x x is computed as

Boundary conditions are set as follows:

Initial conditions are based on

The methodology focuses on minimizing errors and maximizing convergence through different sets of monitoring functions with MMPDE6 [39–41] as the preferred method, which is given as

where 0 < τ < 1 is a relaxation parameter and ϖ ( Θ , Ω , x ) is a curvature-dependent monitoring function targeting regions with significant solution variations:

The user-specified weight parameters g 1 and g 2 allocate additional nodes to areas with high error. The approach includes discretizing the spatial derivative while maintaining a continuous temporal derivative. The coupled equations are solved using FORTRAN software and ODE solvers, employing finite differences within the adaptive mesh framework to achieve accurate solutions.

The IGREM method results, visualized through 3D surface plots, assess the dynamics of the functions Θ ( x , t ) and Ω ( x , t ) in relation to exact wave patterns and chaotic dynamics. Analytical solutions (Figure 3(a) and (c)) demonstrate precise, theoretical wave patterns indicative of stable and predictable system behaviors, which are crucial for verifying the accuracy of numerical methods. The numerical solutions (Figure 3(b) and (d)), utilizing an adaptive moving mesh, reflect the system’s sensitivity to initial conditions, characteristic feature of chaotic dynamics. The variable mesh density is strategically employed to better capture complex behaviors in areas of high gradient. The observable discrepancies between the analytical and numerical solutions suggest the presence of chaotic dynamics, particularly evident in abrupt changes within the numerical plots. This highlights the method’s dual capability to capture both the predictable and intricate behaviors of the modeled system, illustrating the complexities involved in simulating theoretical precision and chaotic dynamics within practical computational frameworks.