Exploring the interaction between lump, stripe and double-stripe, and periodic wave solutions of the Konopelchenko–Dubrovsky–Kaup–Kupershmidt system

1 Introduction

Studying nonlinear differential equations is crucial for understanding the complex dynamics governing physical phenomena in fluid dynamics, plasma, and elastic materials [1–7]. The analytical modeling of several linear and nonlinear, static, and dynamic elasticity models has become increasingly popular in recent years among scholars. Hence, several techniques have been developed for acquiring exact solutions of such models, including Lie group methods [8], Darboux and Bäcklund transformations [9], inverse scattering transform [10], Riemann–Hilbert method [11], the long wave limit method [12], Pfaffian technique [13], Painlevé analysis method [14], and Hirota’s bilinear method [15].

In this contribution, we aim to study the (2+1)-dimensional Konopelchenko–Dubrovsky–Kaup–Kupershmidt system (KDKKS) [16]

where φ = φ ( x , y , t ) , and the coefficients ℓ i ( i = 1, 2 , … , 8 ) are the real parameters.

Note that Eq. (1) can be used to describe various phenomena in various fields, including fluid mechanics, ocean dynamics, and plasma physics. It also covers several well-known equations as its special case. For example:

1. If we take ℓ 1 = 1 , ℓ 2 = − 6 , and ℓ 3 = ℓ 4 = ℓ 5 = ℓ 6 = ℓ 7 = ℓ 8 = 0 , Eq. (1) is reduced to the well-known Korteweg-de Vries equation

   φ t − 6 φ φ x + φ x x x = 0 ,

   which describes long waves in shallow water under gravity and ion-acoustic solitons in plasma [17].

2. If we take ℓ 1 = 1 , ℓ 2 = 6 , and ℓ 3 = ℓ 5 = ℓ 6 = ℓ 7 = ℓ 8 = 0 , Eq. (1) is reduced to a (2+1)-dimensional Kadomtsev–Petviashvili equation

   ( φ t + 6 φ φ x + φ x x x ) x + ℓ 4 φ y y = 0 ,

   which arises in the investigation of the weakly nonlinear dispersive waves in plasma and also in the modulation of the weakly nonlinear long waterwaves [18].

3. If we take ℓ 2 = 6 ℓ 1 , ℓ 6 = 4 ℓ 5 , and ℓ 3 = ℓ 4 = ℓ 7 = ℓ 8 = 0 , Eq. (1) is reduced to as a (2+1)-dimensional generalized breaking soliton equation or Bogoyavlensky–Konoplechenko (BK) model

   φ t + ℓ 1 φ x x x + 6 ℓ 1 φ φ x + ℓ 5 φ x x y + 4 ℓ 5 φ φ y + 4 ℓ 5 φ x ( ∂ x − 1 φ y ) = 0 ,

   for the interaction of a Riemann wave along the y axis and a long wave along the x axis [19,20]. Moreover, the symbol ∂ x − 1 ( . ) is used for the integral operator with respect to x .

4. If we take ℓ 1 = ℓ 2 = 0 , ℓ 3 = 1 , ℓ 4 = ℓ 5 = − 5 , and ℓ 6 = − 15 , ℓ 7 = 15 , ℓ 8 = 45 , Eq. (1) is reduced to a (2+1)-dimensional B-type Kadomtsev–Petviashvili equation

   φ t + φ x x x x x − 5 ( ∂ x − 1 φ y ) y − 5 φ x x y − 15 ( φ φ y + φ x ( ∂ x − 1 φ y ) ) + 15 ( φ x φ x x + φ φ x x x ) + 45 φ 2 φ x = 0 ,

   for the shallow water waves or electrostatic wave potential in plasma [21].

5. If we take ℓ 1 = ℓ 2 = 0 , ℓ 3 = − 1 , ℓ 4 = 5 , and ℓ 5 = ℓ 6 = ℓ 7 = ℓ 8 = − 5 , Eq. (1) is reduced to a (2+1)-dimensional Sawada–Kotera equation

   φ t − φ x x x x x + 5 ( ∂ x − 1 φ y y ) − 5 φ x x y − 5 ( φ φ y + φ x ( ∂ x − 1 φ y ) ) − 5 ( φ x φ x x + φ φ x x x ) − 5 φ 2 φ x = 0 ,

   in the atmosphere, rivers, lakes, and oceans, as well as the conformal field and two-dimensional quantum gravity gauge field [22].

6. Moreover, the Caudrey–Dodd–Gibbon–Sawada–Kotera equation [23] is also another special case of Eq. (1).

Due to its extensive use in real-world applications, the KDKK equation (Eq. (1)) has been the focus of many fascinating investigations in the literature. For instance, Cheng et al. [16] derived a bilinear form and proved Pfaffian solutions of up to N th order of Eq. (1) using the Hirota method, where N is the positive integer. These solutions reveal various soliton and breather structures, including elastic and inelastic interactions.

In the work presented in the study by Feng et al. [24], Bell’s polynomials and multidimensional Riemann theta functions are employed to suggest a methodical way to investigate periodic wave solutions and the asymptotic behavior of the system. Han and Bao [25] proposed two methods to find higher-order mixed localized wave solutions of the (3+1)-dimensional version of the system in Eq. (1). Furthermore, several numerical simulations and bilinear auto-Bäcklund transformations were used to obtain rich localized structures, including bell-shaped waves, periodic breather waves, and lump waves.

By reexamining the problem with a novel method based on bilinear transformation, a number of traveling wave solutions were obtained in the study by Jia [26]. These solutions include soliton molecules, solitons with double and triple-peaked structures, peak-plateau solitons, and few-cycle-pulse solitons. The focus of the research in the study by Ma et al. [27] is the (2+1)-dimensional generalized Konopelchenko–Dubrovsky–Kaup–Kupershmidt (gKDKK) equation. To achieve fission and fusion solutions, also known as “y-type solutions,” a new constraint condition is added based on Hirota’s bilinear method.

In the study by Li et al. [28], several soliton molecules were discovered by applying a velocity resonance ansatz to N -soliton solutions of the (2+1)-dimensional gKDKK equation. These solitons can transform into asymmetric solitons by adjusting specific parameters. Additionally, it was pointed out that a double-peaked lump solution to the system is constructed using a breather degeneration approach. Moreover, based on the Hirota’s bilinear form of this system, lump waves, lump-off waves, and rogue waves are described in the study by Liu et al. [29]. Deng et al. [30] employed the Pfaffian and Wronskian procedures to obtain the N th-order Pfaffian and Wronskian solutions, respectively. Higher-order lump, breather, and hybrid solutions to the system were also presented in the study of Zhang et al. [31]. As noted, the higher-order lump solutions are generalized through the long-wave limit method, while breather and hybrid solutions composed of solitons, breathers, and lumps are derived.

The (2+1)-dimensional gKDKK equation is solved in the study by Zhao et al. [32] using a bilinear form, yielding lump-type, interaction, and breather solutions. The trajectory, peaks, and valleys of the lump wave are analyzed theoretically and displayed visually through graphs. Feng et al. [33] explored a (3+1)-dimensional form of the Konopelchenko–Dubrovsky–Kaup–Kupershmidt equation using the Hirota method to derive bilinear forms, soliton solutions, breather solutions, and lump solutions. This equation has also been investigated in the study by Tariq et al. [34], where its integrability characteristics and wave patterns are explored. This study employs Hirota’s bilinear form and the exp ( − ϕ ( ξ ) ) expansion method to discover novel wave solutions for the model, including bell-shaped bright soliton solutions and periodic wave solutions.

In recent years, bilinear approaches have been widely utilized as essential tools for obtaining soliton wave solutions in higher-order nonlinear differential equations [20,35–37]. In this study, we aim to uncover novel traveling wave solutions for the ( 2 + 1 ) -dimensional KDKK system (1), using Hirota’s bilinear transformation and various combinations of interactions between lumps, stripes, and double stripes, as well as periodic waves. These solutions are new and innovative, having not been previously identified in the literature. To achieve our objectives, we have structured this article as follows: Section 2 presents the bilinear differential operator for the system. Then, we investigate novel wave solutions in Section 3 of our work. Section 3 additionally offers several types of 2D and 3D simulations that demonstrate the time evolution of the obtained solutions. The final section of the article summarizes the key conclusions.

2 Bilinear differential operator for Eq. (1)

To obtain the Hirota’s bilinear form of Eq. (1), we start by applying the following Cole–Hopf transformations:

where φ 0 and ϑ 0 are two constants.

Utilizing the transformations in Eq. (2) along with the following coefficient conditions [16]

which can be expressed in an equivalent form

the Hirota’s bilinear equation corresponding to Eq. (1) is derived as

In the Hirota’s bilinear form of (5), the operator of D x m D y n D t s is symbolized as

Using the bilinear differential operator specified in Eq. (6), the constituent equations of Eq. (5) are computed as

Taking (7) into account, we may reduce the proposed bilinear differential operator in (5) to a more simplified structure, as follows:

3 Interactions of wave solutions for system (6)

The combination of Hirota’s bilinear method and some interaction solutions has always been one of the most interesting and efficient approaches to determining analytical solutions for equations. In this section, we aim to discover several lump solutions of various interactions for Eq. (1). For this purpose, five different structures are considered for the solution of the equation.

3.1 Class 1: Collision among two lump solutions

In this section, we construct the following combination of two lump solutions:

(9) Λ = κ 0 + κ 1 Ξ 1 2 + κ 2 Ξ 2 2 ,

where

Ξ 1 = α 1 x + β 1 y + ρ 1 t , Ξ 2 = α 2 x + β 2 y + ρ 2 t ,

and the unknown constants α i , β i , and ρ i are yet to be determined.

The following solutions for the unknowns are obtained by using symbolic software and plugging Eq. (9) into Eq. (8):

Item 1-1:

β 2 = − α 2 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 , κ 1 = 0 , ρ 2 = α 2 ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) 4 ℓ 8 ,

and α 2 , κ 0 , and κ 2 are the free constants.

Substituting these results into (9) and then using transformation (2), we obtain

(10) Λ 1 = κ 0 + κ 2 t α 2 ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) 4 ℓ 8 + x α 2 − y α 2 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 2 .

Numerical results illustrating the solution, using (10) and (2), and taking ℓ 2 = 0.9 , ℓ 6 = 0.2 , ℓ 7 = 0.2 , ℓ 8 = 0.25 , φ 0 = ϑ 0 = 0.3 , along with κ 0 = 1.4 , κ 2 = 0.8 , α 2 = 1.5 , are depicted in Figure 1.

Figure 1

Time evolution of the obtained solution φ 1 ( x , y , t ) using (10) and (2).

Item 1-2:

β 1 = − α 1 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 , β 2 = − α 2 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 , ρ 1 = α 1 ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) 4 ℓ 8 , ρ 2 = α 2 ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) 4 ℓ 8 ,

and α 1 , α 2 , κ 0 , κ 1 , and κ 2 are the free constants.

Substituting these results in (9) and then using transformation (2), we obtain

(11) Λ 2 = κ 0 + κ 1 t α 1 ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) 4 ℓ 8 + x α 1 − y α 1 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 2 + κ 2 t α 2 ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) 4 ℓ 8 + x α 2 − y α 2 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 2 .

Numerical results illustrating the solution, using (11) and (2), and taking ℓ 2 = 0.2 , ℓ 6 = 0.1 , ℓ 7 = 0.2 , ℓ 8 = 0.05 , φ 0 = ϑ 0 = 0.3 , along with κ 0 = 6.5 , κ 1 = 1.8 , κ 2 = 2.6 , α 1 = 1.75 , and α 2 = 1.5 , are the depicted in Figure 2.

Figure 2

Time evolution of the obtained solution φ 2 ( x , y , t ) using (11) and (2).

Item 1-3:

β 1 = α 1 β 2 α 2 , κ 0 = 0 , κ 2 = − κ 1 α 1 2 α 2 2 , ρ 1 = − ( 2 ℓ 8 2 α 2 2 φ 0 2 + α 2 ( ( 2 ℓ 2 φ 0 + 2 ℓ 6 ϑ 0 ) α 2 + 2 β 2 ℓ 6 φ 0 + ρ 2 ) ℓ 8 − 2 ℓ 6 2 β 2 2 ) α 1 α 2 2 ℓ 8 ,

and κ 1 , α 1 , α 2 , β 2 , ρ 2 are the free constants.

Using these results in (9) and then utilizing transformation (2) yield

Λ 3 = κ 1 − t − 2 ℓ 6 2 β 2 2 + ( 2 β 2 ℓ 6 φ 0 + ρ 2 ) α 2 ℓ 8 ( 2 ℓ 2 φ 0 + 2 ℓ 6 ϑ 0 ) α 2 2 + 2 ℓ 8 2 α 2 2 φ 0 2 α 1 α 2 2 ℓ 8 + x α 1 + y α 1 β 2 α 2 2 − κ 1 α 1 2 ( t ρ 2 + x α 2 + y β 2 ) 2 α 2 2 .

Item 1-4:

κ 0 = 3 ℓ 7 κ 1 ℓ 8 φ 0 + ℓ 2 2 α 1 2 + β 1 ℓ 6 α 1 κ 1 + κ 2 ℓ 8 φ 0 + ℓ 2 2 α 2 + β 2 ℓ 6 α 2 ( κ 1 α 1 2 + κ 2 α 2 2 ) 2 ℓ 6 2 κ 1 κ 2 ( − α 1 β 2 + α 2 β 1 ) 2 , ρ 1 = − β 1 ℓ 6 α 2 κ 2 ( α 2 ℓ 8 φ 0 − 2 β 2 ℓ 6 ) − ℓ 8 κ 1 ( ℓ 8 φ 0 2 + ℓ 2 φ 0 + ℓ 6 ϑ 0 ) α 1 3 − ℓ 6 ℓ 8 α 1 2 β 1 κ 1 φ 0 + α 1 ( − ℓ 8 2 α 2 2 κ 2 φ 0 2 − κ 2 α 2 2 ( ℓ 2 φ 0 + ℓ 6 ϑ 0 ) ℓ 8 + ℓ 6 2 ( κ 1 β 1 2 − κ 2 β 2 2 ) ) ℓ 8 ( κ 1 α 1 2 + κ 2 α 2 2 ) , ρ 2 = β 2 ℓ 6 α 1 κ 1 ( − φ 0 ℓ 8 α 1 + 2 β 1 ℓ 6 ) − ℓ 8 κ 2 ( ℓ 8 φ 0 2 + ℓ 2 φ 0 + ℓ 6 ϑ 0 ) α 2 3 − ℓ 6 ℓ 8 α 2 2 β 2 κ 2 φ 0 + α 2 ( − ℓ 8 2 α 1 2 κ 1 φ 0 2 − α 1 2 κ 1 ( ℓ 2 φ 0 + ℓ 6 ϑ 0 ) ℓ 8 − ℓ 6 2 ( κ 1 β 1 2 − κ 2 β 2 2 ) ) ℓ 8 ( κ 1 α 1 2 + κ 2 α 2 2 ) ,

and α 1 , α 2 , β 1 , β 2 , κ 1 , and κ 2 are the free constants.

Substituting these results into (9) and then using transformation (2), we obtain

Λ 4 = 3 ℓ 7 κ 1 ℓ 8 φ 0 + ℓ 2 2 α 1 2 + β 1 ℓ 6 α 1 κ 1 + κ 2 ℓ 8 φ 0 + ℓ 2 2 α 2 + β 2 ℓ 6 α 2 ( κ 1 α 1 2 + κ 2 α 2 2 ) 2 ℓ 6 2 κ 1 κ 2 ( − α 1 β 2 + α 2 β 1 ) 2 + κ 1 t − β 1 ℓ 6 α 2 κ 2 ( α 2 ℓ 8 φ 0 − 2 β 2 ℓ 6 ) − ℓ 8 κ 1 ( ℓ 8 φ 0 2 + ℓ 2 φ 0 + ℓ 6 ϑ 0 ) α 1 3 − ℓ 6 ℓ 8 α 1 2 β 1 κ 1 φ 0 + α 1 ( − ℓ 8 2 α 2 2 κ 2 φ 0 2 − κ 2 α 2 2 ( ℓ 2 φ 0 + ℓ 6 ϑ 0 ) ℓ 8 + ℓ 6 2 ( κ 1 β 1 2 − κ 2 β 2 2 ) ) ℓ 8 ( κ 1 α 1 2 + κ 2 α 2 2 ) + x α 1 + y β 1 2

+ κ 2 t β 2 ℓ 6 α 1 κ 1 ( − φ 0 ℓ 8 α 1 + 2 β 1 ℓ 6 ) − ℓ 8 κ 2 ( ℓ 8 φ 0 2 + ℓ 2 φ 0 + ℓ 6 ϑ 0 ) α 2 3 − ℓ 6 ℓ 8 α 2 2 β 2 κ 2 φ 0 + α 2 ( − ℓ 8 2 α 1 2 κ 1 φ 0 2 − α 1 2 κ 1 ( ℓ 2 φ 0 + ℓ 6 ϑ 0 ) ℓ 8 − ℓ 6 2 ( κ 1 β 1 2 − κ 2 β 2 2 ) ) ℓ 8 ( κ 1 α 1 2 + κ 2 α 2 2 ) + x α 2 + y β 2 2 .

It is noteworthy to note that although the structure (9) employed in our study bears some parallels to that of the publications [32–34], our obtained results are more thorough and in general distinct from those reported in these articles.

3.2 Class 2: Collision among two lump solutions and stripe soliton

In this section, we construct a test function of

(12) Λ = κ 0 + Ξ 1 2 + κ 1 Ξ 2 2 + κ 2 exp ( Ξ 3 ) ,

where

Ξ 1 = α 1 x + β 1 y + ρ 1 t , Ξ 2 = α 2 x + β 2 y + ρ 2 t , Ξ 3 = α 3 x + β 3 y + ρ 3 t .

The following solutions for the unknowns are obtained by using symbolic software and plugging Eq. (12) into Eq. (8):

Item 2-1:

β 1 = − α 1 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 , β 3 = − α 3 ( 2 ℓ 7 α 3 2 + 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 , κ 1 = 0 , ρ 1 = α 1 ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) 4 ℓ 8 , ρ 3 = 4 ℓ 8 2 φ 0 2 + 2 ( 6 ℓ 7 α 3 2 φ 0 + ℓ 2 φ 0 − 2 ℓ 6 ϑ 0 ) ℓ 8 + 36 ℓ 7 2 α 3 4 5 + 4 ℓ 2 ℓ 7 α 3 2 + ℓ 2 2 α 3 4 ℓ 8

and α 1 , α 3 , κ 0 , and κ 2 are the free constants.

By considering these results in (12) and subsequently applying transformation (2), we obtain

(13) Λ 5 = κ 0 + t α 1 ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) 4 ℓ 8 + x α 1 − y α 1 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 2 + κ 2 exp t 2 ( 6 ℓ 7 α 3 2 φ 0 + ℓ 2 φ 0 − 2 ℓ 6 ϑ 0 ) ℓ 8 + 4 ℓ 8 2 φ 0 2 + 36 ℓ 7 2 α 3 4 5 + 4 ℓ 2 ℓ 7 α 3 2 + ℓ 2 2 α 3 4 ℓ 8 + x α 3 − y α 3 ( 2 ℓ 7 α 3 2 + 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 .

Numerical results illustrating the solution, using (13) and (2), and taking ℓ 2 = 1.4 , ℓ 6 = 2.8 , ℓ 7 = 0.9 , ℓ 8 = 0.3 , φ 0 = ϑ 0 = 0.3 , along with κ 0 = 0.3 , κ 2 = 1.2 , α 1 = 0.35 , α 3 = 0.22 , are depicted in Figure 3.

Figure 3

Time evolution of the obtained solution φ 5 ( x , y , t ) using (13) and (2).

Item 2-2:

β 1 = − α 1 ( 3 ℓ 2 ℓ 7 α 2 α 3 4 + 2 ℓ 7 2 α 2 α 3 6 + 2 ℓ 6 2 β 3 ( − α 3 β 2 + α 2 β 3 ) + ℓ 7 α 3 3 ( ℓ 6 ( α 3 β 2 + 4 α 2 β 3 ) + 6 ℓ 8 α 2 α 3 φ 0 ) ) ℓ 6 α 3 ( ℓ 7 α 2 α 3 3 + 2 ℓ 6 ( α 3 β 2 − α 2 β 3 ) ) ,

κ 0 = 16 α 1 2 ℓ 6 β 3 + 1 2 α 3 ( ℓ 2 − ℓ 7 α 3 2 + 2 ℓ 8 φ 0 ) ℓ 7 2 α 2 2 α 3 6 + ℓ 6 2 ( α 3 β 2 − α 2 β 3 ) 2 + 2 ℓ 7 α 2 2 α 3 3 ℓ 6 β 3 + ℓ 7 α 2 α 3 4 ℓ 6 β 2 + 3 2 α 2 ( ℓ 2 + 2 ℓ 8 φ 0 ) ( ℓ 7 α 2 α 3 4 + 2 ℓ 6 α 3 ( α 3 β 2 − α 2 β 3 ) ) 2 ( 4 ℓ 6 β 3 + α 3 ( 2 ℓ 2 + ℓ 7 α 3 2 + 4 ℓ 8 φ 0 ) ) ,

κ 1 = 3 ( ℓ 7 α 3 3 + 4 ℓ 8 α 3 φ 0 + 2 ℓ 2 α 3 + 4 ℓ 6 β 3 ) ℓ 7 α 1 2 α 3 3 ( − ℓ 7 α 2 α 3 3 + 2 ℓ 6 α 2 β 3 − 2 ℓ 6 α 3 β 2 ) 2 ,

ρ 1 = − α 1 ℓ 6 ℓ 7 α 3 5 ℓ 6 β 2 ( 3 ℓ 2 + 5 ℓ 8 φ 0 ) + ℓ 7 α 3 2 ℓ 6 β 2 − 5 α 2 3 ℓ 2 10 + ℓ 8 φ 0 − ℓ 7 2 α 2 α 3 4 + ℓ 8 α 2 ( ℓ 6 ϑ 0 − φ 0 ( 2 ℓ 2 + 5 ℓ 8 φ 0 ) ) × ( − 2 ℓ 6 α 3 β 3 2 ( ℓ 6 β 2 + ℓ 8 α 2 φ 0 ) + ℓ 7 α 3 4 β 3 ( 8 ℓ 6 β 2 + α 2 ( 3 ℓ 2 + 2 ℓ 8 φ 0 ) ) + 2 ℓ 6 2 α 2 β 3 3 + α 3 2 ( 2 ℓ 8 ( α 3 β 2 − α 2 β 3 ) φ 0 ( ℓ 2 + ℓ 8 φ 0 ) + ℓ 6 ( 3 ℓ 7 α 2 α 3 β 3 2 + 2 ℓ 8 ( ( α 3 β 2 − α 2 β 3 ) ϑ 0 + β 2 β 3 φ 0 ) ) ) ) 4 α 3 2 1 2 ℓ 7 α 2 α 3 3 + ℓ 6 α 3 β 2 − ℓ 6 α 2 β 3 ℓ 8 ,

ρ 2 = − ℓ 7 2 α 2 α 3 6 + 2 ℓ 6 ℓ 7 α 2 α 3 3 β 3 + ℓ 6 2 β 3 ( − 2 α 3 β 2 + α 2 β 3 ) + ℓ 7 α 3 4 ℓ 6 β 2 + 3 2 α 2 ( ℓ 2 + 2 ℓ 8 φ 0 ) + ℓ 8 α 3 2 ( α 2 φ 0 ( ℓ 2 + ℓ 8 φ 0 ) + ℓ 6 ( α 2 ϑ 0 + β 2 φ 0 ) ) ℓ 8 α 3 2 ,

ρ 3 = − 2 ℓ 7 2 α 3 6 − 5 ℓ 7 ( 2 ℓ 8 φ 0 + ℓ 2 ) α 3 4 − 10 ℓ 6 ℓ 7 α 3 3 β 3 − 10 ℓ 8 ( ℓ 8 φ 0 2 + ℓ 2 φ 0 + ℓ 6 ϑ 0 ) α 3 2 − 10 ℓ 6 ℓ 8 α 3 β 3 φ 0 + 10 ℓ 6 2 β 3 2 10 ℓ 8 α 3 ,

and α 1 , α 2 , α 3 , β 2 , β 3 , and κ 2 are the free constants.

Integrating these results in (12) and then applying transformation (2) yield the analytic solution Λ 6 to the model.

Item 2-3:

β 1 = − α 1 ( 2 ℓ 2 + 3 ℓ 7 α 3 2 + 4 ℓ 8 φ 0 ) 4 ℓ 6 , β 2 = − α 2 ( 2 ℓ 2 + 3 ℓ 7 α 3 2 + 4 ℓ 8 φ 0 ) 4 ℓ 6 , β 3 = − α 3 ( α 3 2 ℓ 7 + 4 φ 0 ℓ 8 + 2 ℓ 2 ) 4 ℓ 6 , ρ 1 = α 1 ℓ 8 ( φ 0 ( 9 α 3 2 ℓ 7 + 2 ℓ 2 ) − 4 ℓ 6 ϑ 0 ) + 3 2 α 3 2 ℓ 7 + ℓ 2 2 + 4 φ 0 2 ℓ 8 2 4 ℓ 8 , ρ 2 = α 2 ℓ 8 ( φ 0 ( 9 α 3 2 ℓ 7 + 2 ℓ 2 ) − 4 ℓ 6 ϑ 0 ) + 3 2 α 3 2 ℓ 7 + ℓ 2 2 + 4 φ 0 2 ℓ 8 2 4 ℓ 8 , ρ 3 = 4 ℓ 8 2 φ 0 2 + ( 3 ℓ 7 α 3 2 φ 0 + 2 ℓ 2 φ 0 − 4 ℓ 6 ϑ 0 ) ℓ 8 + 9 ℓ 7 2 α 3 4 20 + ℓ 2 ℓ 7 α 3 2 + ℓ 2 2 α 3 4 ℓ 8 , κ 1 = − α 1 2 α 2 2 ,

and α 1 , α 2 , α 3 , κ 0 , and κ 2 are the free constants.

Taking these results in (12) and then using transformation (2), we obtain

Λ 7 = κ 0 + t 4 ℓ 8 2 φ 0 2 + 3 ℓ 7 α 3 2 2 + ℓ 2 2 + ( 9 ℓ 7 α 3 2 φ 0 + 2 ℓ 2 φ 0 − 4 ℓ 6 ϑ 0 ) ℓ 8 α 1 4 ℓ 8 + x α 1 − y 3 ℓ 7 α 3 2 2 + 2 ℓ 8 φ 0 + ℓ 2 α 1 2 ℓ 6 2 − α 1 2 α 2 2 t 4 ℓ 8 2 φ 0 2 + 3 ℓ 7 α 3 2 2 + ℓ 2 2 + ( 9 ℓ 7 α 3 2 φ 0 + 2 ℓ 2 φ 0 − 4 ℓ 6 ϑ 0 ) ℓ 8 α 2 4 ℓ 8 + x α 2 − y 3 ℓ 7 α 3 2 2 + 2 ℓ 8 φ 0 + ℓ 2 α 2 2 ℓ 6 2 + κ 2 exp t 4 ℓ 8 2 φ 0 2 + 9 ℓ 7 2 α 3 4 20 + ℓ 2 ℓ 7 α 3 2 + ℓ 2 2 ( 3 ℓ 7 α 3 2 φ 0 + 2 ℓ 2 φ 0 − 4 ℓ 6 ϑ 0 ) ℓ 8 α 3 4 ℓ 8 + x α 3 − y α 3 ℓ 7 α 3 2 2 + 2 ℓ 8 φ 0 + ℓ 2 2 ℓ 6 .

Item 2-4:

β 1 = − α 1 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 , β 2 = − α 2 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 , β 3 = − α 3 ( 2 ℓ 7 α 3 2 + 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 , ρ 1 = α 1 ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) 4 ℓ 8 , ρ 2 = ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) α 2 4 ℓ 8 , ρ 3 = 4 ℓ 8 2 φ 0 2 + 2 ( 6 ℓ 7 α 3 2 φ 0 + ℓ 2 φ 0 − 2 ℓ 6 ϑ 0 ) ℓ 8 + 36 ℓ 7 2 α 3 4 5 + 4 ℓ 2 ℓ 7 α 3 2 + ℓ 2 2 α 3 4 ℓ 8 ,

and α 1 , α 2 , α 3 , κ 0 , κ 1 , and κ 2 are the free constants.

Taking these results in (12) and then using transformation (2), we obtain

(14) Λ 8 = κ 0 + t α 1 ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) 4 ℓ 8 + x α 1 − y α 1 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 2 + κ 1 t ( 4 ℓ 8 2 φ 0 2 + 2 ℓ 2 ℓ 8 φ 0 − 4 ℓ 6 ℓ 8 ϑ 0 + ℓ 2 2 ) α 2 4 ℓ 8 + x α 2 − y α 2 ( 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6 2 + κ 2 exp t 2 ( 6 ℓ 7 α 3 2 φ 0 + ℓ 2 φ 0 − 2 ℓ 6 ϑ 0 ) ℓ 8 + 4 ℓ 8 2 φ 0 2 + 36 ℓ 7 2 α 3 4 5 + 4 ℓ 2 ℓ 7 α 3 2 + ℓ 2 2 α 3 4 ℓ 8 + x α 3 − y α 3 ( 2 ℓ 7 α 3 2 + 2 ℓ 8 φ 0 + ℓ 2 ) 2 ℓ 6

Numerical results illustrating the solution, using (14) and (2), and taking ℓ 2 = 0.4 , ℓ 6 = 0.8 , ℓ 7 = 0.9 , ℓ 8 = 0.17 , φ 0 = ϑ 0 = 0.3 , along with κ 0 = 0.5 , κ 1 = 0.2 , κ 2 = 0.6 , α 1 = 0.75 , α 2 = 2.2 , and α 3 = 0.5 , are the depicted in Figure 4.

Figure 4

Time evolution of the obtained solution φ 8 ( x , y , t ) using (14) and (2).