Exact traveling wave and soliton solutions for chemotaxis model and (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation

Natarajan Annapoorani; Nagarajan Keerthana; Kottakkaran Sooppy Nisar

doi:10.1515/phys-2025-0204

Artikel Open Access

Exact traveling wave and soliton solutions for chemotaxis model and (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation

Natarajan Annapoorani , Nagarajan Keerthana und Kottakkaran Sooppy Nisar

Veröffentlicht/Copyright: 17. September 2025

Veröffentlicht von

Veröffentlichen auch Sie bei De Gruyter Brill

Manuskript einreichen Informationen für Autor*innen Erkunden Sie dieses Fachgebiet

Aus der Zeitschrift Open Physics Band 23 Heft 1

Abstract

This study investigates exact traveling wave solutions for two important nonlinear models: the hyperbolic chemotaxis model and the (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation. Using the G ′ G 2 -expansion method, we derive solutions in the form of trigonometric, exponential, and rational functions, showcasing diverse wave behaviors such as periodic solitons, kink-type waves, and singular solitons. By visualizing these solutions with 3D and contour plots, we provide deeper insights into the spatial–temporal evolution of the wave dynamics. The findings highlight the versatility of the G ′ G 2 -expansion method for solving complex nonlinear equations and demonstrate its relevance in both biological systems, like chemotaxis, and physical wave dynamics.

Keywords: chemotaxis; partial differential equation; traveling wave solution; (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation

1 Introduction

Traveling wave solutions of nonlinear partial differential equations (NLPDEs) have attracted considerable attention in recent years, owing to their fundamental role in describing wave-like phenomena and signal propagation across a wide range of physical and biological systems. Over recent decades, researchers have increasingly focused on understanding the inherent properties and quantitative features of NLPDEs to better describe wave phenomena in fields such as fluid dynamics, plasma physics, and biological modeling. In particular, traveling waves play a critical role in explaining the movement and interactions of organisms in chemotaxis.

This growing interest in NLPDEs has led to significant research dedicated to finding exact solutions. Advances in symbolic computation have enabled the development of various powerful techniques for obtaining these solutions. For example, the Hirota direct method [1] is renowned for its effectiveness in solving integrable equations, while the homogeneous balance method [2] is commonly used to derive solitary and periodic wave solutions. The tanh -function method [3] is especially effective for higher-dimensional models, and the auxiliary equation method [4,5] provides a systematic approach for identifying exact solutions in complex systems. Additionally, the three-wave approach [6] has been valuable in constructing multi-soliton solutions, offering insights into wave interactions.

In recent years, significant progress has been made in the development of advanced analytical techniques aimed at uncovering the intricate structures and long-time dynamics of nonlinear evolution equations. Notably, the application of the linear Fokas unified transform method has led to the discovery of a novel long-range instability phenomenon associated with the inhomogeneous linear Schrödinger equation posed on the vacuum spacetime quarter-plane [7]. The Dbar-steepest descent method has proven effective in characterizing the long-time asymptotic behavior of the Wadati–Konno–Ichikawa equation, resulting in a rigorous resolution of the soliton resolution conjecture and the establishment of asymptotic stability for its solutions [8,9]. Moreover, the integration of bi-Hamiltonian structures and recursion operators has facilitated the proof of nonlinear stability for exact smooth multi-soliton solutions within the framework of the two-component Camassa–Holm system [10]. These advancements illustrate the evolving landscape of mathematical methodologies employed in the study of nonlinear wave propagation across diverse physical and geometric contexts.

In this article, we focus on studying traveling wave solutions for the hyperbolic chemotaxis model and the (3+1)-dimensional Boiti–Leon–Manna–Pempinelli (BLMP) equation.

Chemotaxis, the directed movement of organisms in response to chemical gradients, has been extensively studied in mathematical biology, especially concerning bacterial behavior. In 1971, Keller and Segel [11] introduced a foundational model that describes bacterial aggregation and chemotactic pattern formation. Their work formulated a system of partial differential equations, to represent the movement of E. coli in response to nutrient gradients, revealing phenomena such as bacterial band formation. This model highlighted the importance of understanding traveling wave solutions, which has since become a central focus in chemotaxis research.

Expanding upon the Keller–Segel model, Hillen and Stevens [12,13] introduced a hyperbolic model that considers the finite speed of organism movement, particle densities, and interactions with an external chemical signal. Their model incorporates turning rates influenced by the concentration of the chemical signal and its spatial gradient, as well as a reaction term to describe the chemical signal’s production and degradation. This approach offers a more detailed perspective on the dynamics underlying chemotactic behavior:

(1) u t + + ( γ ( S , S t , S x ) u + ) x = − μ + ( S , S t , S x ) u + + μ − ( S , S t , S x ) u − , u t − − ( γ ( S , S t , S x ) u − ) x = μ + ( S , S t , S x ) u + − μ − ( S , S t , S x ) u − , τ S t = D 0 S x x + f ( S , u + + u − ) ,

where u ± ( x , t ) represents the particle densities of moving particles, while γ ( S , S t , S x ) denotes the speed of these particles, and μ ± ( S , S t , S x ) represents the turning rates. The concentration of the external chemical signal is given by S ( x , t ) , and its production and degradation are described by the reaction term f .

In the following years, Rivero et al. [14] and Ford et al. [15] applied Hillen and Stevens models to experimental data, providing numerical solutions that helped validate and refine the theoretical framework. These studies examined the spatiotemporal development of bacterial densities and the formation of chemotactic patterns, offering insights into the stability of these patterns under various conditions. Concurrently, Hillen and Levine [16] investigated finite-time blow-up in chemotaxis systems, where solutions could display blow-up behavior based on specific reaction rates and the characteristics of the turning functions.

The investigation of traveling wave solutions in chemotaxis models remains a prominent focus within mathematical biology. Liu [17] derived explicit solutions for a particular configuration of the system, incorporating reaction terms and turning rates consistent with the Hillen and Stevens framework. This study underscores the importance of constructing traveling waves, which are critical for comprehending the long-term dynamics of bacterial aggregation and chemotactic pattern formation. In a special case of Eq. (1), the system is reformulated into an equivalent structure involving the total particle density u = u + + u − and the particle flux v = u + − u − . The resulting system is as follows:

(2) u t + γ v x = 0 , v t + γ u x = γ D 1 ω w x w u − γ 2 D 1 v , w t = D 0 w x x + w u ,

where the total particle density u ( x , t ) represents the combined density of particles moving in both directions, while the particle flux v ( x , t ) denotes the net flux of particles, with v ( x , t ) = 0 indicating no net movement. The chemical signal concentration w ( x , t ) describes the concentration of the external chemical signal. The particle speed γ reflects the finite speed at which particles move, while the diffusion coefficient D 0 indicates the rate at which the chemical signal diffuses. The particle diffusion constant D 1 characterizes the diffusion behavior of particles. The interaction parameter ω defines the sensitivity of particle movement to the chemical gradients.

The BLMP equation is a notable PDE that has attracted considerable interest due to its broad applications in fields such as plasma physics, fluid dynamics, ocean engineering, astrophysics, and aerodynamics. Initially introduced by Boiti et al. [18], this equation effectively captures the dynamics of complex wave phenomena by incorporating both nonlinear interactions and higher-order effects. As such, it has proven to be a crucial tool for modeling wave propagation in incompressible fluids under various conditions. Particularly, the (3+1)-dimensional form of the BLMP equation plays a significant role in describing wave dynamics within incompressible media. When z = 0 , it provides valuable insights into the propagation of Riemann waves, enhancing our understanding of wave behavior in both the presence and absence of specific boundary conditions [19]. Given its ability to model intricate wave interactions, the BLMP equation has become an essential component of nonlinear wave theory and fluid mechanics, contributing to analyses of complex wave dynamics in systems such as optical fibers, plasma environments, and fluid models.

Darvishi [20] introduced the (3+1)-dimensional BLMP equation given by

(3) q y t + q z t + q x x x y + q x x x z − 3 q x ( q x y + q x z ) − 3 q x x ( q y + q z ) = 0

and studied multi-wave solutions for this equation using the multiple exp-function method. Exact solutions, including rational solutions, soliton solutions, positons, and negatons, were obtained via the Wronskian technique [21]. Lax pairs, Bäcklund transformations, and multi-soliton solutions with two distinct dispersion relations, along with a bilinear form using a general logarithmic transformation, were derived for the equation [22–27]. Additionally, three-wave solutions, such as kinky periodic solitary-wave solutions, periodic soliton solutions, and kink solutions, have also been presented [28].

Special cases of Eq. (3) are identified as follows:

When y = z and q ( x , y , t ) = q ( θ , t ) = q [ x + ψ ( y ) , t ] , with ψ as an analytic function of y , Eq. (3) reduces to the Korteweg–de Vries equation [29]:
(4) q θ t + q θ θ θ θ − 6 q θ q θ θ = 0 ,
which is widely used to model shallow water waves with long wavelengths and small amplitude, as observed in fluid and plasma dynamics [30].
When y = z = ρ , t = τ , ψ = ς , and the variable transformation σ ′ = σ ρ is applied, Eq. (3) takes the form [19]:
(5) σ τ ′ + σ ς ς ς ′ − 3 ( σ ′ σ ″ ) ς = 0 , σ ″ = σ ς ,
where σ ′ and σ ″ represent the components of velocity in an incompressible fluid. Eq. (5) is recognized as a generalized model for incompressible fluids and can be seen as an extension of a (2+1)-dimensional system discussed in previous studies [19,31].

In a more complex adaptation of the (3+1)-dimensional BLMP Eq. (3), Wazwaz [32] introduced an additional spatial derivative q x into the derivative term q y + q z , yielding a new form of the equation:

(6) ( q x + q y + q z ) t + α ( q x + q y + q z ) x x x + β ( q + x ( q x + q y + q z ) ) x = 0 ,

where q ( x , y , z , t ) represents an unknown analytical function defined over the spatial variables x , y , and z and temporal variable t , while α and β are non-zero constants.

Only a few soliton solutions and their interactions have been explored in previous research using methods such as the Painlevé test and Hirota’s direct method [33–37]. This study focuses on deriving a limited set of exact solutions for the (3+1)-dimensional BLMP equation using the G ′ G 2 -expansion method. These solutions, including lump solitons and solitary waves, provide fresh insights into their interactions and dynamics. By analyzing these solutions in the context of chemotaxis and other nonlinear systems, we offer new perspectives on traveling wave behaviors. This work contributes to the understanding of multi-soliton dynamics, presenting specific results through detailed 3D and contour visualizations, thus extending the current literature with these novel findings.

To the best of the author’s knowledge, this is the first study to derive a comprehensive set of exact traveling wave solutions, including all types of soliton solutions, for both the hyperbolic chemotaxis model and the (3+1)-dimensional BLMP equation using the G ′ G 2 -expansion method [38,39]. Unlike traditional methods, which often struggle to capture the full spectrum of solution types, the G ′ G 2 -expansion method enables the derivation of hyperbolic, rational, and singular solitons with greater precision and flexibility. This approach facilitates the construction of explicit and accurate expressions for traveling waves, thereby advancing the understanding of soliton dynamics. In this work, we apply the method to derive the complete set of soliton solutions for Eqs. (2) and (6), significantly broadening the range of exact solutions in nonlinear wave theory. Building upon our previous work [40], which focused on the existence and asymptotic behavior of traveling waves in the chemotaxis model, we extend this analysis to derive a full spectrum of soliton solutions, marking a significant advancement in nonlinear wave modeling.

The structure of this article is as follows: Section 2 provides an in-depth explanation of the G ′ G 2 -expansion method, outlining its theoretical foundation and procedural approach for application to NLPDEs. In Section 3, the method is implemented to derive exact traveling wave solutions for the hyperbolic chemotaxis model and the (3+1)-dimensional BLMP equation. Section 4 discusses the results obtained, providing numerical validation of the theoretical solutions and exploring the behavior of the solutions under various parameter conditions. Finally, Section 5 concludes the article by summarizing the key findings, discussing their implications, and highlighting the broader impact of the G ′ G 2 -expansion method.

2 Methodology: The G ′ G 2 -expansion method

This algorithm details the G ′ G 2 expansion method, a systematic approach tailored for solving NLPDEs of the form:

(7) P 1 ( h , h t , h x , h t t , h x x , h x t , … ) = 0 ,

where P 1 represents a polynomial involving the unknown function h = h ( x , t ) and its various partial derivatives with respect to independent variables x and t . The method outlined follows the subsequent procedure:

Transformation into ordinary differential equations (ODEs): Initiate the process by transforming (7) into an ODE, denoted as follows:

(8) P 2 ( H , H ′ , H ″ , … ) = 0 ,

where P 2 is a polynomial involving H and its derivatives concerning a new variable z . Employ a transformation suited for NLPDEs, expressed as follows:

(9) h ( x , t ) = H ( z ) , z = m x − n t ,

where m and n are the constants. Given the non-linear ODE

(10) G ′ G 2 ′ = μ + λ G ′ G 2 2 ,

where μ and λ are the integers.

Distinct scenarios for solutions: The solutions of (10) can be categorized into three distinct scenarios:

If μ λ > 0 ,
(11) G ′ G 2 = μ λ C cos ( μ λ z ) + D sin ( μ λ z ) D cos ( μ λ z ) − C sin ( μ λ z ) .
If μ λ < 0 ,
(12) G ′ G 2 = − ∣ μ λ ∣ λ C sinh ( 2 ∣ μ λ ∣ z ) + C cosh ( 2 ∣ μ λ ∣ z ) + D C sinh ( 2 ∣ μ λ ∣ z ) + C cosh ( 2 ∣ μ λ ∣ z ) − D .
If μ = 0 and λ ≠ 0 ,
(13) G ′ G 2 = − C λ ( C z + D ) ,
where C and D are the arbitrary non-zero constants.

Series solution in terms of G ′ G 2 polynomial: Utilize the G ′ G 2 expansion method to express the series solution G ( z ) of Eq. (8) in terms of the G ′ G 2 polynomial:
(14) H ( z ) = a 0 + ∑ j = 1 N a j G ′ G 2 j + b j G ′ G 2 − j ,
where a 0 , a j , and b j ( j = 1,2 , … , N ) are the constants.
Determination of positive integer N : Ascertain the positive integer N in Eq. (14) through the application of the homogeneous balance principle.
Systematic polynomial transformation and coefficient equations: Insert the expressions (14) and (10) into (8), yielding a polynomial in terms of G ′ G 2 . Equate the coefficients of the resulting polynomials to zero, establishing a solvable system of algebraic equations for the unknown constants present in the series solution. This system can be effectively solved using symbolic computational packages like Maple or Mathematica. The explicit exact solutions of (7) can be obtained by inserting the values of a 0 , a j , b j ( j = 1 , … , N ) , m , n and the solutions in (11)–(13) into (14) with the transformation in (9).

3 Implementation of the G ′ G 2 -expansion method

3.1 Hyperbolic chemotaxis model

By excluding v from the system of Eq. (2), it becomes apparent that the resulting equation can be formulated as follows:

(15) − 1 γ u t t + γ u x x = ω γ D 1 u w x w x + γ D 1 u t , w t = D 0 w x x + w u .

Considering the transformation

(16) u = u ( z ) , w = e η ( z ) , z = m ( x − n t ) + z 0 ,

with m > 0 , n > 0 , and z 0 representing any constant, Eq. (15) can be expressed as follows:

(17) m γ − n 2 γ u ″ = ω γ m D 1 ( u ′ η ′ + u η ″ ) − γ n D 1 u ′ ,

(18) − m n η ′ = D 0 m 2 ( ( η ′ ) 2 + η ″ ) + u .

From the second equation of (17), we have

(19) u = − m n η ′ − D 0 m 2 ( η ′ ) 2 − D 0 m 2 η ″ ,

(20) u ′ = − m n η ″ − 2 D 0 m 2 η ′ η ″ − D 0 m 2 η ‴ ,

(21) u ″ = − m n η ‴ − 2 D 0 m 2 ( η ″ ) 2 − 2 D 0 m 2 η ′ η ‴ − D 0 m 2 η ″ ″ .

By substituting equations from (19) to (21), in the first equation of (17), we have

(22) − m 2 n γ η ‴ − 2 D 0 m 3 γ ( η ″ ) 2 − 2 D 0 m 3 γ η ′ η ‴ − D 0 m 3 γ η ″ ″ + m 2 n 3 γ η ‴ + 2 γ D 0 m 3 n 2 ( η ″ ) 2 + 2 γ D 0 m 3 n 2 η ′ η ‴ + D 0 γ m 3 n 2 η ″ ″ + ω D 1 γ n m 2 η ′ η ″ + 2 D 1 D 0 ω m 3 γ ( η ′ ) 2 η ″ + D 0 D 1 m 3 ω γ η ′ η ‴ + ω D 1 γ m 2 n η ′ η ″ + D 0 D 1 ω m 3 γ ( η ′ ) 2 η ″ + D 0 D 1 ω m 3 γ ( η ″ ) 2 − m D 1 n 2 γ η ″ − 2 D 1 D 0 γ m 2 n η ′ η ″ − D 0 D 1 m 2 γ n η ‴ = 0 .

By simplifying and applying the homogeneous balance principle to the highest order derivative and the nonlinear term ⇒ N + 3 = N + ( N + 2 ) ⇒ N = 1 . Hence, the specific form of the solution of Eq. (15) is written as

(23) η ( z ) = a 0 + a 1 G ′ G 2 + b 1 G ′ G 2 − 1 ,

where a 0 , a 1 , and b 1 represent the undetermined coefficients. Upon substituting (23) into (22) together with (10), we proceed to group coefficients with identical powers of G ′ G 2 i (where i = 0 , ± 1 , ± 2 , … ). By then setting these coefficients to zero, we obtain the ensuing system of algebraic equations:

G ′ G 2 − 4 : − 6 D 0 m 3 γ b 1 2 μ 2 + 6 D 0 m 3 γ b 1 μ 3 + 2 γ D 0 m 3 n 2 b 1 2 μ 2 − 6 γ D 0 b 1 m 3 n 2 μ 3 − 3 D 1 D 0 ω m 3 γ b 1 3 μ + 3 D 1 D 0 ω m 3 γ b 1 2 μ 2 + 4 D 1 D 0 m 3 n 2 b 1 2 μ 2 ,

G ′ G 2 4 : − 6 D 0 m 3 γ a 1 2 λ 2 − 6 D 0 m 3 γ a 1 λ 3 + 2 γ D 0 m 3 n 2 a 1 2 λ 2 + 6 γ D 0 m 3 n 2 a 1 λ 3 + 4 D 1 D 0 m 3 n 2 a 1 2 λ 2 + 3 D 1 D 0 ω m 3 γ a 1 3 λ + 3 D 1 D 0 ω m 3 γ a 1 2 λ 2 ,

G ′ G 2 − 3 : − 4 D 0 m 3 γ a 0 b 1 μ 2 − 4 D 1 D 0 m 3 n 2 a 0 b 1 μ 2 − 6 D 1 D 0 m 3 ω γ a 0 b 1 2 μ + 2 D 1 D 0 m 3 ω γ a 0 b 1 μ 2 + 2 D 1 D 0 γ m 2 n b 1 2 μ + 2 D 1 D 0 m 2 γ n b 1 μ 2 − 2 m 2 n γ b 1 μ 2 + 2 γ m 2 n 3 b 1 μ 2 − 2 D 1 ω γ n m 2 b 1 2 μ ,

G ′ G 2 3 : − 4 D 0 m 3 γ a 0 a 1 λ 2 + 4 D 1 D 0 m 3 n 2 a 0 a 1 λ 2 + 6 D 1 D 0 m 3 ω γ a 0 a 1 2 λ + 2 D 1 D 0 m 3 ω γ a 0 a 1 λ 2 − 2 D 1 D 0 γ m 2 n a 1 2 λ − 2 D 1 D 0 m 2 γ n a 1 λ 2 − 2 m 2 n γ a 1 λ 2 + 2 γ m 2 n 3 a 1 λ 2 + 2 D 1 ω γ n m 2 a 1 2 λ ,

G ′ G 2 − 2 : − 8 D 0 m 3 γ b 1 2 λ μ + 8 D 0 m 3 γ b 1 λ μ 2 + 4 γ D 0 m 3 n 2 b 1 2 λ μ − 4 γ D 0 m 3 n 2 a 1 b 1 μ 2 + 4 D 1 D 0 m 3 n 2 b 1 2 μ λ − 8 γ D 0 m 3 n 2 b 1 λ μ 2 − 3 D 1 D 0 ω m 3 γ a 0 2 b 1 μ − 3 D 1 D 0 m 3 ω γ a 1 b 1 2 μ + 4 D 1 D 0 m 3 ω γ b 1 2 μ λ + 2 D 1 D 0 γ m 3 n a 0 b 1 μ − 2 D 1 ω γ m 2 n a 0 b 1 μ + m D 1 γ n 2 b 1 μ + 4 D 1 D 0 m 3 n 2 a 1 b 1 μ 2 − 3 D 1 D 0 ω m 3 γ b 1 3 λ ,

G ′ G 2 2 : − 8 D 0 m 3 γ a 1 2 λ μ − 8 D 0 m 3 γ μ λ 2 a 1 + 4 γ D 0 m 3 n 2 a 1 2 λ μ − 4 γ D 0 m 3 n 2 a 1 b 1 λ 2 + 4 D 1 D 0 m 3 n 2 a 1 b 1 λ 2 + 8 γ D 0 m 3 n 2 a 1 λ 2 μ + 3 D 1 D 0 m 3 ω γ a 0 2 a 1 λ + 3 D 1 D 0 ω γ m 3 a 1 2 b 1 λ + 3 D 1 D 0 ω γ m 3 a 1 3 μ + 4 D 1 D 0 m 3 ω γ a 1 2 μ λ − 2 D 1 D 0 γ m 3 n a 0 b 1 λ + 2 D 1 ω γ n m 2 a 0 b 1 λ − m D 1 n 2 γ a 1 λ + 4 D 1 D 0 m 3 n 2 a 1 2 μ λ ,

G ′ G 2 − 1 : − 4 D 0 m 3 γ a 0 b 1 λ μ + 4 D 1 D 0 m 3 n 3 a 0 b 1 λ μ − 6 D 1 D 0 ω m 3 γ a 0 b 1 2 λ − 2 m 2 n γ b 1 λ μ + 2 D 1 D 0 m 3 ω γ a 0 b 1 λ μ + 2 D 1 D 0 γ m 3 n b 1 2 λ − 2 D 1 D 0 m 2 γ n b 1 λ μ − 2 D 1 ω γ m 2 n b 1 2 λ + 2 γ m 2 n 3 b 1 λ μ ,

G ′ G 2 1 : − 4 D 0 m 3 γ a 0 b 1 λ μ + 4 D 1 D 0 m 3 n 2 a 0 a 1 λ μ + 6 D 1 D 0 ω m 3 γ a 0 a 1 2 μ − 2 m 2 n γ a 1 λ μ + 2 D 1 D 0 m 3 ω γ a 0 a 1 λ μ − 2 D 1 D 0 γ m 2 n a 1 2 μ − 2 D 1 D 0 m 2 γ n a 1 λ μ − 2 D 1 ω γ m 2 n a 1 2 μ + 2 γ m 2 n 3 a 1 λ μ ,

G ′ G 2 0 : − 2 D 0 m 3 γ a 1 2 μ 2 − 2 D 0 m 3 γ b 1 2 λ 2 − 2 D 0 m 3 γ a 1 λ μ 2 + 2 D 0 m 3 γ b 1 μ λ 2 + 2 γ D 0 m 3 n 2 a 1 2 μ 2 + 2 γ D 0 m 3 n 2 b 1 2 λ 2 − 8 γ D 0 m 3 n 2 a 1 b 1 λ μ + 8 D 1 D 0 m 3 n 2 a 1 b 1 μ λ + 2 γ D 0 m 3 n 2 a 1 λ μ 2 − 2 γ D 0 m 3 n 2 b 1 μ λ 2 + 3 D 1 D 0 ω m 3 γ a 0 2 a 1 μ − 3 D 1 D 0 ω m 3 γ a 0 2 b 1 λ + 3 D 1 D 0 ω m 3 γ a 1 2 b 1 μ − 3 D 1 D 0 ω m 3 γ a 1 b 1 2 λ − 2 D 1 D 0 γ m 2 n a 0 a 1 λ − m D 1 n 2 γ a 1 μ + 2 D 1 D 0 γ m 2 n a 0 b 1 μ + D 0 D 1 ω m 3 γ a 1 2 μ 2 + D 0 D 1 ω m 3 γ b 1 2 λ 2 + m D 1 n 2 γ b 1 λ + 2 D 1 ω γ n m 2 a 0 a 1 μ − 2 D 1 ω γ n m 2 a 0 a 1 λ .

By solving the aforementioned algebraic equations using Maple software, we obtained three sets of values for a 0 , a 1 , b 1 , n , and m . Utilizing these values, along with Eqs. (11) through (13) and substituting them into the solution form (23), we derive three distinct results.

Result 1:

When λ μ > 0 :
The solution involving a trigonometric function can be expressed as follows:
(24) u 1 1 ( x , t ) = a 0 ∓ 32 D 0 λ μ n 1 3 μ 2 D 1 + 3 μ λ 2 + n 2 μ × C cos ( λ μ z ) + D sin ( λ μ z ) D cos ( λ μ z ) − C sin ( λ μ z ) − 1 ,
where a 0 is a constant, and C and D are the parameters that determine the specific trigonometric behavior of the solution. The term λ μ reflects the square root of the product of λ and μ , affecting the oscillatory nature of the solution.
When λ μ < 0 :
The solution involving an exponential function is given by
(25) u 2 1 ( x , t ) = a 0 ∓ 64 D 0 λ μ n 1 3 μ 2 D 1 + 3 μ λ 2 + n 2 μ × 2 ∣ λ μ ∣ − 4 C ∣ λ μ ∣ e 2 z ∣ λ μ ∣ C e 2 z ∣ λ μ ∣ − D − 1 ,
where ∣ λ μ ∣ indicates the absolute value of the product λ μ , and the solution incorporates exponential terms, which imply growth or decay depending on the sign and magnitude of λ μ .
When λ ≠ 0 and μ = 0 :
The solution simplifies to a constant:
(26) u 3 1 ( x , t ) = a 0 .
This indicates that when μ is zero and λ is non-zero, the solution does not vary with x and t , leading to a steady-state or equilibrium solution.

Result 2:

When λ μ > 0 :
The alternative trigonometric function solution is
(27) u 1 2 ( x , t ) = a 0 ± 32 D 0 λ μ n 1 3 μ 2 D 1 + 3 μ λ 2 + n 2 μ × C cos ( λ μ z ) + D sin ( λ μ z ) D cos ( λ μ z ) − C sin ( λ μ z ) .
This expression provides a trigonometric solution without inversion, which means it directly involves the ratio of trigonometric functions.
When λ μ < 0 :
The solution involving an exponential function is given by
(28) u 2 2 ( x , t ) = a 0 ± 16 D 0 n 1 3 μ 2 D 1 + 3 μ λ 2 + n 2 μ × 2 ∣ λ μ ∣ − 4 C ∣ λ μ ∣ e 2 z ∣ λ μ ∣ C e 2 z ∣ λ μ ∣ − D .
This solution simplifies the previous exponential function by removing the inversion, leading to a direct exponential term involving e 2 z ∣ λ μ ∣ .
When λ ≠ 0 and μ = 0 : The rational function solution is
(29) u 3 2 ( x , t ) = a 0 ∓ 32 D 0 n D 1 C C z + D .
This shows a rational function that varies linearly with z , where the term C z + D appears in the denominator, affecting the shape of the solution.

Result 3:

When λ μ > 0 :
The trigonometric function solution for an alternative case is
(30) u 1 3 ( x , t ) = a 0 ± 32 D 0 λ μ n 1 7 D 1 + 9 λ μ 2 + 4 n × C cos ( λ μ z ) + D sin ( λ μ z ) D cos ( λ μ z ) − C sin ( λ μ z ) ∓ 32 D 0 λ μ n 1 7 D 1 + 9 λ μ 2 + 4 n × C cos ( λ μ z ) + D sin ( λ μ z ) D cos ( λ μ z ) − C sin ( λ μ z ) − 1 .
This result combines both the trigonometric function and its reciprocal, indicating a more complex solution structure.
When λ μ < 0 :
The solution involving an exponential function is given by
(31) u 2 3 ( x , t ) = a 0 ± 16 D 0 n 1 7 D 1 + 9 λ μ 2 + 4 n × 2 ∣ λ μ ∣ − 4 C ∣ λ μ ∣ e 2 z ∣ λ μ ∣ C e 2 z ∣ λ μ ∣ − D ∓ 64 D 0 μ λ n 1 7 D 1 + 9 λ μ 2 + 4 n × 2 ∣ λ μ ∣ − 4 C ∣ λ μ ∣ e 2 z ∣ λ μ ∣ C e 2 z ∣ λ μ ∣ − D − 1 .
This solution includes both direct and inverse exponential terms, reflecting more intricate behavior under the condition λ μ < 0 .
When λ ≠ 0 and μ = 0 :
The rational function solution is:
(32) u 3 3 ( x , t ) = a 0 ∓ 32 D 0 n D 1 C C z + D .
This result provides a rational function form similar to Result 2.

3.2 (3+1)-dimensional BLMP equation

By selecting α = 1 and β = − 3 in Eq. (6), we reduce the model to

(33) ( q x + q y + q z ) t + ( q x + q y + q z ) x x x − 3 ( q + x ( q x + q y + q z ) ) x = 0 .

To further analyze the traveling wave solutions, we introduce the transformation q ( x , y , z , t ) = Q ( ξ ) , where ξ = k 1 x + k 2 y + k 3 z + c t . This substitution transforms Eq. (6) into an ODE given by

(34) c ( k 1 + k 2 + k 3 ) Q ″ + k 1 3 ( k 1 + k 2 + k 3 ) Q 4 − 6 k 1 2 ( k 1 + k 2 + k 3 ) Q ′ Q ″ = 0 .

By integrating Eq. (34) with respect to ξ , we obtain the reduced form:

(35) c ( k 1 + k 2 + k 3 ) Q ′ + k 1 3 ( k 1 + k 2 + k 3 ) Q ‴ − 3 k 1 2 ( k 1 + k 2 + k 3 ) ( Q ′ ) 2 = 0 .

To solve Eq. (35), we apply the homogeneous balance principle to determine the order of the polynomial solution. This method involves balancing the highest-order nonlinear term with the highest-order derivative term. Through this process, we determine that N = 1 , indicating that the solution will have a specific polynomial form.

The explicit form of the solution to Eq. (33) is constructed in a manner analogous to the expression provided in Eq. (23). Substituting the solution form given in Eq. (23) into Eq. (35), along with the conditions derived from Eq. (10), we expand the resulting expression in terms of G ′ G 2 i , where i = 0 , ± 1 , ± 2 , … . By setting the coefficients of each power of G ′ G 2 i to zero, we generate a system of algebraic equations. Solving this system allows us to obtain the parameters required for the exact traveling wave solutions:

G ′ G 2 − 4 : − 6 k 1 4 b 1 μ 3 − 6 k 1 3 k 2 b 1 μ 3 − 6 k 1 3 k 3 b 1 μ 3 − 3 k 1 3 b 1 2 μ 2 − 3 k 1 2 k 2 b 1 2 μ 2 − 3 k 1 2 k 3 b 1 2 μ 2 ,

G ′ G 2 4 : 6 k 1 4 a 1 λ 3 + 6 k 1 3 k 2 a 1 λ 3 + 6 k 1 3 k 3 a 1 λ 3 − 3 k 1 3 a 1 2 λ 2 − 3 k 1 2 k 2 a 1 2 λ 2 − 3 k 1 2 k 3 a 1 2 λ 2 ,

G ′ G 2 − 2 : − c k 1 b 1 μ − μ c k 2 b 1 − c k 3 b 1 μ + 4 k 1 4 b 1 μ 2 λ − 2 k 1 4 b 1 μ 2 λ − 18 k 1 4 b 1 μ 2 λ + 8 k 1 4 b 1 μ 2 λ + 4 k 1 3 k 2 b 1 μ 2 λ − 2 k 1 3 k 2 b 1 μ 2 λ − 18 k 1 3 k 2 b 1 μ 2 λ + 8 k 1 3 k 2 b 1 μ 2 λ + 4 k 1 3 k 3 b 1 μ 2 λ − 2 b 1 k 1 3 k 3 μ 2 λ − 18 k 1 3 k 3 b 1 μ 2 λ + 8 k 1 3 k 3 b 1 μ 2 λ − 6 k 1 3 b 1 2 μ λ + 6 k 1 3 a 1 b 1 μ 2 − 6 k 1 2 k 2 b 1 2 μ λ + 6 k 1 2 k 2 a 1 b 1 μ 2 − 6 k 1 2 k 3 b 1 2 μ λ + 6 k 1 2 k 3 a 1 b 1 μ 2 ,

G ′ G 2 2 : c k 1 a 1 λ + c k 2 a 1 λ + c k 3 a 1 λ + 2 k 1 4 a 1 μ λ 2 + 6 k 1 4 a 1 μ λ 2 + 4 k 1 4 b 1 λ 3 − 6 k 1 4 b 1 λ 3 − 6 k 1 4 b 1 λ 3 + 8 k 1 4 b 1 λ 3 + 2 k 1 3 k 2 a 1 μ λ 2 + 6 k 1 3 k 2 a 1 λ 2 μ + 4 k 1 3 k 2 b 1 λ 3 − 6 k 1 3 k 2 b 1 λ 3 − 6 k 1 3 k 2 b 1 λ 3 + 8 k 1 3 k 2 b 1 λ 3 + 2 k 1 3 k 3 a 1 μ λ 2 + 6 k 1 3 k 3 a 1 μ λ 2 + 4 k 1 3 k 3 b 1 λ 3 − 6 b 1 k 1 3 k 3 λ 3 − 6 k 1 3 k 3 b 1 λ 3 + 8 k 1 3 k 3 b 1 λ 3 − 6 k 1 3 a 1 2 μ λ + 6 k 1 3 a 1 b 1 λ 2 − 6 k 1 2 k 2 a 1 2 μ λ + 6 k 1 2 k 2 a 1 b 1 λ 2 − 6 k 1 2 k 3 a 1 2 μ λ + 6 k 1 2 k 3 a 1 b 1 λ 2 ,

G ′ G 2 0 : c k 1 a 1 μ − c k 1 b 1 λ + c k 2 a 1 μ − c k 2 b 1 λ + c k 3 a 1 μ − c k 3 b 1 λ + 2 k 1 4 a 1 μ 2 λ + 4 k 1 4 b 1 μ λ 2 + 4 k 1 4 b 1 λ 2 μ − 2 k 1 4 b 1 μ λ 2 − 6 k 1 4 b 1 μ λ 2 − 18 k 1 4 b 1 μ λ 2 + 8 k 1 4 b 1 λ 2 μ + 8 k 1 4 b 1 λ 2 μ + 2 k 1 3 k 2 a 1 μ 2 λ + 4 k 1 3 k 2 b 1 μ λ 2 + 4 k 1 3 k 2 b 1 μ λ 2 − 2 k 1 3 k 2 b 1 μ λ 2 − 6 k 1 3 k 2 b 1 λ 2 μ − 18 k 1 3 k 2 b 1 μ λ 2 + 8 k 1 3 k 2 b 1 μ λ 2

+ 8 k 1 3 k 2 b 1 μ λ 2 + 2 k 1 3 k 3 b 1 μ 2 λ + 4 k 1 2 k 3 b 1 μ λ 2 + 4 k 1 2 k 3 b 1 μ λ 2 − 2 b 1 k 1 3 k 3 μ λ 2 − 6 b 1 k 1 3 k 3 μ λ 2 − 18 k 1 3 k 3 b 1 μ λ 2 + 8 k 1 3 k 3 b 1 μ λ 2 + 8 k 1 3 k 3 b 1 μ λ 2 − 3 k 1 3 a 1 2 μ 2 − 3 k 1 3 b 1 2 λ 2 + 6 k 1 3 a 1 b 1 μ λ + 6 k 1 3 a 1 b 1 μ λ − 3 k 1 2 k 2 a 1 2 μ 2 − 3 k 1 2 k 3 b 1 2 λ 2 + 6 k 1 2 k 2 a 1 b 1 μ λ + 6 k 1 2 k 2 a 1 b 1 μ λ − 3 k 1 2 k 3 a 1 2 μ 2 − 3 k 1 2 k 3 b 1 2 λ 2 + 6 k 1 2 k 3 a 1 b 1 μ λ + 6 k 1 2 k 3 a 1 b 1 λ μ .

By solving the algebraic equations derived from the methodology, we utilized Maple to compute three distinct sets of parameter values for a 0 , a 1 , b 1 , k 1 , k 2 , k 3 , and c . Each set of these parameters corresponds to a specific configuration of the traveling wave solution. By substituting these parameter values into the expressions outlined in equations (11) to (13) and further integrating them into the solution form provided in Eq. (23), we successfully constructed three separate and well-defined solutions. These solutions highlight the diverse nature of the wave structures that can emerge under the given conditions.

Result 1:

When λ μ > 0 :
The solution in terms of a trigonometric function is given by
(36) q 1 1 ( x , y , z , t ) = a 0 − 2 μ c μ λ × C cos ( λ μ ξ ) + D sin ( λ μ ξ ) D cos ( λ μ ξ ) − C sin ( λ μ ξ ) − 1 .
When λ μ < 0 :
The solution that includes an exponential function can be represented as
(37) q 2 1 ( x , y , z , t ) = a 0 − μ c λ × 2 ∣ λ μ ∣ − 4 C ∣ λ μ ∣ e 2 ξ ∣ λ μ ∣ C e 2 ξ ∣ λ μ ∣ − D − 1 .
When λ ≠ 0 and μ = 0 :
The solution reduces to a constant expression:
(38) q 3 1 ( x , y , z , t ) = a 0 .

Result 2:

When λ μ > 0 :
The solution involving trigonometric functions can be written as
(39) q 1 2 ( x , y , z , t ) = a 0 + 2 λ c μ λ × C cos ( λ μ ξ ) + D sin ( λ μ ξ ) D cos ( λ μ ξ ) − C sin ( λ μ ξ ) .
Here, this expression directly includes trigonometric ratios without inversion.
When λ μ < 0 :
The corresponding solution in terms of an exponential function is
(40) q 2 2 ( x , y , z , t ) = a 0 + c × 2 ∣ λ μ ∣ − 4 C ∣ λ μ ∣ e 2 ξ ∣ λ μ ∣ C e 2 ξ ∣ λ μ ∣ − D .
When λ ≠ 0 and μ = 0 :
A rational function solution is obtained as follows:
(41) q 3 2 ( x , y , z , t ) = a 0 − c C C ξ + D .

Result 3:

When λ μ > 0 :
The trigonometric solution for another scenario is given by
(42) q 1 3 ( x , y , z , t ) = a 0 + 2 λ c μ λ C cos ( λ μ ξ ) + D sin ( λ μ ξ ) D cos ( λ μ ξ ) − C sin ( λ μ ξ ) − 2 μ c μ λ C cos ( λ μ ξ ) + D sin ( λ μ ξ ) D cos ( λ μ ξ ) − C sin ( λ μ ξ ) − 1 .
When λ μ < 0 :
The exponential function solution is expressed as
(43) q 2 3 ( x , y , z , t ) = a 0 + c 2 ∣ λ μ ∣ − 4 C ∣ λ μ ∣ e 2 ξ ∣ λ μ ∣ C e 2 ξ ∣ λ μ ∣ − D − μ c λ 2 ∣ λ μ ∣ − 4 C ∣ λ μ ∣ e 2 ξ ∣ λ μ ∣ C e 2 ξ ∣ λ μ ∣ − D − 1 .
When λ ≠ 0 and μ = 0 :
The rational form of the solution is
(44) q 3 3 ( x , y , z , t ) = a 0 − c C C ξ + D .

4 Discussion and results

This section presents the graphical representations of the obtained results, illustrating the dynamical properties of the solutions through 3D surface plots and contour graphs. These visualizations reveal the breadth of solutions derived using the G ′ G 2 -expansion method, demonstrating its efficacy in identifying diverse soliton structures that have not been previously explored. For both the hyperbolic chemotaxis model and the (3+1)-dimensional BLMP equation, the method effectively produces a range of distinct solutions, including periodic solitons, singular periodic solitons, kink-type waves, and singular solitons, achieved through appropriate selection of free parameters. This diversity underscores the robustness of the G ′ G 2 -expansion method in addressing the complexity of nonlinear systems, offering novel and significant contributions to the study of soliton dynamics in mathematical physics and applied nonlinear science.

For the hyperbolic chemotaxis model, the periodic soliton solutions derived from Eqs. (24) and (42) are computed with the parameter values λ = 1.1 , μ = 2 , a 0 = 1 , D 0 = 0.5 , D 1 = 0.3 , n = 1 , C = 1 , and D = 2 within the spatial and temporal range x = t = − 10 to 10. These solutions illustrate the interaction between particle density and the chemical signal, revealing a periodic oscillatory behavior typical of soliton dynamics.

The graphical representations in Figures 1 and 2, which include both 3D and contour views, illustrate the periodic soliton solutions showing stable oscillations. These oscillations reflect recurring patterns of concentration and flux, arising from the dynamic interaction with the chemical signal. The periodic nature of the system’s evolution is highlighted, with wave-like behavior emphasizing the stable oscillatory dynamics that govern the movement of particles and the diffusion of the chemical signal.

$Figure 1 3D representation (a) and contour graph (b) of the periodic wave solution for u 1 1 ( x , t ) {u}_{1}^{1}\left(x,t) .$

Figure 1

3D representation (a) and contour graph (b) of the periodic wave solution for u 1 1 ( x , t ) .

$Figure 2 3D representation (a) and contour graph (b) of the periodic wave solution for u 1 3 ( x , t ) {u}_{1}^{3}\left(x,t) .$

Figure 2

3D representation (a) and contour graph (b) of the periodic wave solution for u 1 3 ( x , t ) .

For Eq. (27), the solution shown in Figure 3, with the parameters λ = 1.5 , μ = 2 , a 0 = 1 , D 0 = 0.6 , D 1 = 0.4 , n = 2 , C = 1 , and D = 2 within the spatial and temporal range x = t = − 10 to 10, demonstrates a singular periodic behavior. This solution is characterized by sharp transitions and discontinuities in the concentration and flux profiles. Unlike conventional periodic solutions that show smooth oscillations, this one reveals an intricate structure with abrupt changes. These features indicate strong nonlinear interactions within the system, where the dynamic interplay between particle movement and chemical signaling leads to such discontinuities. The observed behavior highlights the complexity and intensity of the underlying dynamics in the chemotaxis model.

$Figure 3 3D representation (a) and contour graph (b) of the singular periodic wave solution for u 1 2 ( x , t ) {u}_{1}^{2}\left(x,t) .$

Figure 3

3D representation (a) and contour graph (b) of the singular periodic wave solution for u 1 2 ( x , t ) .

For Eqs. (28) and (43), the singular solution is obtained by considering the parameters λ = − 1.2 , μ = 1.5 , D 0 = 0.7 , D 1 = 0.5 , n = 3 , C = 3 , and D = 4 for Eq. (28), and λ = − 0.8 , μ = 1.3 , D 0 = 0.6 , D 1 = 0.4 , n = 2 , C = 2 , and D = 5 for Eq. (43). These specific values lead to a singular solution, characterized by discontinuities or sharp transitions in the concentration and flux patterns, differing significantly from smooth or periodic solutions. The singularity reflects abrupt changes in the system, such as jumps in particle density or chemical signal strength, and highlights the intensity of non-linear interactions within the model. Graphical results in Figures 4 and 5 further emphasize these features, showcasing the distinct behavior of singular solutions in chemotaxis models. These solutions shed light on the dynamic processes governing particle movement and chemical signaling in biological systems, particularly in scenarios with rapid transitions and non-continuous patterns.

$Figure 4 3D representation (a) and contour graph (b) of the singular wave soliton for u 2 2 ( x , t ) {u}_{2}^{2}\left(x,t) .$

Figure 4

3D representation (a) and contour graph (b) of the singular wave soliton for u 2 2 ( x , t ) .

$Figure 5 3D representation (a) and contour graph (b) of the singular wave soliton for u 2 3 ( x , t ) {u}_{2}^{3}\left(x,t) .$

Figure 5

3D representation (a) and contour graph (b) of the singular wave soliton for u 2 3 ( x , t ) .

For Eq. (25), the solution shown in Figure 6, with the parameters λ = − 0.3 , μ = 0.5 , D 0 = 0.7 , D 1 = 0.4 , n = 2 , C = 3 , and D = 1 , exhibits a kink solution. Kink solutions feature a sharp, non-smooth transition, where the particle density shifts abruptly between two stable states. This distinct profile emerges due to the specific interaction of the chemical signal and particle density under the given conditions, emphasizing the system’s nonlinear nature. The kink solution is particularly notable in systems where localized changes in concentration occur over a narrow region, reflecting a critical response to parameter changes such as λ and μ . This behavior is characteristic of systems displaying strong nonlinearity, where the steady-state solutions correspond to abrupt shifts in concentration, similar to what is observed in phase transitions or other systems with well-defined interfaces.

$Figure 6 3D representation (a) and contour graph (b) of the kink wave soliton for u 2 1 ( x , t ) {u}_{2}^{1}\left(x,t) .$

Figure 6

3D representation (a) and contour graph (b) of the kink wave soliton for u 2 1 ( x , t ) .

In the case of the (3+1) dimensional BLMP equation, the periodic solutions derived from equations (36) and (42), with the parameters a 0 = 1 , μ = 0.5 , k = 3 , λ = 1 , C = 1 , D = 2 , k 1 = 1 , k 2 = 1.5 , and k 3 = 0.5 , depict wave-like phenomena that exhibit periodic behavior over both spatial and temporal dimensions. Graphically, these solutions display oscillatory patterns, where the interaction between the cosine and sine functions (dependent on λ μ ) governs the periodicity and amplitude of the waves. The specific parameters chosen control the wave speed, intensity, and modulation, resulting in periodic structures that repeat at regular intervals. These solutions illustrate the dynamic nature of wave propagation and interference within the system, where the interplay of the nonlinear terms yields a complex but predictable periodic pattern, as demonstrated by the periodic oscillations observed in Figures 7 and 8.

$Figure 7 3D representation (a) and contour graph (b) of the periodic wave solution for q 1 1 ( x , y , z , t ) {q}_{1}^{1}\left(x,y,z,t) .$

Figure 7

3D representation (a) and contour graph (b) of the periodic wave solution for q 1 1 ( x , y , z , t ) .

$Figure 8 3D representation (a) and contour graph (b) of the periodic wave solution for q 1 3 ( x , y , z , t ) {q}_{1}^{3}\left(x,y,z,t) .$

Figure 8

3D representation (a) and contour graph (b) of the periodic wave solution for q 1 3 ( x , y , z , t ) .

For Eq. (39), the solution exhibits singular periodic behavior when the following parameters are used: a 0 = 1 , λ = 2 , μ = 2.5 , C = 1 , D = 1.5 , c = 1 , k 1 = 0.5 , k 2 = 1 , and k 3 = 1.5 . The singular periodic solutions shown in Figure 9 exhibit periodic oscillations interspersed with singularities, where the wave becomes undefined or infinitely large at specific points. These singularities result from the interaction of nonlinear terms with trigonometric functions, leading to sharp features or discontinuities in the wave structure. This behavior is a characteristic of nonlinear wave dynamics, where the wave retains its periodic nature but exhibits singular growth at certain locations.

$Figure 9 3D representation (a) and contour graph (b) of the singular periodic wave solution for q 1 2 ( x , y , z , t ) {q}_{1}^{2}\left(x,y,z,t) .$

Figure 9

3D representation (a) and contour graph (b) of the singular periodic wave solution for q 1 2 ( x , y , z , t ) .

In the context of Eqs. (37) and (43), the kink solutions as shown in Figures 10 and 11, display exponential growth and decay, marked by sharp transitions between two distinct states. With parameters a 0 = 1 , λ = − 1 , μ = − 0.5 , C = 2 , D = 1 , c = 1 , k 1 = 0.3 , k 2 = 1 , and k 3 = 0.7 , these solutions model localized transitions that can arise in systems exhibiting nonlinearity.

$Figure 10 3D representation (a) and contour graph (b) of the kink wave soliton for q 2 1 ( x , y , z , t ) {q}_{2}^{1}\left(x,y,z,t) .$

Figure 10

3D representation (a) and contour graph (b) of the kink wave soliton for q 2 1 ( x , y , z , t ) .

$Figure 11 3D representation (a) and contour graph (b) of the kink wave soliton for q 2 3 ( x , y , z , t ) {q}_{2}^{3}\left(x,y,z,t) .$

Figure 11

3D representation (a) and contour graph (b) of the kink wave soliton for q 2 3 ( x , y , z , t ) .

For Eq. (40), the singular solution can be obtained with the following parameter values: a 0 = 1 , λ = − 1 , μ = 0.5 , C = 1 , D = 1 , c = 1 , k 1 = 0.3 , k 2 = 1 , and k 3 = 0.7 . Figure 12 illustrates a localized wave-like structure with an amplitude that decays exponentially at both ends. Unlike oscillatory waves, it does not propagate periodically but remains confined to a specific region. The wave shape reflects a balance between exponential growth and decay, influenced by the nonlinear terms in the equation and the selected parameters.

$Figure 12 3D representation (a) and contour graph (b) of the singular wave soliton for q 2 2 ( x , y , z , t ) {q}_{2}^{2}\left(x,y,z,t) .$

Figure 12

3D representation (a) and contour graph (b) of the singular wave soliton for q 2 2 ( x , y , z , t ) .

5 Conclusion

In this study, we applied the G ′ G 2 -expansion method to derive exact traveling wave solutions for both the hyperbolic chemotaxis model and the (3+1)-dimensional BLMP equation. By selecting appropriate parameter values, we demonstrated the method’s capability to generate a range of soliton solutions, including periodic solitons, singular periodic solitons, kink-type waves, and singular solitons, each with distinct characteristics reflecting the complex dynamics of nonlinear systems.

The periodic solutions exhibit oscillatory behavior, while singular periodic and kink solutions represent abrupt transitions and localized disturbances, which are critical for understanding sharp variations in concentration or flux. These solutions emphasize the diversity of nonlinear wave phenomena, driven by the interplay of nonlinearities in the models.

This study highlights the effectiveness of the G ′ G 2 -expansion method in solving NLPDEs, proving its versatility in extracting traveling wave solutions from complex mathematical models. The broad spectrum of solution types obtained underscores the method’s ability to address various nonlinear wave behaviors, further enhancing our understanding of soliton dynamics.

The findings have important implications for both theoretical and applied fields, particularly in biological modeling, where traveling wave solutions are essential for describing processes like chemotaxis, as well as in other complex physical systems governed by nonlinear dynamics. This research contributes to the expanding field of NPDEs, offering new perspectives on wave propagation and opening up new directions for future exploration.

Acknowledgments

This study was supported via funding from Prince Sattam bin Abdulaziz University Project Number (PSAU/2025/R/1446).

Funding information: The funding from Prince Sattam bin Abdulaziz University Project Number (PSAU/2025/R/1446).
Author contributions: Conceptualization: NK; formal analysis: KSN; investigation: NK, NA; methodology: NA; software: NK, NA, KSN; validation: KSN; writing – original draft: NK, NA, KSN; writing – review editing: NA, KSN. All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The authors state no conflict of interest.
Data availability statement: All data generated or analysed during this study are included in this published article.

References

[1] Wazwaz A-M. Soliton solutions for two (3+1)-dimensional non-integrable KdV-type equations. Math Comput Model. 2012;55:1845–8. 10.1016/j.mcm.2011.11.082. Suche in Google Scholar

[2] Fan E. Two new applications of the homogeneous balance method. Phys Lett A. 2000;265:353–7. 10.1016/S0375-9601(00)00010-4. Suche in Google Scholar

[3] Parkes EJ, Duffy BR. An automated tanh-function method for finding solitary wave solutions to non-linear evolution equations. Comput Phys Commun. 1996;98:288–300. 10.1016/0010-4655(96)00104-X. Suche in Google Scholar

[4] Yin X, Zuo D. Soliton, lump and hybrid solutions of a generalized (2+1)-dimensional Benjamin-Ono equation in fluids. Nonlinear Dyn. 2024;1–19. 10.1007/s11071-024-10552-8. Suche in Google Scholar

[5] Zainab I, Akram G. Effect of β-derivative on time fractional Jaulent-Miodek system under modified auxiliary equation method and exp(-g(Ω))-expansion method. Chaos Solitons Fractals. 2023;168:113147. 10.1016/j.chaos.2023.113147. Suche in Google Scholar

[6] Yin Y-H, Lu X, Ma W-X. Bäcklund transformation, exact solutions and diverse interaction phenomena to a (3+1)-dimensional nonlinear evolution equation. Nonlinear Dyn. 2022;108:4181–94. 10.1007/s11071-021-06531-y. Suche in Google Scholar

[7] Chatziafratis A, Ozawa T, Tian SF. Rigorous analysis of the unified transform method and long-range instabilities for the inhomogeneous time-dependent Schrödinger equation on the quarter-plane. Math An. 2024;389:3535–602. 10.1007/s00208-023-02698-4. Suche in Google Scholar

[8] Li ZQ, Tian SF, Yang JJ. On the soliton resolution and the asymptotic stability of N-soliton solution for the Wadati-Konno-Ichikawa equation with finite density initial data in space-time solitonic regions. Adv Math. 2022;409:108639. 10.1016/j.aim.2022.108639. Suche in Google Scholar

[9] Li ZQ, Tian SF, Yang JJ. Soliton resolution for the Wadati-Konno-Ichikawa equation with weighted Sobolev initial data. Ann Henri Poincaré. 2022;23:2611–55. 10.1007/s00023-021-01143-z. Suche in Google Scholar

[10] Wu ZJ, Tian SF, Liu Y, Wang Z. Stability of smooth multisolitons for the two-component Camassa-Holm system. J London Math Soc. 2025;111:e70158. 10.1112/jlms.70158. Suche in Google Scholar

[11] Keller EF, Segel LA. Traveling bands of chemotactic bacteria: a theoretical analysis. J Theoret Biol. 1971;30(2):235–48. 10.1016/0022-5193(71)90051-8. Suche in Google Scholar PubMed

[12] Hillen T, Painter KJ. A user’s guide to PDE models for chemotaxis. J Math Biol. 2009;58(1):183–217. 10.1007/s00285-008-0201-3. Suche in Google Scholar PubMed

[13] Hillen T, Stevens A. Hyperbolic models for chemotaxis in 1-D. Nonlinear Anal Real World Appl. 2000;1(3):409–33. 10.1016/S0362-546X(99)00284-9. Suche in Google Scholar

[14] Rivero MA, Tranquillo RT, Buettner HM, Lauffenburger DA. Transport models for chemotactic cell populations based on individual cell behavior. Chem Eng Sci. 1989;44(12):2881–97. 10.1016/0009-2509(89)85098-5. Suche in Google Scholar

[15] Ford RM, Phillips BR, Quinn JA, Lauffenburger DA. Measurement of bacterial random motility and chemotaxis coefficients: I. Stopped-flow diffusion chamber assay. Biotech Bioeng. 1991;37(7):647–60. 10.1002/bit.260370707. Suche in Google Scholar PubMed

[16] Hillen T, Levine HA. Blow-up and pattern formation in hyperbolic models for chemotaxis in 1-D. Zeitschrift für angewandte Mathematik und Physik ZAMP. 2003;54:839–68. 10.1007/s00033-003-3206-1. Suche in Google Scholar

[17] Liu FG. Exact solutions to hyperbolic models for chemotaxis in 1-D. J Biomath. 2008;23:656–60. Suche in Google Scholar

[18] Boiti M, Leon JJ-P, Manna M, Pempinelli F. On the spectral transform of a Korteweg-de Vries equation in two spatial dimensions. Inverse Problems. 1986;2(3):271. 10.1088/0266-5611/2/3/005. Suche in Google Scholar

[19] Bai C-J, Zhao H. New solitary wave and Jacobi periodic wave excitations in (2+1)-dimensional Boiti–Leon–Manna–Pempinelli system. Int J Mod Phys B. 2008;22(15):2407–20. 10.1142/S021797920803954X. Suche in Google Scholar

[20] Darvishi MTA, Najafi Mb, Kavitha LC, Venkatesh MC. Stair and step soliton solutions of the integrable (2+1) and (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equations. Commun Theor Phys. 2012;58(6):785. 10.1088/0253-6102/58/6/01. Suche in Google Scholar

[21] Ma H, Bai Y, et al. Wronskian determinant solutions for the (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation. J Appl Math Phys. 2013;1(5):18–24. 10.4236/jamp.2013.15004. Suche in Google Scholar

[22] Mabrouk S-M, Rashed A-S. Analysis of (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation via Lax pair investigation and group transformation method. Comput Math Appl. 2017;74(10):2546–56. 10.1016/j.camwa.2017.07.033. Suche in Google Scholar

[23] Zuo DW, Gao YT, Yu X, Sun YH, Xue L. On a (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation. Z Naturforsch A. 2015;70(5):309–16. 10.1515/zna-2014-0340. Suche in Google Scholar

[24] Wu W, Manafian J, Ali KK, Karakoç SBG, Taqi AH, Mahmoud MA. Numerical and analytical results of the 1D BBM equation and 2D coupled BBM-system by finite element method. Int J Modern Phys B. 2022;36:2250201. 10.1142/S0217979222502010. Suche in Google Scholar

[25] Karakoç SBG, Ali KK. New exact solutionsand numerical approximations of the generalized kdv equation. Comput Methods Differ Equ. 2021;9:670–91. 10.22034/cmde.2020.36253.1628. Suche in Google Scholar

[26] Ali KK, Karakoç SBG, Rezazadeh H. Optical soliton solutions of the fractional perturbed nonlinear schrodinger equation. TWMS J Appl Eng Math. 2020;10:930–9. http://jaem.isikun.edu.tr/web/index.php/archive/108-vol10no4/604. Suche in Google Scholar

[27] Saha A, Karakoç BG, Ali KK. New exact soliton solutions, bifurcation and multistability behaviors of traveling waves for the (3+1)-dimensional modified Zakharov-Kuznetsov equation with higher order dispersion. Authorea. 2021. 10.22541/au.163533971.16788809/v1. Suche in Google Scholar

[28] Liu JG, Du JQ, Zeng ZF, Nie B. New three-wave solutions for the (3+ 1)-dimensional Boiti–Leon–Manna–Pempinelli equation. Nonlinear Dyn 2017;88:655–61. 10.1007/s11071-016-3267-2. Suche in Google Scholar

[29] Delisle L, Mosaddeghi M. Classical and SUSY solutions of the Boiti–Leon–Manna–Pempinelli equation. J Phys A Math Theor. 2013;46(11):115203. 10.1088/1751-8113/46/11/115203. Suche in Google Scholar

[30] Ablowitz MJ, Segur H. Solitons and the inverse scattering transform. Philadelphia: SIAM; 1981. 10.1137/1.9781611970883Suche in Google Scholar

[31] Hirota R, Satsuma J. N-soliton solutions of model equations for shallow water waves. J Phys Soc Japan. 1976;40(2):611–2. 10.1143/JPSJ.40.611. Suche in Google Scholar

[32] Wazwaz A-M. Painlevé analysis for new (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equations with constant and time-dependent coefficients. Int J Numer Methods Heat Fluid Flow. 2020;30(9):4259–66. 10.1108/HFF-10-2019-0760. Suche in Google Scholar

[33] Yuan N. Rich analytical solutions of a new (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation. Results Phys. 2021;22:10392. 10.1016/j.rinp.2021.103927. Suche in Google Scholar

[34] Liu JG, Du JQ, Zeng ZF, Nie B. New three-wave solutions for the (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation. Nonlinear Dyn. 2017;88:655–61. 10.1007/s11071-016-3267-2. Suche in Google Scholar

[35] Liu JG, Tian Y, Hu JG. New non-traveling wave solutions for the (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation. Appl Math Lett. 2018;79:162–8. 10.1016/j.aml.2017.12.011. Suche in Google Scholar

[36] Liu JG, Wazwaz AM. Breather wave and lump-type solutions of new (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation in incompressible fluid. Math Methods Appl Sci. 2021;44(2):2200–8. 10.1002/mma.6931. Suche in Google Scholar

[37] Rani M, Ahmed N, Dragomir SS, Mohyud-Din ST. Traveling wave solutions of 3+1-dimensional Boiti–Leon–Manna–Pempinelli equation by using improved tanh (ϕ 2)-expansion method. Partial Differ Equ Appl Math. 2022;6:100394. 10.1016/j.padiff.2022.100394. Suche in Google Scholar

[38] Akram G, Sadaf M, Zainab I. Observations of fractional effects of β-derivative and M-truncated derivative for space time fractional Phi-4 equation via two analytical techniques. Chaos Solitons Fractals. 2022;154:111645. 10.1016/j.chaos.2021.111645. Suche in Google Scholar

[39] Behera S, Aljahdaly NH, Virdi JPS. On the modified G′G2-expansion method for finding some analytical solutions of the traveling waves. J Ocean Eng Sci. 2022;7(4):313–20. 10.1016/j.joes.2021.08.013. Suche in Google Scholar

[40] Keerthana N, Saranya R, Annapoorani N. Dynamics and diffusion limit of traveling waves in a two-species chemotactic model with logarithmic sensitivity. Math Comput Simulat. 2024;222:311–29. 10.1016/j.matcom.2023.08.035. Suche in Google Scholar

Received: 2024-12-13

Revised: 2025-07-17

Accepted: 2025-08-05

Published Online: 2025-09-17

This work is licensed under the Creative Commons Attribution 4.0 International License.

Artikel in diesem Heft

https://doi.org/10.1515/phys-2025-0204

Schlagwörter für diesen Artikel

chemotaxis; partial differential equation; traveling wave solution; (3+1)-dimensional Boiti–Leon–Manna–Pempinelli equation

Creative Commons

BY 4.0