Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) Symbolic Computational Algorithm for Hirota Bilinear Form to Higher-dimensional Nonlinear Partial Differential Equations in Nonlinear Sciences Ankit Dadhich1, Hari Pratap2, Brij Mohan3 1Department of Mathematics, Central University of Haryana, Haryana-123031, India 2*Department of Mathematics, PGDAV(E) College, University of Delhi, Delhi 110065, India 3Department of Mathematics, Hansraj College, University of Delhi, Delhi -110007, India Corresponding Author: pratap.pgdavcollege.hari@gmail.com Received: 01-12-2024 Revised: 06-12-2024 Accepted: 12-12-2024 Published: 20-12-2024 Abstract: In the physical world, many real systems are governed by nonlinear partial differential equations from fluid dy- namics and plasma phyasics to shallow-water waves and oceanographic systems. There is no uniform approach for solving nonlinear partial differential equations; consequently, we consider each equation as a separate problem. The most effective technique for building multi-soliton solutions of an integrable nonlinear PDE is the Hirota direct method from the several methods used to explore nonlinear PDEs to obtain solitons, lumps, rogue waves, breathers, and kink waves. To derive the multi-solitons of the non-linear PDEs, this method requires first transforming the equation into the bilinear form proposed by Hirota, and then applying the dependent variable transformation. Converting nonlinear evolution equations into Hirota bilinear form is the primary goal of the research study. This work investigated several well-known nonlinear PDEs, including the KdV Equation, Boussinesq Equation, GS Equation, KP equation, and other equations in order to comprehend and apply the concept. We design the algorithm using symbolic software Mathematica. These equations are widely recognized for multi-solitons and their integrability, which having applications in diverse fileds oceanography, fluid dynamics, plasma physics, mechanics, and other nonlinear sciences. Keyword: Symbolic Computation, Bilinear Form, Nonlin- ear PDE, Symbolic Algorithm AMS Subject Classification: 35G20, 47A07, 68W30, 1 Introduction In mathematics and science, a partial differential equation that has nonlinear components is called a nonlinear evolution equation or partial differential equation (PDE). They describe distinct physical systems from water dynamics to gravity, and have been used in mathematics to solve important problems such as the Poincaré and Calabi conjectures. They are challenging to analyze since there aren’t many general techniques that work for all of these equations, and usually, each one needs to be looked at separately. The Hirota method, https://internationalpubls.com 1212 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) which Ryogo Hirota introduced in 1971 [1–5], is a popular and powerful mathematical tool for identifying soliton solutions after changing the PDE to bilinear form. In several disciplines, such as oceanography, fluid dynamics, optics, plasma physics, and engineering sciences, soliton is a highly regular solution that has wave-like properties. Therefore bilinearization of nonlinear PDEs is the most necessary step. In physics and mathematics, a soliton is a confined wave packet that is extremely stable, nonlinear, and self-reinforcing. After colliding with other localized wave packets, it maintains its form even when moving freely at a constant speed. Its exceptional stability results from the medium’s dispersive and nonlinear effects being balancedly cancelled. Weakly nonlinear dispersive PDEs have wide group to describe nonlinear systems was eventually demonstrated to have stable solutions in soliton theory [6]. For integrable non-linear PDEs, precise N-soliton or multisoliton solutions may be found using Hirota’s bilinear approach. However, Hirota asserts that converting a non-linear PDE into its bilinear form is the most crucial step in this process. Even when the correct transformation of the dependant variable is known, the process of converting a nonlinear evolution equation into a bilinear equation becomes time-consuming. As a result, creating a method to determine a nonlinear PDE’s bilinear form is crucial. To do these calculations, system software like Maple Matlab, and Mathematica, are useful. Hirota first applied this method to the KdV equation χt + 6χχx + χxxx = 0, converting it into the bilinear form (D4 x +DxDt)h · h = 0, by using the transformation χ = 2(lnh)xx, which resulted in a simple soliton solution development [3]. Many scholars have been interested in the nonlinearity of PDEs and have used a number of systematic techniques to obtain multiple lump, breather, and soliton solutions. There are several different methods that are to be used such as Bäcklund transformation [7,8], Darboux transformation [9–12], Lie symmetry analysis [13,14], simplified Hirota method [15,16], Pfaffian technique [17,18], and other methods. In this work, we investigate how several nonlinear PDEs, including the Korteweg-de Vries (KdV) [4], Kadomtsev–Petviashvili (KP) equations [5], and other equations that are transformed into Hirota bilinear form. The most effective technique for building multi-soliton solutions of an integrable nonlinear PDE is the Hirota direct approach. Solitons are crucial to the analysis of shallow water waves because they are created when the nonlinearity and dispersion effect are ignored. Moreover, they are found in a number of disciplines, including fluid dynamics, dusty plasma, oceanography, marine engineering, and plasma physics. The following section 2 explores the general structure of the bilinear form in D-operator for a general nonlinear PDE. It describe the Cole-Hopf transformation, Hirota’s D-operator and algorithm for bilinear form. In section 3, we show the application of algorithm to several well-known equations to obtain the bilinear forms. Section 4 discusses the results of the bilinear form for the studied equations, and the ending section concludes the research work. 2 Hirota Bilinear Form of a nonlinear PDE The algorithm for obtaining the Hirota Bilinear Form and its application to differential equation solving will be the main topics of the next part. Hirota provided an algebraic method for finding precise soliton solutions, provided that the nonlinear PDE could be transformed into bilinear form. Initially, dependent variable transformation was used to convert the PDE to a bilinear equation for an auxiliary function [1]. As we will see, the concept is based on the properties of a certain bilinear differential operator known as the Hirota bilinear operator (D-operator) [19]. https://internationalpubls.com 1213 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 2.1 Cole–Hopf transformations A mathematical method for investigating partial differential equations (PDEs), especially nonlinear PDEs, is the Cole–Hopf transformation. Cole and Hopf [20,21] created the transformation in the 1950s to simplify and sometimes linearize certain types of nonlinear PDEs.The Cole–Hopf transformation has been found to be a great mathematical tool to explore solitons and the integrable equations. It helps scientists to investigate the wave structures of nonlinear equations to explore the existence of solitary waves, lump solutions, and several others that can persist in specific nonlinear systems. We have Cole-Hopf transformation u = K (log h)xp , where p is the order of partial derivative in x based on the balance between the PDE’s higher-order and nonlinear terms. The phase variable must be used to obtain the dispersion in order to create the mentioned transformation [22], which we’ll discuss through examples in a further section. 2.2 The Hirota Bilinear Operator Hirota gave the D-operator, a binary form that outputs a new function after receiving two functions as input. It has numerous qualities that make it helpful for differential equation analysis. In particular, it enables us to identify soliton solutions, which are analytic solutions to them. A pair of functions provide the operator’s definition as (h, k) of a real variable a [19], is Da(h, k) = lim a′→a ( ∂ ∂a − ∂ ∂a′ ) h(a)k(a′). (1) When the operator is applied repeatedly and to distinct variables (a and b in this case), but it is obviously generalizable to any number of real variables, we expand the definition [4] to Dp aD q b(h, k) = lim a′→a, b′→b ( ∂ ∂a − ∂ ∂a′ )p( ∂ ∂b − ∂ ∂b′ )q h(a, b)k(a′, b′). Below are some brief outcomes to help the D-operator get some intuition. Dy(h · k) = hyk − hky, D4 x(h · k) = hxxxxk − 4.hxxxkx + 6.hxxkxx − 4.hxkxxx + hkxxxx, D2 y(h.k) = hyyk + hkyy + 2hyky Actually, we may observe the following formula, which is created by applying binomial expansion to powers of (1), for the expansion of the nth power of the operator., Dn xf1 · f2 = n∑ k=0 (−1)k ( n k ) ∂n−kx f1(x) ∂kxf2(x). (2) 2.3 Algorithm for bilinear Form We begin with a nonlinear PDE of the form in (n+ 1) dimensions as P (u, ux1 , ux2 , . . .) = 0, (3) where the function u = u(x1, x2, . . . , xn, t) depends on the spatial variables x1, x2, . . . , xn and the temporal variable t, and the operator P involves u together with its partial derivatives. https://internationalpubls.com 1214 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) Step 1: Analysis of the dispersion relation We introduce the phase variable ξi as ξi = k1ix1 + k2ix2 + k3ix3 + · · · + knixn + ωit, (4) where ωi represents the dispersion relation and the coefficients kNi (1 ≤ N ≤ n) are constant wave numbers. The above form of the phase variable is a standard choice, illustrating the adaptability of the method to different classes of nonlinear PDEs. Depending on the structure of the PDE, however, this form may vary. We substitute the trial solution u = eξi in linear terms of Eq. (3) and solve for dispersion ωi. Step 2: Determining the transformation constant R We now apply the logarithmic transformation u = K ∂η ∂χη (lnh), (5) where η is chosen so that the order of the highest derivative in Eq. (3) balances with the nonlinear terms. For testing purposes, we consider h = 1 + eξ1 . Substituting this into Eq. (5) provides different possible values of the parameter K, from which a suitable one can be selected. Step 3: Deriving a quadratic or quartic equation in Φ Substituting the transformation (5) into the original PDE (3) and integrating in x (at least once) leads to an algebraic relation in terms of h. By setting the integration constant to zero at the lowest admissible order, we obtain either a quadratic or a quartic form, Q(h, hx1 , hx2 , . . .) = 0, (6) where Q involves h together with its derivatives with respect to t and the spatial variables x1, x2, . . . , xn. Step 4: Reformulation into Hirota’s bilinear form The quadratic structure obtained in Step 3 can be expressed in terms of Hirota’s bilinear operators. The Hirota derivative is defined by Dp aD q b(h, k) = lim a′→a, b′→b ( ∂ ∂a − ∂ ∂a′ )p( ∂ ∂b − ∂ ∂b′ )q h(a, b)k(a′, b′). Using this representation, Eq. (6) can be rewritten in bilinear form, thereby yielding the Hirota bilinear expression corresponding to the original nonlinear PDE (3). 3 Applications of Hirota Bilinear form The process for obtaining Hirota bilinear forms is often illustrated using classical examples of integrable equations, such as the KdV and the KP equation. These equations are essential to the theory of solitons and provide excellent examples of how well Hirota’s approach works. 3.1 KdV equation in (1+1)-dimensions The nonlinear KdV equation models the evolution of long, weak nonlinear waves in one spatial dimension [22]. Originally, it was described in hydro-dynamics to discuss the wave motion in shallow water. This equation is ut + 6uux + uxxx = 0, (7) https://internationalpubls.com 1215 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) where x denotes the spatial coordinate, t represents time, and u as dependent variable for wave amplitude. A remarkable property of Eq. (7) is the existence of solitons. The persistence of solitons is a direct con- sequence of the balance between the nonlinear and dispersive contributions in the equation. Because of this, the KdV equation has become a central model in the study of nonlinear wave dynamics, with applications ranging across oceanography, plasma physics, and optics. To explore the dispersion relation in the KdV equation (7), we introduce the phase variable ξi = µix+ dit, where µi (i = 1, 2, . . .) are constants and di is the dispersion parameter. Substituting the trial solution u = eξi into the linear part of Eq. (7) gives u = eξi , ut = die ξi , uxxx = µ3i e ξi . Thus, from the relation ut + uxxx = 0, we obtain die ξi + µ3i e ξi = eξi ( µ3i + di ) = 0, which implies di = −µ3i . Hence, the phase variable becomes ξi = µix− µ3i t. Next, we construct solutions of Eq. (7) using the Cole–Hopf transformation, which introduces a logarithmic representation of the dependent variable: u = K(log h)xx, (8) where h = h(x, t). It may be written as u = wxx, with w = K(log h). (9) To determine the value of K, we choose the auxiliary function in the Cole–Hopf transformation as h(x, t) = 1 + eξ1 = 1 + eµ1x+d1t = 1 + eµ1x−µ3 1t. (10) Substituting Eq. (10) into Eq. (7) and get the K as K = 2. Therefore, the logarithmic transformation reduces to u = 2(logK)xx. (11) Now, from (9), we have ut = wxxt, ux = wxxx and uxxx = wxxxxx putting the above expressions into Eq.(7), we get wxxt + 6wxxwxxx + wxxxxx (12)= 0 https://internationalpubls.com 1216 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) on integrating w.r.t. x wxt + 6 ∫ wxxwxxx ∂x+ wxxxx = 0. (13) We compute integral term in equation (13) with constant of integration as zero I = ∫ wxxwxxx ∂x = 1 2 ∫ 2wxxwxxx ∂x = 1 2 w2 xx, (14) substituting the value of I in equation (13), we get wxt + 3w2 xx + wxxxx = 0. (15) As we have w = 2(log f), we can get the followings: wx = 2 hx h , wxt = 2 hhxt − hxht h2 , wxx = 2 hhxx − h2x h2 , wxxx = 2 2h3x − 3hhxhxx + h2hxxx h3 , wxxxx = 2 −6h4x + 12hh2xhxx − 3h2h2xx − 4hxh 2hxxx h4 , putting all the above values in equation (15), we get a quadratic equation in h as −2 hxht h2 + 2 hxt h + 6 h2xx h2 − 8 hxhxxx h2 + 2 hxxxx h = 0, or hhxt − hxht + 3h2xx − 4hxhxxx + hhxxxx = 0, (16) Since Hirota bilinear operator is defined as Dp aD q b(h, k) = lim a′→a, b′→b ( ∂ ∂a − ∂ ∂a′ )p( ∂ ∂b − ∂ ∂b′ )q h(a, b)k(a′, b′). let p = 1 and q then= 1, DxDt(h.k) = hxtk − htkx − hxkt + hkxt, DxDt(h.h) = 2 (hhxt − hxht) let p = 4 and q then= 0, D4 x(h.k) = hxxxxk − 4.hxxxkx + 6.hxxkxx − 4.hxkxxx + hkxxxx D4 x(h.h) = 2 ( 3h2xx − 4hxhxxx + hhxxxx ) Equation (16) can be rewritten in terms of the Hirota D operator as (DxDt +D4 x)h.h (17)= 0 This equation (17) is called the Hirota bilinear form for KDV equation (7). https://internationalpubls.com 1217 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 3.2 Boussinesq Equation in (1+1)-dimensions The integrable Boussinesq equation [23] is structured as utt − 3(u2)xx − uxx − uxxx = 0. (18) To investigate the dispersion relation for Eq. (18), we introduce the phase variable ξi = µix+ dit, (19) where µi (i = 1, 2, . . .) are constants and di is the dispersion parameter. Substituting the trial solution u = eξi into the linear part of Eq. (18), namely utt − uxx − uxxx = 0, leads to the relation di = √ µ2i + µ4i . Thus, the phase variable becomes ξi = µix+ √ µ4i + µ2i t. Next, we apply a logarithmic-type transformation of the dependent variable: u = K(lnh)xx, (20) which may be equivalently expressed as u = wxx, where w = K(lnh). (21) Choosing the function h as h = 1 + eξ1 , with ξ1 = µ1x+ √ µ21 + µ41 t, and substituting into Eq. (18), we obtain K = 2. Hence, the transformation (20) simplifies to u(x, t) = 2(lnh)xx. Now, from (21), we have utt = wxxtt, uxx = wxxxx and uxxx = wxxxxx putting the above expressions into Eq.(18), we get wxxtt + wxxxx − 3 (w2 xx)xx − wxxxxx (22)= 0 on integrating w.r.t. x wxtt + wxxx − 3 (w2 xx)x − wxxxx (23)= 0 on again integrating w.r.t. x wtt + wxx − 3 (w2 xx) − wxxx (24)= 0 As we have w = 2(log h), we can get the followings: wx = 2 hx h , https://internationalpubls.com 1218 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) wxx = 2 hhxx − h2x h2 , wxxx = 2 2h3x − 3hhxhxx + h2hxxx h3 , wtt = 2 hhtt − h2t h2 , which changes Eq.(18) into a bilinear equation in h as hhtt − h2t − hhxx + h2x − hhxxxx + 4hxhxxx − 3h2xx (25)= 0 Since Hirota bilinear operator is defined as Dp aD q b(h, k) = lim a′→a, b′→b ( ∂ ∂a − ∂ ∂a′ )p( ∂ ∂b − ∂ ∂b′ )q h(a, b)k(a′, b′). let p = 0 and q then= 2, D2 t (h.k) = httk − 2htkt + hktt, D2 t (h.h) = 2 ( hhtt − h2t ) let p = 2 and q then= 0, D2 x(h.k) = hxxk − 2hxkx + hkxx, D2 x(h.h) = 2 ( hhxx − h2x ) let p = 4 and q then= 0, D4 x(h.k) = hxxxxk − 4.hxxxkx + 6.hxxkxx − 4.hxkxxx + hkxxxx D4 x(h.h) = 2 ( 3.h2xx − 4.hxhxxx + hhxxxx ) Equation (25) can be written in terms of the operator D as (D2 t −D2 x −D4 x)h.h (26)= 0 This equation (26) is called the Hirota bilinear form for Boussinesq equation (18). 3.3 KP equation in (2+1)-dimensions The KP equation extends the KdV model, which was originally formulated in two dimensions [24]. Its mathematical form is (ut + 6uux + uxxx)x − uyy = 0, (27) where x and y and t are spatial and temporal variable, and u represents the wave amplitude. The KP equation is a fundamental model in the theory of solitons and integrable systems, as it admits solu- tions that describe nonlinear wave interactions in two dimensions. Similar to the KdV equation, it supports soliton solutions; however, owing to its higher dimensionality, it captures more intricate behaviors, such as two-soliton interactions with nontrivial structures. Applications of the KP equation arise in diverse areas, including oceanography, plasma physics, and fluid mechanics, where it provides insight into multidimensional nonlinear wave dynamics. For the KP equation (27), we define the phase variable as ξi = µix+ viy − dit, (28) https://internationalpubls.com 1219 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) where µi, vi (i = 1, 2, . . .) are constants, and di denotes the dispersion parameter. Substituting the trial solution u = eξi in linear part of Eq. (27), we obtain the dispersion as di = µ4i − v2i µi . (29) Next, we introduce the transformation u(x, y, t) = K(lnh)xx, (30) which may also be written as u = wxx, with w = K(lnh). (31) Choosing the function h as h = 1 + eξ1 , where ξ1 = µ1x+ v1y − ( µ41 − v21 µ1 ) t, and substituting into Eq. (27), the value of K is determined to be K = 2. Thus, the logarithmic transformation (30) reduces to u = 2(lnh)xx. (32) Now, from (31), we have ut = wxxt, ux = wxxx, uxxx = wxxxxx and uyy = wxxyy putting the above expressions into Eq.(27), we get (wxxt + 6wxxwxxx + wxxxxx)x − wxxyy (33)= 0 on integrating w.r.t. x wxxt + 6wxxwxxx + wxxxxx − wxyy = 0. (34) on again integrating w.r.t. x wxt + 6 ∫ wxxwxxx ∂x+ wxxxx − wyy = 0. (35) We compute integral term in equation (35) with constant of integration as zero I = ∫ wxxwxxx ∂x = 1 2 ∫ 2wxxwxxx ∂x = 1 2 w2 xx, (36) substituting the value of I in equation (35), we get wxt + 3w2 xx + wxxxx − wyy = 0. (37) As we have w = 2(log h), we can get the followings: wx = 2 hx h , wxt = 2 hhxt − hxht h2 , https://internationalpubls.com 1220 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) wxx = 2 hhxx − h2x h2 , wxxx = 2 2h3x − 3hhxhxx + h2hxxx h3 , wxxxx = 2 −6h4x + 12hh2xhxx − 3h2h2xx − 4hxh 2hxxx h4 , wyy = 2 hhyy − h2y h2 , which converts Eq.(27) into a bilinear equation in h as hhxt − hxht + 3h2xx − 4hxhxxx + hhxxxx − hhyy + h2y = 0. (38) Since Hirota bilinear operator is defined as Dp aD q bD r c(h, k) = lim a′→a, b′→b, c′→c ( ∂ ∂a − ∂ ∂a′ )p( ∂ ∂b − ∂ ∂b′ )q ( ∂ ∂c − ∂ ∂c′ )r h(a, b, c)k(a′, b′, c′). let p = 1,q = 0 and r then= 1, DxDt(h.k) = hxtk − htkx − hxkt + hkxt, DxDt(h.h) = 2 (hhxt − hxht) let p = 4, q = 0 and r = 0, then D4 x(h.k) = hxxxxk − 4hxxxkx + 6hxxkxx − 4hxkxxx + hkxxxx D4 x(h.h) = 2 ( 3h2xx − 4hxhxxx + hhxxxx ) let p = 0, q = 2 and r = 0, then D2 y(h.k) = hyyk + hkyy + 2hyky D2 y(h.h) = 2 ( hhyy − h2y ) Therefore, using bilinear differentials D, the bilinear Eq. (38) can be expressed in Hirota’s bilinear form as[ DxDt +D4 x −D2 y ] h.h (39)= 0 This equation (39) is called the Hirota bilinear form for KP equation (27). 3.4 KP equation with variable coefficients in (2+1)-dimensions The nonlinear KP equation with a variable coefficient represents a generalization of the standard KP equation [25]. Its mathematical form is (ut + uux + uxxx)x + g(t)uxy + 3uyy = 0, (40) where u denotes the wave and x, y, t are spatial and temporal coordinates, and g(t) is a time-dependent coefficient that modulates the coupling between the x and y directions. To determine the dispersion relation, we introduce the phase variable ξi = µix+ viy − di(t), (41) https://internationalpubls.com 1221 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) where µi, vi (i = 1, 2, . . .) are constants and di(t) is the time-dependent dispersion function. Substituting u = eξi into the linear terms of Eq. (40) yields (ut + uxxx)x + g(t)uxy + 3uyy = 0. The relevant partial derivatives are ut = −d′i(t)eξi , ux = µieξi , uxxx = µ3i e ξi , uyy = v2i e ξi , uxy = µivieξi . Substituting these expressions into the linearized equation gives ( −d′i(t)eξi + µ3i e ξi ) x + 3v2i e ξi + g(t)µivie ξi = 0, which simplifies to (−d′i(t) + µ3i )µi + 3v2i + g(t)µivi = 0. Solving for d′i(t) leads to d′i(t) = µ3i + 3v2i µi + g(t)vi, and integrating over time gives the dispersion relation di(t) = ∫ ( µ3i + g(t)vi + 3v2i µi ) dt. (42) To construct solutions, we apply the transformation u(x, y, t) = K(lnh)xx, (43) which can equivalently be expressed as u = wxx, where w = K(lnh). (44) Choosing the function h as h = 1 + eξ1 , with ξ1 = µ1x+ v1y − ∫ ( µ31 + g(t)v1 + 3v21 µ1 ) dt, and substituting into Eq. (40), we find that K = 12. Therefore, the transformation (43) reduces to u = 12(lnh)xx. (45) Now, from (44), we have ut = wxxt, ux = wxxx, uxxx = wxxxxx, uyy = wxxyy and uxy = wxxxy putting the above expressions into Eq.(40), we get (wxxt + wxxwxxx + wxxxxx)x + 3wxxyy + g(t)wxxxy (46)= 0 https://internationalpubls.com 1222 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) on integrating w.r.t. x wxxt + wxxwxxx + wxxxxx + 3wxyy + g(t)wxxy = 0. (47) on again integrating w.r.t. x wxt + ∫ wxxwxxx ∂x+ wxxxx + 3wyy + g(t)wxy = 0. (48) We compute integral term in equation (48) with constant of integration as zero I = ∫ wxxwxxx ∂x = 1 2 ∫ 2wxxwxxx ∂x = 1 2 w2 xx, (49) substituting the value of I in equation (48), we get wxt + w2 xx 2 + wxxxx + 3wyy + g(t)wxy = 0. (50) As we have w = 12(log h), we can get the followings: wx = 12 hx h , wxt = 12 hhxt − hxht h2 , wxx = 12 hhxx − h2x h2 , wxxx = 12 2h3x − 3hhxhxx + h2hxxx h3 , wxxxx = 12 −6h4x + 12hh2xhxx − 3h2h2xx − 4hxh 2hxxx h4 , wyy = 12 hhyy − h2y h2 , wxy = 12 hhxy − hxhy h2 , which converts Eq.(40) into a bilinear equation in h as hhxt − hxht + 3h2xx − 4hxhxxx + hhxxxx + 3hhyy − 3h2y − g(t)hyhx + g(t)hhxy = 0. (51) Since Hirota bilinear operator is defined as Dp aD q bD r c(h, k) = lim a′→a, b′→b, c′→c ( ∂ ∂a − ∂ ∂a′ )p( ∂ ∂b − ∂ ∂b′ )q ( ∂ ∂c − ∂ ∂c′ )r h(a, b, c)k(a′, b′, c′). let p = 1,q = 0 and r then= 1, DxDt(h.k) = hxtk − htkx − hxkt + hkxt, DxDt(h.h) = 2 (hhxt − hxht) let p = 4, q = 0 and r = 0, then Dx 4(h.k) = hxxxxk − 4hxxxkx + 6hxxkxx − 4hxkxxx + hkxxxx https://internationalpubls.com 1223 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) D4 x(h.h) = 2 ( 3h2xx − 4hxhxxx + hhxxxx ) let p = 0, q = 2 and r then= 0, D2 y(h.k) = hyyk + hkyy + 2hyky D2 y(h.h) = 2 ( hhyy − h2y ) let p = 1,q = 1 and r = 0, then DxDy(h.k) = hxyk − hykx − hxky + hkxy, DxDy(h.h) = 2 (hhxy − hxhy) Therefore, using bilinear differentials D, the Eq. (51) gives the Hirota’s bilinear form as[ DxDt +D4 x + 3D2 y + g(t)DxDy ] h.h (52)= 0 This equation (52) is called the Hirota bilinear form for KP equation with variable coefficient (40). 3.5 Graphene-Sheets Equation in (2+1)-dimensions Let us consider a variable-coefficient equation that models the thermophoretic waves in graphene sheets [26]. The equation is expressed as uxt + ( uux + uxxx + (α(t) + β)ux ) x + γ(t)uyy = 0, (53) where u = u(x, y, t) represents the thermophoretic displacement, x and y are the longitudinal and lateral coordinates, t is time, α(t) and β denote coefficient of thermal conductivity and γ(t) as lateral dispersion coefficient. We take the phase variable as ξi = µix+ viy − di(t), (54) where µi, vi (i = 1, 2, . . .) are constants and di(t) is a time-dependent dispersion function. Substituting u = eξi into the linear part of Eq. (53) gives uxt + uxxxx + (α(t) + β)uxx + γ(t)uyy = 0. The relevant derivatives are ut = −di′(t)u, uxxxx = µi 4u, uxx = µi 2u, uyy = vi 2u. Substituting these into the linearized equation leads to −di′(t)u+ µi 4u+ (α(t) + β)µi 2u+ γ(t)vi 2u = 0, which simplifies to −di′(t) + µi 4 + µi 2α(t) + µi 2β + γ(t)vi 2 = 0. Solving for di ′(t) gives di ′(t) = µi 4 + µi 2α(t) + µi 2β + γ(t)vi 2, https://internationalpubls.com 1224 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) and integrating over time yields the dispersion relation di(t) = ∫ µ4i + µ2iα(t) + γ(t)v2i + µ2iβ µi dt. (55) To construct solutions, we employ thetransformation u(x, y, t) = K(lnh)xx, (56) which can equivalently be written as u = wxx, where w = K(lnh). (57) Choosing the auxiliary function h(x, y, t) = 1 + eξ1 , with ξ1 = µ1x+ v1y − ∫ µ41 + µ21α(t) + γ(t)v21 + µ21β µ1 dt, and substituting into Eq. (53), we find K = 12. Therefore, the transformation (56) reduces to u = 12(lnh)xx. (58) Now, from (57), we have uxt = wxxxt, ux = wxxx, uxxx = wxxxxx and uyy = wxxyy putting the above expressions into Eq.(53), we get wxxxt + (wxxwxxx + wxxxxx + (α(t) + β)wxxx)x + γ(t)wxxyy (59)= 0 on integrating w.r.t. x wxxt + wxxwxxx + wxxxxx + (α(t) + β)wxxx + γ(t)wxyy (60)= 0 on again integrating w.r.t. x wxt + ∫ wxxwxxx dx+ wxxxx + (α(t) + β)wxx + γ(t)wyy = 0. (61) We compute integral term in equation (61) with constant of integration as zero I = ∫ wxxwxxx ∂x = 1 2 ∫ 2wxxwxxx ∂x = 1 2 w2 xx, (62) substituting the value of I in equation (61), we get wxt + w2 xx 2 + wxxxx + (α(t) + β)wxx + γ(t)wyy = 0. (63) As we have w = 12(log h), we can get the followings: wx = 12 hx h , wxt = 12 hhxt − hxht h2 , https://internationalpubls.com 1225 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) wxx = 12 hhxx − h2x h2 , wxxx = 12 2h3x − 3hhxhxx + h2hxxx h3 , wxxxx = 12 −6h4x + 12hh2xhxx − 3h2h2xx − 4hxh 2hxxx h4 , wyy = 12 hhyy − h2y h2 , which converts Eq.(53) into a bilinear equation in h as (hhxt − hxht) + (β + α(t))(hhxx + h2x) + γ(t)(hhyy − h2y) + (hhxxxx − 4hxxxhx + 3h2xx) = 0. (64) Since Hirota bilinear operator is defined as Dp aD q bD r c(h, k) = lim a′→a, b′→b, c′→c ( ∂ ∂a − ∂ ∂a′ )p( ∂ ∂b − ∂ ∂b′ )q ( ∂ ∂c − ∂ ∂c′ )r h(a, b, c)k(a′, b′, c′). let p = r = 1 and q = 0 then DxDt(h.k) = hxtk − htkx − hxkt + hkxt, DxDt(h.h) = 2 (hhxt − hxht) let p = 2 and q = r = 0then D2 x(h.k) = hxxk + hkxx + 2hxkx D2 x(h.h) = 2 ( hhxx − h2x ) let p = r = 0 and q = 2 then D2 y(h.k) = hyyk + hkyy + 2hyky D2 y(h.h) = 2 ( hhyy − h2y ) let p = 4 and q = r = 0 then D4 x(h.k) = hxxxxk − 4.hxxxkx + 6.hxxkxx − 4.hxkxxx + hkxxxx D4 x(h.h) = 2 ( hhxxxx − 4.hxhxxx + 3.h2xx ) Therefore, using bilinear differentials D, the bilinear Eq. (64) can be expressed in Hirota’s bilinear form as[ DxDt + (β + α(t))D2 x + γ(t)D2 y +D4 x ] h.h (65)= 0 This equation (65) is called the Hirota bilinear form for Graphene sheets equation (53). https://internationalpubls.com 1226 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 3.6 Bogoyavlenskii–Kadomtsev–Petviashvili (BKP) Equation in (3+1)-dimensions The integrable BKP equation is given by [27] uyt + 3uxz − 3uxuxy − 3uxxuy − uxxxy = 0, (66) where x, y, z, t are spatial and temporal coordinates, and u represents the wave amplitude. To analyze the dispersion properties, we define the phase variable ξi = µix+ viy + wiz − dit, (67) where µi, vi, wi (i = 1, 2, . . .) are constants, and di is the dispersion coefficient. Substituting u = eξi into the linear terms of Eq. (66), namely uyt + 3uxz − uxxxy = 0, yields the dispersion relation di = −3µiwi + µ3i vi vi . (68) Next, we employ the transformation u(x, y, z, t) = K(lnh)x, (69) and select the h function as h = 1 + eξ1 , with ξ1 = µ1x+ v1y + w1z − ( −3µ1w1 + µ31v1 v1 ) t. Substituting into Eq. (66) and simplifying gives K = 2. Therefore, the transformation (69) reduces to u = 2(lnh)x. (70) Now, from (70), we have ux = 2 ( hxxh− h2x h2 ) , uy = 2 ( hxyh− hxhy h2 ) , uz = 2 ( hxzh− hxhz h2 ) , ut = 2 ( hxth− hxht h2 ) . The higher-order derivatives become uyt = 2 (hxyth+ hxyht − hxthy − hxhyt)h 2 − 2hht(hxyh− hxhy) h4 , uxz = 2 (hxxzh+ hxxhz − 2hxhxz)h 2 − 2hhz(hxxh− h2x) h4 , uxx = 2 hxxxh 3 − 3hxhxxh 2 + 2h3xh h4 , uxy = 2 (hxxyh+ hxxhy − 2hxhxy)h 2 − 2hhy(hxxh− h2x) h4 , uxxxy = 2 hxxxy h − 2 hxxxxhy h2 − 8 hxxxyhx h2 − 8 hxxxhxy h2 + 16 hxxxhxhy h3 − 12 hxxhxxy h2 + 12 h2xxhy h3 + 36 hxhxyhxx h3 + 18 h2xhxxy h3 − 54 h2xhxxhy h4 − 16 h3xhxy h4 + 16 hx 4hy h5 . https://internationalpubls.com 1227 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) Substituting these expressions into Eq. (66), we get hhyt − hyht + 3 (hhxz − hxhz) − hxxxyh+ 3hxxyhx − 3hxyhxx + hyhxxx = 0. (71) Since the Hirota bilinear operator is defined as Dp aD q bD r cD s d(h, k) = lim a′→a, b′→b, c′→c, d′→d ( ∂ ∂a − ∂ ∂a′ )p( ∂ ∂b − ∂ ∂b′ )q × ( ∂ ∂c − ∂ ∂c′ )r ( ∂ ∂d − ∂ ∂d′ )s h(a, b, c, d)k(a′, b′, c′, d′), , we evaluate the following cases: DyDt(h.h) = 2 (hhyt − hyht) D, xDz(h.h) = 2 (hhxz − hxhz) , (D3 xDy)(h.h) = 2 (hxxxyh− 3hxxyhx + 3hxyhxx − hyhxxx) . Thus, Eq. (71) gives Hirota bilinear form as( DyDt + 3DxDz −D3 xDy ) h · h = 0. (72) This equation (72) is called the Hirota bilinear form for the BKP equation (66). 3.7 Boiti–Leon–Manna–Pempinelli equation in (3+1)-dimensions The integrable BLMP equation is expressed as [28] uyt + uzt + uxxxy + uxxxz − 3uxuxy − 3uxuxz − 3uxxuy − 3uxxuz = 0, (73) where x, y, z, t are spatial and temporal coordinates, and u represents the wave amplitude. To determine the dispersion relation, we introduce the phase variable ξi = µix+ viy + wiz − dit, (74) where µi, vi, wi (i = 1, 2, . . .) are constants and di is the dispersion coefficient. Substituting u = eξi into the linear part of Eq. (73), namely uyt + uzt + uxxxy + uxxxz = 0, the relevant derivatives are uyt = vi(−di)eξi , uzt = wi(−di)eξi , uxxxy = µ3i vie ξi , uxxxz = µ3iwie ξi . Substituting these expressions into the linearized equation yields −(vi + wi)di + µ3i (vi + wi) = 0, which leads to the dispersion relation di = µ3i . (75) https://internationalpubls.com 1228 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) Next, we employ the transformation u(x, y, z, t) = K(lnh)x, (76) and take h function as h = 1 + eξ1 , with ξ1 = µ1x+ v1y + w1z − µ31t. Substituting into Eq. (73) and simplifying, we find K = −2. Therefore, the transformation (76) reduces to u = −2(lnh)x. (77) Now, from (77), we have ux = 2 ∂ ∂x ( hx h ) = 2 ( hxxh− h2x h2 ) uy = 2 ∂ ∂y ( hy h ) = 2 ( hxyh− hxhy h2 ) uz = 2 ∂ ∂z ( hz h ) = 2 ( hxzh− hxhz h2 ) ut = 2 ∂ ∂t ( ht h ) = 2 ( hxth− hxht h2 ) uyt = 2 (hxyth+ hxyht − hxthy − hxhyt)h 2 − 2hht(hxyh− hxhy) h4 uxz = 2 (hxxzh+ hxxhz − 2hxhxz)h 2 − 2hhz(hxxh− h2x) h4 uxx = 2 hxxxh 3 − 3hxhxxh 2 + 2h3xh h4 uxy = 2 (hxxyh+ hxxhy − 2hxhxy)h 2 − 2hhy(hxxh− h2x) h4 uxxxy = −2 hxxxxy h + 2 hxxxxhy h2 + 8 hxxxyhx h2 + 8 hxxxhxy h2 − 16 hxxxhxhy h3 +12 hxxhxxy h2 − 12 h2xxhy h3 − 36 hxhxyhxx h3 − 18 h2xhxxy h3 + 54 h2xhxxhy h4 +16 h3xhxy h4 − 16 h4xhy h5 uxxxz = −2 hxxxxz h + 2 hxxxxhz h2 + 8 hxxxzhx h2 + 8 hxxxhxz h2 − 16 hxxxhxhz h3 +12 hxxhxxz h2 − 12 h2xxhz h3 − 36 hxhxzhxx h3 − 18 h2xhxxz h3 + 54 h2xhxxhz h4 +16 h3xhxz h4 − 16 h4xhz h5 putting the above expressions into Eq.(73), we get hhyt−hyht+hhzt−hzht+hhxxxy−3hxhxxy+3hxyhxx−hyhxxx+hhxxxz−3hxhxxz+3hxzhxx−hzhxxx = 0 (78) https://internationalpubls.com 1229 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) Since Hirota bilinear operator is defined as Dp aD q bD r cD s d(h, k) = lim a′→a, b′→b, c′→c, d′→d ( ∂ ∂a − ∂ ∂a′ )p( ∂ ∂b − ∂ ∂b′ )q ( ∂ ∂c − ∂ ∂c′ )r ( ∂ ∂d − ∂ ∂d′ )s h(a, b, c, d)k(a′, b′, c′, d′). Let p = 0,q = 1,r = 0 and s then= 1, DyDt(h.k) = hytk − htky − hykt + hkyt, DyDt(h.h) = 2 (hhyt − hyht) let p = 0,q = 0,r = 1 and s then= 1, DzDt(h.k) = hztk − htkz − hzkt + hkyz, DzDt(h.h) = 2 (hhzt − hzht) let p = 3,q = 1,r = 0 and s = 0, then (D3 xDy)(hk) = hxxxyk − 3.hxxykx + 3.hxykxx − hykxxx − hxxxky + 3.hxxkxy − 3.hxkxxy + hkxxxy (D3 xDy)(hh) = 2 (hxxxyh− 3.hxxyhx + 3.hxyhxx − hyhxxx) let p = 3,q = 0,r = 1 and s then= 0, (D3 xDz)(hk) = hxxxzk − 3.hxxzkx + 3.hxzkxx − hzkxxx − hxxxkz + 3.hxxkxz − 3.hxkxxz + hkxxxz (D3 xDz)(hh) = 2 (hxxxzh− 3.hxxzhx + 3.hxzhxx − hzhxxx) Equation (78) makes D-operator form as (DyDt +DzDt +D3 xDy +D3 xDz)h.h (79)= 0 This equation (79) is called the Hirota bilinear form for BLMP equation (73). 4 Results and Analysis In above examples, each nonlinear equation was analyzed for Hirota’s bilinear equation and its bilinear form. First, the original equation was transformed using a Cole-Hopf transformation that converts the equation to a simplified bilinear equation. Thereafter, Hirota’s D-operators were applied to express the equation in bilinear operator notation. This approach works for equations with multiple dimensions, variable coefficients, or higher-order derivatives. The bilinear forms provide a convenient framework for obtaining exact solutions with the symbolic computation. 4.1 KdV Equation The equation in u(x, t) as ut + 6uux + uxxx = 0. Using the Cole-Hopf transformation, the solution is expressed as u = 2(lnh)xx. that gives a bilinear equation as hhxt − hxht + 3h2xx − 4hxhxxx + hhxxxx = 0, which gives Hirota’s bilinear operator form as (DxDt +Dx 4)h · h = 0. https://internationalpubls.com 1230 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 4.2 Boussinesq Equation The Boussinesq equation is utt − uxx − 3(u2)xx − uxxxx = 0. Using the Cole-Hopf transformation, the solution is expressed as u = 2(ln h)xx, that gives a bilinear equation as hhtt − h2t − hhxx + h2x − hhxxxx + 4hxhxxx − 3h2xx = 0, which gives Hirota’s bilinear operator form as (D2 t −D2 x −D4 x)h · h = 0. 4.3 KP Equation The equation is (ut + 6uux + uxxx)x − uyy = 0. Using the Cole-Hopf transformation, the solution is expressed as u = 2(lnh)xx, that gives a bilinear equation as hhxt − hxht + 3h2xx − 4hxhxxx + hhxxxx − hhyy + h2y = 0. which gives Hirota’s bilinear operator form as (DxDt +D4 x −D2 y)h · h = 0. 4.4 KP Equation with Variable Coefficient The vc-KP equation is (ut + uux + uxxx)x + g(t)uxy + 3uyy = 0. Using the transformation, the solution is expressed as u = 12(lnh)xx, that gives a bilinear equation as hhxt − hxht + 3h2xx − 4hxhxxx + hhxxxx + 3hhyy − 3h2y − g(t)hyhx + g(t)hhxy = 0, which gives Hirota’s bilinear operator form as (DxDt +D4 x + 3D2 y + g(t)DxDy)h · h = 0. https://internationalpubls.com 1231 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 4.5 Graphene-Sheets Equation The Graphene-Sheets equation is uxt + (uux + uxxx + (α(t) + β)ux)x + γ(t)uyy = 0. Using the transformation, the solution is expressed as u = 12(lnh)xx, that gives a bilinear equation as (hhxt − hxht) + (β + α(t))(hhxx + h2x) + γ(t)(hhyy − h2y) + (hhxxxx − 4hxhxxx + 3h2xx) = 0, which gives Hirota’s bilinear operator form as (DxDt + (β + α(t))D2 x + γ(t)D2 y +D4 x)h · h = 0. 4.6 BKP Equation The BKP equation is uyt + 3uxz − 3uxuxy − 3uxxuy − uxxxy = 0. Using the transformation, the solution is expressed as u = 2(lnh)x, that gives a bilinear equation as hhyt − hyht + 3(hhxz − hxhz) − hhxxxy + 3hxxyhx − 3hxyhxx + hyhxxx = 0, which gives Hirota’s bilinear operator form as (DyDt + 3DxDz −D3 xDy)h · h = 0. 4.7 BLMP Equation The equation is uyt + uzt + uxxxy + uxxxz − 3uxuxy − 3uxuxz − 3uxxuy − 3uxxuz = 0. Using the transformation, the solution is expressed as u = −2(lnh)x, that gives a bilinear equation as hhyt−hyht+hhzt−hzht+hhxxxy−3hxhxxy +3hxyhxx−hyhxxx+hhxxxz−3hxhxxz +3hxzhxx−hzhxxx = 0, which gives Hirota’s bilinear operator form as (DyDt +DzDt +D3 xDy +D3 xDz)h · h = 0. https://internationalpubls.com 1232 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 5 Conclusions In this manuscript, we explored the symbolic algorithm for creating the Hirota’a bilinear form for nonlinear PDEs in higher-dimensions. Algorithm found the Cole-Hopf transformation by balancing nonlinear and higher order terms. It constructed the bilinear equations using the obtained transformations. Next, it con- structed the Hirota’s D-operators so that it could convert the bilinear equation to bilinear D-operator form. After converting nonlinear PDEs into bilinear forms,we can obtain exact solutions such as lumps, breathers, and solitons. We successfully used the algorithm to obtain the Hirota bilinear form for several well-known equations such as the KdV, Boussinesq, KP, and BLMP equations and other nonlinear equations. The resulting bilinear forms serve as a basis for additional research, such as the development of multi-soliton solutions and the investigation of interaction phenomena. The algorithm can be further streamlined and made accessible for more complex systems by utilizing symbolic computation tools such as Maple or Matlab. Declarations Competing interests There is no conflict of interest, according to the authors. Authors’ contributions Each author made an equal contribution to the final draft of the work. The authors would have consented and approved the final work. Data availability statement Not applicable to this research as no data were analyzed and created in this work. References [1] S. Kumar and B. Mohan, “A novel and efficient method for obtaining hirota’s bilinear form for the non- linear evolution equation in (n+ 1) dimensions,” Partial Differential Equations in Applied Mathematics, vol. 5, p. 100274, 2022. [2] A. Parihar and D. Malik, “Application of hirota’s direct method to nonlinear partial differential equa- tions: Bilinear form and soliton solutions,” 2022. [3] R. Hirota, “Exact solution of the korteweg—de vries equation for multiple collisions of solitons,” Physical Review Letters, vol. 27, no. 18, p. 1192, 1971. [4] R. Hirota, The direct method in soliton theory. No. 155, Cambridge university press, 2004. [5] W.-X. Ma, “N-soliton solutions and the hirota conditions in (2+ 1)-dimensions,” Optical and Quantum Electronics, vol. 52, no. 12, p. 511, 2020. [6] A. C. Newell, Solitons in mathematics and physics. SIAM, 1985. [7] Z. Du, B. Tian, X.-Y. Xie, J. Chai, and X.-Y. Wu, “Bäcklund transformation and soliton solutions in terms of the wronskian for the kadomtsev–petviashvili-based system in fluid dynamics,” Pramana, vol. 90, pp. 1–6, 2018. https://internationalpubls.com 1233 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) [8] X.-W. Yan, S.-F. Tian, M.-J. Dong, and L. Zou, “Bäcklund transformation, rogue wave solutions and interaction phenomena for a (3+ 1)(3+ 1)-dimensional b-type kadomtsev–petviashvili–boussinesq equation,” Nonlinear Dynamics, vol. 92, pp. 709–720, 2018. [9] X. Guan, W. Liu, Q. Zhou, and A. Biswas, “Darboux transformation and analytic solutions for a generalized super-nls-mkdv equation,” Nonlinear Dynamics, vol. 98, no. 2, pp. 1491–1500, 2019. [10] X. Wang and J. Wei, “Three types of darboux transformation and general soliton solutions for the space- shifted nonlocal pt symmetric nonlinear schrödinger equation,” Applied Mathematics Letters, vol. 130, p. 107998, 2022. [11] Y. Shen, B. Tian, T.-Y. Zhou, and C.-D. Cheng, “Localized waves of the higher-order nonlinear schrödinger-maxwell-bloch system with the sextic terms in an erbium-doped fiber,” Nonlinear Dy- namics, vol. 112, no. 2, pp. 1275–1290, 2024. [12] Y. Shen, B. Tian, T.-Y. Zhou, and C.-D. Cheng, “Complex kraenkel-manna-merle system in a ferrite: N- fold darboux transformation, generalized darboux transformation and solitons,” Mathematical Modelling of Natural Phenomena, vol. 18, p. 30, 2023. [13] I. Hamid and S. Kumar, “Symbolic computation and novel solitons, traveling waves and soliton-like solutions for the highly nonlinear (2+ 1)-dimensional schrödinger equation in the anomalous dispersion regime via newly proposed modified approach,” Optical and Quantum Electronics, vol. 55, no. 9, p. 755, 2023. [14] S. Kumar, W.-X. Ma, S. K. Dhiman, and A. Chauhan, “Lie group analysis with the optimal system, gen- eralized invariant solutions, and an enormous variety of different wave profiles for the higher-dimensional modified dispersive water wave system of equations,” The European Physical Journal Plus, vol. 138, no. 5, p. 434, 2023. [15] A.-M. Wazwaz, “The simplified hirota’s method for studying three extended higher-order kdv-type equations,” Journal of Ocean Engineering and Science, vol. 1, no. 3, pp. 181–185, 2016. [16] W. Hereman and A. Nuseir, “Symbolic methods to construct exact solutions of nonlinear partial differ- ential equations,” Mathematics and Computers in Simulation, vol. 43, no. 1, pp. 13–27, 1997. [17] Y. Shen, B. Tian, C.-D. Cheng, and T.-Y. Zhou, “Pfaffian solutions and nonlinear waves of a (3+ 1)-dimensional generalized konopelchenko–dubrovsky–kaup–kupershmidt system in fluid mechanics,” Physics of Fluids, vol. 35, no. 2, 2023. [18] C.-D. Cheng, B. Tian, Y. Shen, and T.-Y. Zhou, “Bilinear form, auto-bäcklund transformations, pfaf- fian, soliton, and breather solutions for a (3+ 1)-dimensional extended shallow water wave equation,” Physics of Fluids, vol. 35, no. 8, 2023. [19] P. Capetillo and J. Hornewall, “Introduction to the hirota direct method,” 2021. [20] E. Hopf, “The partial differential equation,” 1950. [21] J. D. Cole, “On a quasi-linear parabolic equation occurring in aerodynamics,” Quarterly of applied mathematics, vol. 9, no. 3, pp. 225–236, 1951. [22] R. Hirota, The direct method in soliton theory. No. 155, Cambridge university press, 2004. [23] A.-M. Wazwaz, “Multiple-soliton solutions for the boussinesq equation,” Applied Mathematics and Computation, vol. 192, no. 2, pp. 479–486, 2007. https://internationalpubls.com 1234 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) [24] W.-X. Ma, “N-soliton solutions and the hirota conditions in (2+ 1)-dimensions,” Optical and Quantum Electronics, vol. 52, no. 12, p. 511, 2020. [25] S. Kumar and B. Mohan, “A study of multi-soliton solutions, breather, lumps, and their interactions for kadomtsev-petviashvili equation with variable time coeffcient using hirota method,” Physica Scripta, vol. 96, no. 12, p. 125255, 2021. [26] R. M. El-Shiekh and M. Gaballah, “Bilinear form and n-soliton thermophoric waves for the variable coefficients (2+ 1)-dimensional graphene sheets equation,” Optical and Quantum Electronics, vol. 56, no. 5, p. 872, 2024. [27] W.-X. Ma and Z. Zhu, “Solving the (3+ 1)-dimensional generalized kp and bkp equations by the multiple exp-function algorithm,” Applied Mathematics and Computation, vol. 218, no. 24, pp. 11871– 11879, 2012. [28] M. R. Ali and W.-X. Ma, “New exact solutions of nonlinear (3+ 1)-dimensional boiti-leon-manna- pempinelli equation,” Advances in Mathematical Physics, vol. 2019, no. 1, p. 9801638, 2019. https://internationalpubls.com 1235 Introduction Hirota Bilinear Form of a nonlinear PDE Cole–Hopf transformations The Hirota Bilinear Operator Algorithm for bilinear Form Applications of Hirota Bilinear form KdV equation in (1+1)-dimensions Boussinesq Equation in (1+1)-dimensions KP equation in (2+1)-dimensions KP equation with variable coefficients in (2+1)-dimensions Graphene-Sheets Equation in (2+1)-dimensions Bogoyavlenskii–Kadomtsev–Petviashvili (BKP) Equation in (3+1)-dimensions Boiti–Leon–Manna–Pempinelli equation in (3+1)-dimensions Results and Analysis KdV Equation Boussinesq Equation KP Equation KP Equation with Variable Coefficient Graphene-Sheets Equation BKP Equation BLMP Equation Conclusions