Next Article in Journal / Special Issue
Lie Symmetry Analysis, Closed-Form Solutions, and Conservation Laws for the Camassa–Holm Type Equation
Previous Article in Journal
Modeling of Sedimentation of Particles near Corrugated Surfaces by the Meshless Method of Fundamental Solutions
Previous Article in Special Issue
Periodic and Axial Perturbations of Chaotic Solitons in the Realm of Complex Structured Quintic Swift-Hohenberg Equation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Using Symmetries to Investigate the Complete Integrability, Solitary Wave Solutions and Solitons of the Gardner Equation

by
Willy Hereman
1,*,† and
Ünal Göktaş
2,†
1
Department of Applied Mathematics and Statistics, Colorado School of Mines, Golden, CO 80401-1887, USA
2
Department of Computer Science and Engineering, Texas A&M University, College Station, TX 77843-3112, USA
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Math. Comput. Appl. 2024, 29(5), 91; https://doi.org/10.3390/mca29050091
Submission received: 26 July 2024 / Revised: 26 September 2024 / Accepted: 27 September 2024 / Published: 3 October 2024
(This article belongs to the Special Issue Symmetry Methods for Solving Differential Equations)

Abstract

:
In this paper, using a scaling symmetry, it is shown how to compute polynomial conservation laws, generalized symmetries, recursion operators, Lax pairs, and bilinear forms of polynomial nonlinear partial differential equations, thereby establishing their complete integrability. The Gardner equation is chosen as the key example, as it comprises both the Korteweg–de Vries and modified Korteweg–de Vries equations. The Gardner and Miura transformations, which connect these equations, are also computed using the concept of scaling homogeneity. Exact solitary wave solutions and solitons of the Gardner equation are derived using Hirota’s method and other direct methods. The nature of these solutions depends on the sign of the cubic term in the Gardner equation and the underlying mKdV equation. It is shown that flat (table-top) waves of large amplitude only occur when the sign of the cubic nonlinearity is negative (defocusing case), whereas the focusing Gardner equation has standard elastically colliding solitons. This paper’s aim is to provide a review of the integrability properties and solutions of the Gardner equation and to illustrate the applicability of the scaling symmetry approach. The methods and algorithms used in this paper have been implemented in Mathematica, but can be adapted for major computer algebra systems.
PACS:
Primary 35Q51, 35Q53; Secondary 37K10, 37K40

1. Introduction

Several physically important nonlinear partial differential equations (PDEs) are completely integrable by the inverse scattering transform (IST) method, which can be viewed as the nonlinear analog of the Fourier transform method (see [1,2,3] and references therein). Arguably, the most well-known completely integrable PDEs are the Korteweg–de Vries (KdV), modified KdV (mKdV), nonlinear Schrödinger, and sine-Gordon equations. They model wave phenomena in a wide range of applications in modern theoretical and mathematical physics, including fluid dynamics, nonlinear optics, and plasma physics, to name a few.
“Complete integrability” is an elusive term [4], but completely integrable equations have remarkable properties and a rich mathematical structure. For instance, they possess infinitely many conservation laws and high-order symmetries. They admit a Lax pair where the nonlinear PDE is replaced by a pair of linear equations whose compatibility only holds on solutions of the nonlinear PDE. They can be viewed as infinite dimensional Hamiltonian systems with two or three different Hamiltonians. By applying a change of dependent variable, completely integrable PDEs can be transformed into equations that are homogeneous of degree (second or higher) in that new variable and be recast in bilinear form in terms of Hirota operators [5,6]. Most importantly, completely integrable PDEs have solitary waves solutions that maintain their shape and speed while propagating at a constant velocity, and soliton solutions made up of such solitary waves that collide elastically. Nowadays, the terms solitary wave and soliton are used interchangeably, without reference to the elastic scattering property.
To obtain an initial idea about the possible complete integrability of a PDE, one should check if the equation has the Painlevé property [1,2], meaning that its solutions have no worse singularities than movable poles. To do so, one runs the so-called Painlevé test [7], which is algorithmic and can be performed with our Mathematica code PainleveTest.m [8]. If a PDE passes the Painlevé test, there is no guarantee that it is completely integrable, but it is more likely than not to have special properties, viz., conservation laws, generalized symmetries, and so on.
In this paper we use the concept of scaling homogeneity to symbolically compute conserved densities, fluxes, higher-order symmetries, Lax pairs, bilinear forms, and Miura-type transformations for polynomial systems of PDEs, thereby testing their complete integrability in a variety of ways. To keep the paper accessible to non-experts, we only cover PDEs in ( 1 + 1 ) dimensions, i.e., one space variable ( x ) in addition to time ( t ) . To illustrate the steps of the computations we focus on the Gardner equation, which comprises the KdV and mKdV equations as special cases. Therefore, the results for these two important soliton equations are obtained without extra effort. This paper’s purpose is to provide a review of the integrability properties and solutions of the Gardner equation and illustrate the applicability of the scaling symmetry approach for investigating the complete integrability of polynomial nonlinear PDEs.
Scaling homogeneity is a feature common to many nonlinear PDEs, and it provides an elegant way to find, e.g., densities and higher-order symmetries. Indeed, these scalar quantities can be derived as linear combinations with undetermined (constant) coefficients of scaling homogeneous polynomial terms in the dependent variables and their derivatives. Since the defining equations for these quantities are linear, the method comes down to solving linear systems for the undetermined coefficients.
For conservation laws, the time derivative of a candidate density must be the total derivative of an unknown flux. To test this type of “exactness” we use the Euler operator from the calculus of variations. An automated computation of the corresponding fluxes requires integration by parts on the so-called jet space. In essence, one needs to integrate expressions involving arbitrary functions. Unfortunately, computer algebra systems (CAS) often fail at that task. To circumvent this shortcoming, we use the homotopy operator from differential geometry, which reduces the problem to a one-dimensional integral with respect to an auxiliary parameter.
The existence of a sufficiently large number of non-trivial densities suffices to establish complete integrability of a given PDE. In most cases, the corresponding fluxes are not needed. However, for some equations, e.g., the Kadomtsev–Petviashvili equation [1,2], it is necessary to swap the independent variables to be able to rewrite the PDE as a system of evolution equations. In doing so, the roles of density and flux become interchanged, and that is why we also show the computation of the flux.
The computation of higher-order symmetries is performed with the same methodology without having to rely on the relationship between symmetries and conservation laws, as expressed through Noether’s theorem. The scaling symmetry argument and the method of undetermined coefficients also apply to recursion operators, which connect the generalized symmetries. Of course, the recursion operator is hard to compute because it involves total derivatives and anti-derivatives. Fortunately, the defining equation for the recursion operator is linear, which again reduces the problem to solving a linear system. By applying the recursion operator to a low-order generalized symmetry, one can generate all higher-order symmetries one after another. Thus, the existence of a recursion operator provides proof that a nonlinear PDE is completely integrable. With the recursion operator, one can build a hierarchy of completely integrable equations which have the same properties as the original nonlinear PDE. For example, the well-known Lax equation [5,9], which is of fifth order in x, is the first higher-order generalized symmetry of the KdV equation [10]. When encountering a allegedly “new” nonlinear PDE of high-order in the literature, one should verify that it does not belong to the hierarchy of symmetries of a well-studied PDE of lower order.
The computation of Lax pairs based on scaling homogeneity is more challenging because the defining equation is no longer linear, and thus the method of undetermined coefficients leads to an algebraic system with quadratic terms. It is still quite straightforward to solve; however, if parameters are present, finding the conditions of the parameters for which solutions exist is substantially harder. The knowledge of a Lax pair is essential to apply the IST, which allows one to solve the initial value problem for a nonlinear PDE and, consequently, find solitary wave solutions and solitons. Obtaining a Lax pair is the first step in the application of the IST and Riemann–Hilbert method to find soliton solutions. Using the Lax pair one can also derive conservation laws and many other properties of the nonlinear PDE.
The same thing happens when applying the scaling homogeneity method to compute the Miura and Gardner transformations. Here, one has to solve quadratic systems for the undetermined coefficients. The Miura transformation, which connects the KdV and mKdV equations, had a profound impact on the discovery of conservation laws of both equations and, in general, the early development of the notion of complete integrability of PDEs and soliton theory.
To apply Hirota’s direct method for finding solitary waves and solitons, the given PDE must be cast in bilinear form. Finding that bilinear representation is a non-trivial task which requires some insight, experience, and often ingenuity as explained in [5], our scaling-homogeneity approach is still helpful, but has its limitations— in particular, if finding a bilinear representation requires splitting an expression into two parts in such a way that each part can be written in bilinear form. This is exactly what must be conducted with regard to the mKdV and Gardner equations and, as far as we know, there is no algorithmic way to achieve this. Once the bilinear representation is found, special solutions, including solitary waves and solitons, can be computed algorithmically.
As we will also show in Section 9, the Gardner equation can be transformed into the mKdV equation via a Galilean transformation. Consequently, any time new solutions of the mKdV equation are discovered, one obtains additional solutions to the Gardner equation, and vice versa [11]. The search for new solutions, in particular for the defocusing mKdV equation, is still an active area of research (see, e.g., refs. [12,13,14,15]). Over the last two decades, our research team has designed and implemented fast algorithms [16,17,18] to test the integrability of nonlinear PDEs based on the concept of scaling symmetry [9,10,19,20,21,22,23]. In addition, we have implemented a simplified version of Hirota’s method [5,24] and other direct methods [25,26] to find exact solutions for nonlinear PDEs. As a matter of fact, the computations in this paper have been largely performed with our software packages, which are written in Mathematica syntax, but could be adapted for other computer algebra systems such as Maple and REDUCE.
The paper is organized as follows. In Section 2 we introduce the Gardner equation and mention some of its applications. As shown in Section 3, the Gardner equation passes the Painlevé test, which indicates that the equation likely has many interesting properties. In Section 4, we discuss the scaling symmetry of the Gardner equation, as well as the KdV and mKdV equations. Section 5 is devoted to the computation of conservation laws. Section 6 and Section 7 cover the computation of generalized symmetries and the recursion operator that connects them. In Section 8, we turn our attention to the computation of a Lax pair for the Gardner equation. The derivation of a bilinear representation is covered in Section 9. The computation of the Gardner transformation (connecting solutions of the Gardner equation and KdV equations) is shown in Section 10.
With the important quantities of the Gardner equation being computed, we show in Section 11 and Section 12 how our results translate to the KdV and mKdV equations. Solitary wave solutions and solitons are covered in Section 13 and Section 14, where we show that the nature of the solutions of the Gardner and mKdV equations depends on the sign of the cubic nonlinearity. A brief discussion of our software packages used in this study is given in Section 15. Finally, in Section 16, we draw some conclusions and mention a few topics for future work.

2. The Gardner Equation

The Gardner equation [27] is the nonlinear PDE
u t + α u u x + β u 2 u x + u 3 x = 0 ,
where u ( x , t ) is the dependent variable (or field) which is a function of space variable x and time t. The subscripts denote partial derivatives, i.e., u t = u t , u x = u x , and u 3 x = 3 u x 3 . The parameters α and β are real numbers for which the values ± 1 and ± 6 are frequently used in the literature.
The Gardner equation has the KdV and mKdV equations as special cases. Therefore, (1) is sometimes called the combined or mixed KdV-mKdV equation. Indeed, for  β = 0 , (1) reduces to the ubiquitous Korteweg–de Vries equation [28],
u t + α u u x + u 3 x = 0 ,
which models, e.g., shallow water waves [29], ion-acoustic waves in plasmas [1,2,30], and many other nonlinear wave phenomena where solitons arise. When α = 0 , (1) becomes the mKdV equation [1,27],
u t + β u 2 u x + u 3 x = 0 ,
which models internal ocean waves, electromagnetic surface waves, waves in plasmas, and more [1].
Equation (1) appears for the first time in [31,32] in the context of the Miura transformation (see Section 10), which connects the mKdV and KdV equations. The Gardner equation has a long history [27,33] and many applications [34] in fluid dynamics: in particular, for modeling long internal water waves [35,36], the dynamics of undular bores [37], and waves in multi-species plasma physics [38,39,40,41]. There is a wealth of information available about the Gardner equation. The equation appears in most books about solitons and integrability (see, e.g., ref. [29], which includes a list of soliton books).
Without loss of generality, we take α 0 in (1), because if α < 0 , replacing u by u would make the coefficient of u u x positive again. However, no discrete symmetries of x , t , or u will flip the sign of the coefficient of u 2 u x , so the cases for positive and negative β will have to be treated separately. Exact solutions of (1) with β > 0 , called the focusing Gardner equation, are quite different from these of the defocusing case, where β < 0 .
We focus on (1) because it is a typical example of a scalar ( 1 + 1 ) -dimensional evolution equation,
u t = F ( x , t , u , u x , u x x , , u n x ) ,
of order n in x and first order in t, with  n = 3 in (1). More importantly, many of the results for (1) lead to the corresponding results for the KdV and mKdV equations by setting β = 0 or α = 0 , respectively. Like the latter two equations, the Gardner equation is known to be completely integrable for both signs of β , but the solutions are quite different for β > 0 and β < 0 . The latter case is the hardest to deal with. The solitary wave solutions and solitons for the focusing Gardner equation follow from those of a focusing mKdV equation to which the Gardner equation can be reduced. The so-called “table-top” or “flat-top” soliton, corresponding to a large amplitude, only exists for the defocusing Gardner equation.
In some applications, the coefficients α and β are time dependent [35]. In that case, (1) has only a couple of conservation laws (conservation of mass and wave action flux, see Section 5) and is non-integrable. The treatment of (1) with varying coefficients is outside the scope of this paper.

3. Painlevé Analysis

In this Section we check if (1) has the Painlevé property. Performing the Painlevé test is straightforward but usually involves lengthy computations best performed with a software package such as PainleveTest.m [8]. In essence, one investigates a Laurent series solution for (1),
u ( x , t ) = g ( x , t ) σ k = 0 u k ( x , t ) g ( x , t ) k ,
in which the coefficients u k ( x , t ) are analytic functions in a neighborhood of the singular manifold g ( x , t ) . Furthermore, g ( x , t ) = 0 determines the poles and g ( x , t ) is assumed to be non-characteristic (i.e., g x ( x , t ) 0 ). The negative integer σ determines the leading order term, u 0 g σ , in (5). We summarize the main steps of the Painlevé test below; refer to [7] for details.
  • Step 1: Compute the leading order term. To determine σ and u 0 , substitute u 0 g σ into (1). Balance the minimal power terms in g, namely, g 3 σ 1 and g σ 3 , to get σ = 1 . Next, require that the leading terms (in g 4 ) vanish, yielding
    u 0 = ± 6 β g x if β < 0 , or
    u 0 = ± i 6 β g x if β > 0 .
  • Step 2: Compute the resonances. In this step, one determines which functions u k ( x , t ) in (5) will remain arbitrary. That happens at specific values of k, called resonances, denoted by r. To find the resonances, substitute u 0 g σ + u r g σ + r , with  σ = 1 and u 0 given in (6) or (7) into (1), and equate the coefficients of the dominant terms (in g r 4 ) that are linear in u r to get the characteristic equation
    ( r 4 ) ( r 3 ) ( r + 1 ) u r g x 3 = 0 .
    Since g x 0 , the resonances are r 1 = 1 , r 2 = 3 , and  r 3 = 4 . Resonance r 1 = 1 corresponds to the arbitrariness of g.
  • Step 3: Check the compatibility conditions. To do so, substitute
    u = 1 g k = 0 4 u k g k ,
    into (1) and verify that u 1 and u 2 can be unambiguously determined. Verify also that u 3 and u 4 are indeed arbitrary functions, meaning that no compatibility conditions arise. For (1), one readily determines the real functions,
    u 1 = α 2 β ± 1 2 6 β g x x g x ,
    u 2 = ± 4 β g x g t α 2 g x 2 6 β g x x 2 + 4 β g x g 3 x 4 β 6 β g x 3 ,
    when β < 0 (and complex expressions when β > 0 ). At resonance r = 3 , the compatibility condition (arising as the coefficient of g 3 ),
    u 0 2 α + 2 β u 1 6 ( u 0 ) x g x + u 0 g x x g x β u 0 2 ( u 0 ) x = 0 ,
    is identically satisfied upon substitution of (6) and (10). Likewise, a long but straightforward computation shows that the compatibility condition at r = 4 (appearing at order g 2 ),
    u 0 u 1 α + β u 1 + β u 0 2 u 2 + 3 ( u 0 ) x x g x + u 0 g t u 0 ( u 0 ) x α + 2 β u 1 β u 0 2 ( u 1 ) x + 3 ( u 0 ) x g x x + u 0 g 3 x = 0 ,
    is identically satisfied upon substitution of (6), (10), and (11).
Thus, at least in the neighborhood of g ( x , t ) = 0 , solution (9) is free of algebraic and logarithmic movable branch points. Apart from possible essential singularities (which the test is unable to detect), the movable singularities of its general solution are poles determined by g ( x , t ) = 0 . Note that (5) serves as a general solution because, as required by the Cauchy–Kovalevskaya theorem, (1) is of order 3 in x and that number agrees with 3 arbitrary functions g ( x , t ) , u 3 ( x , t ) , and u 4 ( x , t ) in (5).
Moreover, notice that truncating the Laurent series at the constant level term yields an auto-Bäcklund transformation,
u ( x , t ) = ± i 6 β g x g + u 1 ( x , t ) = ± i 6 β ln g x + u 1 ( x , t ) ,
because for an arbitrary g ( x , t ) , both u ( x , t ) and u 1 ( x , t ) must then be solutions of (1) with β > 0 . Setting u 1 ( x , t ) = 0 in (14) motivates the Hirota transformation discussed in Section 9, where it is shown that only for β > 0 one can find soliton solutions of (1).
In conclusion, the Gardner equation passes the Painlevé test and, therefore, has the Painlevé property, which is a good predictor that the equation has conservation laws, generalized symmetries, and so on. The investigation and actual computation of these quantities is based on a scaling symmetry of (1), which we will discuss next.

4. Scaling Symmetry

As stated, when β = 0 , the Gardner equation reduces to the KdV Equation (2), which is scaling homogeneous under the scaling (dilation) symmetry
( x , t , u ) ( κ 1 x , κ 3 t , κ 2 u ) = ( x ˜ , t ˜ , u ˜ ) ,
where κ 0 is an arbitrary scaling parameter. Indeed, replacing ( x , t , u ) in the KdV equation in terms of ( x ˜ , t ˜ , u ˜ ) yields
κ 5 ( u ˜ t ˜ + α u ˜ u ˜ x ˜ + u ˜ 3 x ˜ ) = 0 .
A fast way to compute (15) is to introduce the notions of weight, rank, and uniformity of rank. The weight, W, of a variable is the exponent of κ that multiplies that variable. Thus, based on (15), W ( x ) = 1 , W ( t ) = 3 , and  W ( u ) = 2 . Equivalently, W ( x ) = 1 and W ( t ) = 3 . The rank of a monomial is defined as the total weight of the monomial. For example, α u u x has rank 5 since the parameter α has no weight. An expression (or equation) is uniform in rank if all its monomial terms have equal rank.
Now, if (15) were not known yet, it could quickly be found by requiring that the KdV equation is uniform in rank, yielding
W ( u ) + W ( t ) = 2 W ( u ) + W ( x ) = W ( u ) + 3 W ( x ) .
Solving (17) gives W ( t ) = 3 W ( x ) and W ( u ) = 2 W ( x ) . Since κ is arbitrary, without loss of generality one can set W ( x ) = 1 , resulting in W ( t ) = 3 and W ( u ) = 2 . Thus, requiring uniformity in rank of a PDE allows one to compute the weights of the variables (and, hence, the scaling symmetry) with linear algebra.
Setting α = 0 , the Gardner equation becomes the mKdV equation. Requiring uniformity in rank readily yields W ( t ) = 3 and W ( u ) = 1 . Hence, the mKdV equation is scaling homogeneous under the transformation
( x , t , u ) ( κ 1 x , κ 3 t , κ u ) .
Because W ( u ) is different for the KdV and mKdV equations, the Gardner equation (1) will not be uniform in rank unless we use a trick. We will let the parameter α scale with some power of κ . Doing so, we must solve
W ( u ) + W ( t ) = W ( α ) + 2 W ( u ) + W ( x ) = 3 W ( u ) + W ( x ) = W ( u ) + 3 W ( x ) ,
yielding
W ( x ) = 1 , W ( t ) = 3 , W ( u ) = 1 , W ( α ) = 1 ,
which expresses that (1) has the scaling symmetry
( x , t , u , α ) ( κ 1 x , κ 3 t , κ u , κ α ) .
Starting with conservation laws, in what follows, we will use scaling symmetry to compute important quantities related to (1), thereby establishing its complete integrability.

5. Conservation Laws

A conservation law of (4) is an equation of the form
D t ρ + D x J = ˙ 0 ,
where the dot means that the equality should only hold on solutions u ( x , t ) of (4). ρ is called a conserved density (or charge), and J is the corresponding flux (or current).
The scalar functions ρ and J are functions of u and its partial derivatives with respect to x. All subsequent computations are performed on the jet space, which means that u and its partial derivatives with respect to x are treated as independent, alongside all monomials such as u 3 , u x 2 , etc. The density and flux could also explicitly depend on x and t, but we will not cover such exceptional cases.
In (22), D t is the total derivative with respect to t , defined by
D t ρ = ρ t + ρ [ u t ] = ρ t + k = 0 M ρ u k x D x k u t ,
where ρ [ u t ] is the Fréchet derivative of ρ in the direction of u t , and M is the highest order of ρ in x. In practice, one simply applies the chain rule for differentiation with respect to t, treating u , u x , u x x , etc., as independent functions.
Likewise, D x is the total derivative with respect to x,
D x J = J x + k = 0 N J u k x D x ( u k x ) = J x + k = 0 N J u k x u ( k + 1 ) x ,
where N is the order of J.
Since (22) is linear in the densities (and fluxes), a linear combination of densities with constant coefficients is still a density. The matching flux would, of course, be a linear combination of the corresponding fluxes with the same constant coefficients.
Returning to (1), one can readily verify that
D t ( u ) + D x ( 1 2 α u 2 + 1 3 β u 3 + u x x ) = 0 ,
D t ( u 2 ) + D x ( 2 3 α u 3 + 1 2 β u 4 u x 2 + 2 u u x x ) = 0 ,
D t ( α u 3 + 1 2 β u 4 3 u x 2 ) + D x 3 4 α 2 u 4 + α β u 5 + 1 3 β 2 u 6 6 α u u x 2 6 β u 2 u x 2 + 3 α u 2 u x x + 2 β u 3 u x x + 3 u x x 2 6 u x u 3 x = 0 .
Indeed, (25) is (1) written as a conservation law. Next, (26) straightforwardly follows after multiplication of (1) by 2 u and a bit of integration by parts. Clearly, (27) is far less obvious and will require a computational strategy [9,23] and the use of codes like InvariantsSymmetries.m [17] or ConservationLawsMD.m [18].
Notice that the above densities are uniform in rank. Indeed, ρ ( 1 ) = u has rank 1, ρ ( 2 ) = u 2 has rank 2, and ρ ( 3 ) is of rank 4. The corresponding fluxes are also uniform, with ranks 3 , 4 , and 6 . As a matter of fact, the entire conservation laws themselves are uniform in rank, with ranks 4 , 5 , and 7, respectively. This comes as no surprise because the defining Equation (22) is only non-trivial if evaluated on solutions of the PDE, and therefore the densities, fluxes, and conservation laws “inherit” (or adopt) the scaling homogeneity of the given PDE (and all its other continuous and discrete symmetries, for that matter).
It turns out that the list of conservation laws of (1) continues ad infinitum. The Gardner equation has infinitely many conservation laws, which is a clear indication that the PDE is completely integrable.
Using the scaling symmetry (21), we will now show how to compute (27), which is the shortest possible density of rank 4, since it is free of any terms that could be moved into the flux.
  • Step 1: Construct a candidate density of rank 4 as follows: Make a list of all monomials in u and α of rank 4 or less, i.e., { { u 4 , α u 3 , α 2 u 2 , α 3 u , α 4 } , { u 3 , α u 2 , α 2 u , α 3 } , { u 2 , α u , α 2 } , { u , α } , { 1 } } . For the construction of candidate densities, the constant terms α 4 , α 3 , , 1 can be removed. Then, for each monomial in that list, apply the correct number of x-derivatives so that the resulting term has exactly rank 4. The terms in the first sublist need no derivatives. Those in the second sublist need a single derivative. The next set of terms need two derivatives, etc. For example, for the first element in the third sublist, D x 2 u 2 = 2 u x 2 + 2 u u x x . Obviously, if we carry out partial integration, the highest derivative term u u x x only differs from u x 2 by the x-derivative of 1 2 u 2 and therefore can be ignored. Likewise, terms like u 2 u x , α u u x , α 2 u x , α u x x , and u 3 x can be neglected because they are x-derivatives of single-term monomials ( 1 3 u 3 , 1 2 α u 2 , etc.). There is no need to put terms like u u x x , u 2 u x in the density because they can be moved to the flux.
Gather the resulting monomials after stripping off numerical factors and removing scalar multiples of single-term densities of lower rank (with regard to (25) and (26) these are α 3 u and α 2 u 2 ). Finally, linearly combine the remaining terms with constant coefficients, yielding
ρ = c 1 u 4 + c 2 α u 3 + c 3 u x 2 ,
which is of first order in x ( M = 1 ).
  • Step 2: Compute the undetermined coefficients. Using (23), first compute
    D t ρ = ( 4 c 1 u 3 I + 3 α c 2 u 2 I + 2 c 3 u x D x ) [ u t ] = ( 4 c 1 u 3 + 3 α c 2 u 2 ) u t + 2 c 3 u x u x t ,
    where D x 0 = I is the identity operator. Replace u t and u x t = u t x from (1) to get
    E = ( 4 c 1 u 3 + 3 α c 2 u 2 ) ( α u u x + β u 2 u x + u 3 x ) + ( 2 c 3 u x ) ( α u u x + β u 2 u x + u 3 x ) x = α ( 4 c 1 + 3 β c 2 ) u 4 u x + 4 β c 1 u 5 u x + 4 c 1 u 3 u 3 x + 3 α 2 c 2 u 3 u x + 3 α c 2 u 2 u 3 x + 2 α c 3 u x 3 + 2 α c 3 u u x u x x + 4 β c 3 u u x 3 + 2 β c 3 u 2 u x u x x + 2 c 3 u x u 4 x .
    Next, find the constants c 1 , c 2 , and c 3 so that E = D t ρ matches D x J for some flux J (to be computed in Step 3 below). Mathematically, this means that E must be exact. The Euler operator (variational derivative),
    L u ( x ) = k = 0 K ( D x ) k u k x = u D x u x + D x 2 u x x D x 3 u 3 x + D x 4 u 4 x + ,
    allows one to test exactness [21,22,23,42], where K is the order of the expression the Euler operator is applied to. Therefore, K = 4 for E in (30). Consequently, (31) will terminate after five terms. E will be exact if L u ( x ) E 0 on the jet space (treating u , u x , u x x , etc., and also all monomials in such variables as independent). The computation of the terms in (31) involves nothing more than partial differentiations followed by (total) differentiations with respect to x. Of a total of 30 terms (not listed) generated, many terms are canceled, and one is left with
    L u ( x ) E = 4 ( 6 c 1 + β c 3 ) ( u x 3 + 3 u u x u x x ) 6 α ( 3 c 2 + c 3 ) u x u x x ,
    which must vanish identically, yielding the linear system 6 c 1 + β c 3 = 0 and 3 c 2 + c 3 = 0 with solution c 1 = 1 2 β , c 2 = 1 , c 3 = 3 . Substitute these constants into (28) to obtain
    ρ ( 3 ) = α u 3 + 1 2 β u 4 3 u x 2 ,
    the same expression as in (27). If one were only interested in the density, the computation would finish here. To continue with the computation of the flux (in the next step), substitute the constants into (30), yielding
    E = 5 α β u 4 u x + 2 β 2 u 5 u x + 2 β u 3 u 3 x + 3 α 2 u 3 u x + 3 α u 2 u 3 x 6 α u x 3 12 β u u x 3 6 α u u x u x x 6 β u 2 u x u x x 6 u x u 4 x .
  • Step 3: Compute the flux. Since D x J ( 3 ) = D t ρ ( 3 ) = E , to get J ( 3 ) one must integrate E with respect to x and reverse the sign. There is a tool from differential geometry, called the homotopy operator [43] p. 372, to carry out integration by parts on the jet space. As will be shown below, application of the homotopy operator reduces the integration on the jet space to a standard one-dimensional integration with respect to an auxiliary variable which will be denoted by λ .
The homotopy operator for variable u ( x ) , acting on a exact expression E of order K, is given [21,22,23,42,44] by
H u ( x ) E = 0 1 I u E [ λ u ] d λ λ ,
with integrand
I u E = k = 1 K i = 0 k 1 u i x ( D x ) k ( i + 1 ) E u k x = ( u I ) ( E u x ) + ( u x I u D x ) ( E u x x ) + ( u x x I u x D x + u D x 2 ) ( E u 3 x ) + ( u 3 x I u x x D x + u x D x 2 u D x 3 ) ( E u 4 x ) + .
In (35), ( I u E ) [ λ u ] means that once I u E is computed one must replace u by λ u , u x by λ u x , u x x by λ u x x , etc. Use (34), to get
I u E = ( u I ) ( 5 α β u 4 2 β 2 u 5 3 α 2 u 3 + 18 α u x 2 + 36 β u u x 2 + 6 α u u x x + 6 β u 2 u x x + 6 u 4 x ) + ( u x I u D x ) ( 6 α u u x + 6 β u 2 u x ) + ( u x x I u x D x + u D x 2 ) ( 3 α u 2 2 β u 3 ) + ( u 3 x I u x x D x + u x D x 2 u D x 3 ) ( 6 u x ) = 3 α 2 u 4 + 5 α β u 5 + 2 β 2 u 6 18 α u u x 2 24 β u 2 u x 2 + 9 α u 2 u x x + 8 β u 3 u x x + 6 u x x 2 12 u x u 3 x ,
which already has the terms of J ( 3 ) , but still with incorrect coefficients (and the opposite sign). Finally, using (35),
J ( 3 ) = H u ( x ) ( E ) = 0 1 ( I u ( E ) ) [ λ u ] d λ λ = 0 1 3 α 2 λ 3 u 4 + 5 α β λ 4 u 5 + 2 β 2 λ 5 u 6 18 α λ 2 u u x 2 24 β λ 3 u 2 u x 2 + 9 α λ 2 u 2 u x x + 8 β λ 3 u 3 u x x + 6 λ u x x 2 12 λ u x u 3 x d λ = 3 4 α 2 u 4 + α β u 5 + 1 3 β 2 u 6 6 α u u x 2 6 β u 2 u x 2 + 3 α u 2 u x x + 2 β u 3 u x x + 3 u x x 2 6 u x u 3 x ,
which is exactly the flux in (27).
We close this section with a remark about the use of conserved densities and constants of motion. If J vanishes at infinity (because u and its x-derivatives decay to zero), then integration of (22) with respect to x yields
d d t ρ d x = 0 .
Hence,
P = ρ d x
is constant in time and often referred to as a conserved quantity or constant of motion. Depending on the physical setting, the first few constants of motion express conservation of mass, momentum, and energy. Preserving these types of quantities plays an important role in testing the accuracy of numerical integrators. For example, some symmetric time-stepping methods [45] and explicit finite-difference schemes [46,47] preserve two of these three conserved quantities of low degree, while the more general symplectic integrators preserve the Hamiltonian structure [48,49]. The reader is referred to, e.g.,  [50,51,52,53] for an in-depth discussion of this subject.

6. Generalized Symmetries

A second criterion to establish the complete integrability of (1) is to show that the PDE has infinitely many generalized (higher-order) symmetries. As we will see, the computation of a hierarchy of such symmetries [10] is algorithmic and easier than for conservation laws.
A scalar function G ( x , t , u , u x , u x x , , u m x ) is called a generalized symmetry of (4) if and only if it leaves (4) invariant under the replacement u u + ϵ G within order ϵ , that is,
D t ( u + ϵ G ) = ˙ F ( u + ϵ G )
must hold up to order ϵ on the solutions of (4). Consequently, G must satisfy the linearized equation
D t G = G t + F [ G ] ,
where F is the Fréchet derivative of F in the direction of G:
F [ G ] = ϵ F ( u + ϵ G ) | ϵ = 0
= k = 0 n F u k x D x k G = F u I G + F u x D x G + F u x x D x 2 G + F u 3 x D x 3 G + ,
where n is the (highest) order of F. In (41) and (43), one must not only replace u by u + ϵ G , but also u x by u x + ϵ D x G ,   u x x by u x x + ϵ D x 2 G , etc. As before, I is the identity operator. The total derivative operators, D t and D x , were defined in (23) and (24), respectively. Since the defining Equation (42) is evaluated on solutions of the PDE, the higher-order symmetries inherit the scaling homogeneity of (4).
As we will show, (1) has infinitely many generalized symmetries, starting with
G ( 1 ) = u x ,
G ( 2 ) = α u u x + β u 2 u x + u 3 x ,
G ( 3 ) = 5 6 α 2 u 2 u x + 5 3 α β u 3 u x + 5 6 β 2 u 4 u x + 5 3 β u x 3 + 10 3 α u x u x x + 20 3 β u u x u x x + 5 3 α u u 3 x + 5 3 β u 2 u 3 x + u 5 x .
With regard to the weights in (20), the above symmetries are of ranks 2 , 4 , and 6, respectively.
Except for G ( 1 ) , each of these symmetries leads to a completely integrable nonlinear PDE, u t + G ( j ) = 0 , j = 2 , 3 , , where the weight of t increases as j increases. Notice that u t + G ( 2 ) = 0 corresponds to (1) itself.
Although the computation of symmetries is algorithmic [10] and technically simpler than for densities, it still requires the computation of a plethora of terms, and is best handled with symbolic software such as InvariantsSymmetries.m [10]. As an example, we show how to compute G ( 3 ) of rank 6.
  • Step 1: Construct a candidate symmetry of rank 6 as follows: List all monomials in u and α of rank 6 or less, i.e.,  { { u 6 , α u 5 , α 2 u 4 , α 3 u 3 , α 4 u 2 , α 5 u } , { u 5 , , α 4 u } , { u 4 , , α 3 u } , { u 3 , α u 2 , α 2 u } , { u 2 , α u } , { u } } , without the constant terms α 6 , α 5 , , α , 1 .
Then, apply the correct number of x-derivatives to the monomials in each of the six sublists so that the resulting terms have exactly rank 6. The elements in the first sublist need no derivatives. Those in the second sublist need a single derivative, etc. For example, for the first element in the fourth sublist, compute 3 u 3 x 3 = 6 u x 3 + 18 u u x u x x + 3 u 2 u 3 x . Do this for each element in all sublists and gather the resulting monomials after stripping off numerical factors. To avoid lower-rank symmetries being recomputed, remove scalar multiples of single-term symmetries of lower rank, as well as scalar multiples of the highest-derivative term in multiple-term symmetries of lower rank. Thus, with regard to (45) and (46), the monomials α 4 u x and α 2 u 3 x can be removed. Finally, linearly combine the remaining terms with constant coefficients, yielding
G = c 1 u 6 + α c 2 u 5 + α 2 c 3 u 4 + α 3 c 4 u 3 + α 4 c 5 u 2 + α 5 c 6 u + c 7 u 4 u x + α c 8 u 3 u x + α 2 c 9 u 2 u x + α 3 c 10 u u x + c 11 u 2 u x 2 + c 12 u 3 u x x + α c 13 u u x 2 + α c 14 u 2 u x x + α 2 c 15 u x 2 + α 2 c 16 u u x x + α 3 c 17 u x x + c 18 u x 3 + c 19 u u x u x x + c 20 u 2 u 3 x + α c 21 u x u x x + α c 22 u u 3 x + c 23 u x x 2 + c 24 u x u 3 x + c 25 u u 4 x + α c 26 u 4 x + c 27 u 5 x
which is of fifth order ( m = 5 ).
  • Step 2: Find the undetermined coefficients. Compute D t G and use (4) to remove u t , u x t , etc. This produces an expression with 176 terms (not listed). Next, compute (44) using
    F = ( α u u x + β u 2 u x + u 3 x )
    and G from (48), yielding 193 terms (not listed). Then, D t G F [ G ] , which has 138 terms, must vanish identically on the jet space, yielding a linear system of 20 equations for the nine nonzero coefficients ( c 7 , c 8 , c 9 , c 18 through c 22 , and  c 27 ) :
    2 ( 3 c 20 5 β c 27 ) = 0 , α ( 3 c 22 5 c 27 ) = 0 , , , α ( 3 c 22 5 c 27 ) = 0 , α ( 3 c 22 5 c 27 ) = 0 3 α ( 12 c 8 c 19 2 β c 21 4 β c 22 ) = 0 , 6 ( 3 c 18 + 2 c 19 + c 20 20 β c 27 ) = 0 .
    For brevity, we have shown only a couple of the shortest equations (coming from u u x u 5 x and u x u 5 x , respectively) and two of the longest equations (coming from u u x 2 u x x and u x u x x u 3 x , respectively). The 18 coefficients c 1 through c 6 , c 10 through c 17 , and c 23 through c 26 are all zero. Solving the system yields
    c 7 = 5 6 β 2 c 27 , c 8 = 5 3 β c 27 , c 9 = 5 6 c 27 , c 18 = 5 3 β c 27 , c 19 = 20 3 β c 27 , c 20 = 5 3 β c 27 , c 21 = 10 3 c 27 , c 22 = 5 3 c 27 .
    Set c 27 = 1 and substitute (51) into (48) to obtain
    G = 5 6 α 2 u 2 u x + 5 3 α β u 3 u x + 5 6 β 2 u 4 u x + 5 3 β u x 3 + 10 3 α u x u x x + 20 3 β u u x u x x + 5 3 α u u 3 x + 5 3 β u 2 u 3 x + u 5 x .
    which matches G ( 3 ) in (47).
The code InvariantsSymmetries.m can also be used to verify if an evolution equation, e.g., u t + G ( 3 ) = 0 , with  G ( 3 ) of rank 6, belongs to the hierarchy of completely integrable equations of a PDE of lower order. When asked to compute a symmetry of rank 4 for u t + G ( 3 ) = 0 , the software returns G ( 2 ) in (46), confirming that u t + G ( 3 ) = 0 corresponds to a generalized symmetry of u t + G ( 2 ) = 0 .

7. Recursion Operator

To prove that there are infinitely many higher-order symmetries, we will compute the recursion operator, which generates those symmetries sequentially, starting from the lowest order symmetry G ( 1 ) in (45).
As expected, the recursion operator for the Gardner equation is a combination of the well-known recursion operators for the KdV and mKdV equations (see p. 312 of [43] and [54]):
R = D x 2 + 2 3 α u + β u 2 I + 1 3 α u x D x 1 + 2 3 β u x D x 1 ( u I ) ,
where D x 1 denotes the anti-derivative (or integral) operator. R connects the symmetries sequentially (without gaps in this case; for examples with gaps, we refer to [55]):
G ( j + 1 ) = R G ( j ) , j = 1 , 2 , .
For example,
R u x = D x 2 + 2 3 α u + β u 2 I + 1 3 α u x D x 1 + 2 3 β u x D x 1 ( u I ) u x = u 3 x + 2 3 α u + β u 2 u x + 1 3 α u x D x 1 ( u x ) + 2 3 β u x D x 1 ( 1 2 u 2 ) x = α u u x + β u 2 u x + u 3 x ,
which is G ( 2 ) , and
R G ( 2 ) = D x 2 + 2 3 α u + β u 2 I + 1 3 α u x D x 1 + 2 3 β u x D x 1 ( u I ) α u u x + β u 2 u x + u 3 x = u 5 x + 5 3 α u u 3 x + 5 3 β u 2 u 3 x + + α 3 u x D x 1 1 2 α u 2 + 1 3 β u 3 + u x x x + 2 3 β u x D x 1 1 3 α u 3 + 1 4 β u 4 + u u x x 1 2 u x 2 x = 5 6 α 2 u 2 u x + 5 3 α β u 3 u x + 5 6 β 2 u 4 u x + 5 3 β u x 3 + 10 3 α u x u x x + 20 3 β u u x u x x + 5 3 α u u 3 x + 5 3 β u 2 u 3 x + u 5 x ,
which is G ( 3 ) .
Even for scalar equations like (1), the computation of recursion operators [55] is lengthy, but can be performed with specialized software such as PDERecursionOperator.m [16]. One can also use computer algebra to verify that R is a hereditary operator [56].
If R is a recursion operator for (4), then the Lie derivative [19,43,57,58] of R is zero, yielding
D t R + [ R , F ] = R t + R [ F ] + R F F R = ˙ 0 ,
where [ , ] and ∘ denote the commutator and composition of operators, respectively. F is the Fréchet derivative operator, defined as
F = k = 0 n F u k x D x k ,
where n is the order of F in (4). R [ F ] is the Fréchet derivative of R in the direction of F. For recursion operators of the form
R = j = 1 T U j ( u , u x , u x x , , u m 1 x ) S j ( D x , D x 2 , , D x 1 ) ( V j ( u , u x , u x x , , u m 2 x ) I ) ,
where T is the total number of terms in R , and  m 1 and m 2 are the orders of U and V, respectively, one has
R [ F ] = j = 0 m 1 ( D x j F ) U j u j x S j ( V j I ) + j = 0 m 2 U j S j ( D x j F ) V j u j x I .
With regard to (53), R is a linear integro-differential operator [59,60] which naturally splits into two pieces [55,58,61],
R = R 0 + R 1 ,
where R 0 is a local differential operator and R 1 is a non-local integral operator. R 0 is a linear combination of monomials involving D x , u, and parameters with weight (if applicable) so that each monomial has the correct rank. Note also that D x will always be “propagated” to the right in R 0 . For example,
D x 2 ( u I ) = D x u x I + u D x = u x x I + 2 u x D x + u D x 2 .
Based on the theory of recursion operators, R 1 is a linear combination with constant coefficients of terms of the form
G ( i ) D x 1 L u ( ρ ( j ) ) ,
where G ( i ) is a symmetry and L u ( ρ ( j ) ) is a cosymmetry (Euler operator applied to a density ρ ( j ) ) , selected such that R 1 has the correct rank [58,62]. To standardize the form of R 1 , one propagates D x to the left. For example, D x 1 u x D x = u x I D x 1 u x x I . The local and non-local operators are then added to obtain a candidate recursion operator.
We will now show how (53) is computed. Since the ranks G ( 1 ) and G ( 2 ) (as well as G ( 2 ) and G ( 3 ) ) differ by 2, R must have rank 2. Based on (20), W ( D x 1 ) = 1 , and one can readily verify that all the terms in (53) have rank 2.
  • Step 1: Compute the candidate recursion operator. Using (20), list all monomials in D x , u , and  α of rank 2 or less, i.e.,  { { D x 2 , u D x , α D x , u 2 I , α u I } , { D x , u I } } , where the trivial terms α I and α 2 I have been removed. Apply the correct number of x-derivatives to the monomials in each of the two sublists to assure that each term has rank 2. No action is needed on the first sublist. The elements in the second sublist need a single derivative, yielding { D x 2 , u x I , u D x } . After stripping off numerical factors and removing duplicates, linearly combine the resulting monomials with constant coefficients to obtain a candidate local operator
    R 0 = c 1 D x 2 + c 2 u D x + α c 3 D x + c 4 u 2 I + α c 5 u I + c 6 u x I .
    Using symmetry G ( 1 ) = u x and densities ρ ( 1 ) = u and ρ ( 2 ) = u 2 (all of low rank), compute L u ρ ( 1 ) = 1 and L u ρ ( 2 ) = 2 u . Then, with terms of type (63), make the candidate non-local operator
    R 1 = α c 7 u x D x 1 + c 8 u x D x 1 ( u I ) ,
    so that each term has rank 2. Add both pieces to get the candidate recursion operator
    R = c 1 D x 2 + ( c 2 u + α c 3 ) D x + ( c 4 u 2 + α c 5 u + c 6 u x ) I + α c 7 u x D x 1 + c 8 u x D x 1 ( u I ) .
  • Step 2: Compute the undetermined coefficients. Separately compute all the pieces in (57), beginning with
    D t R = c 2 u t D x + ( 2 c 4 u u t + α c 5 u t + c 6 u x t ) I + α c 7 u x t D x 1 + c 8 u x t D x 1 ( u I ) + c 8 u x D x 1 ( u t I ) = c 2 F D x + ( 2 c 4 u F + α c 5 F + c 6 D x F ) I + α c 7 ( D x F ) D x 1 + c 8 ( D x F ) D x 1 ( u I ) + c 8 u x D x 1 ( F I ) ,
    which can also be computed using (60). With F given in (49), insert
    D x F = F x = ( α u x 2 + α u u x x + 2 β u u x 2 + β u 2 u x x + u 4 x )
    into (67), expand, and simplify. This yields an expression of 27 terms (not listed). Some of the terms have D x and I ; others involve D x 1 , D x 1 ( u I ) , D x 1 ( u u x I ) , D x 1 ( u 2 u x I ) , and D x 1 ( u 3 x I ) . Next, using (58), compute
    F = D x 3 + ( α + β u ) u D x + ( α + 2 β u ) u x I .
    With (66) and (69), then compute
    R F = c 1 D x 2 + ( c 2 u + α c 3 ) D x + ( c 4 u 2 + α c 5 u + c 6 u x ) I + α c 7 u x D x 1 + c 8 u x D x 1 ( u I ) D x 3 + ( α + β u ) u D x + ( α + 2 β u ) u x I ,
    using, for example,
    D x 2 ( u 2 D x ) = D x ( 2 u u x D x + u 2 D x 2 ) = 2 u x 2 D x + 2 u u x x D x + 4 u u x D x 2 + u 2 D x 3 ,
    to consistently move the operator D x from the left to the right. To be able to match the integral terms in (67) which are mentioned below (68), repeated integration by parts is needed [55]. For example,
    D x 1 ( u D x 3 ) = u D x 2 D x 1 ( u x D x 2 ) = u D x 2 u x D x + D x 1 ( u x x D x ) = u D x 2 u x D x + u x x I D x 1 ( u 3 x I ) ,
    converts the integral D x 1 ( u D x 3 ) into D x 1 ( u 3 x I ) by moving the D x operator under the D x 1 operator from the right to the left. After integration by parts, R F has 58 terms (not listed). Next, compute
    F R = D x 3 + ( α + β u ) u D x + ( α + 2 β u ) u x I c 1 D x 2 + ( c 2 u + α c 3 ) D x + ( c 4 u 2 + α c 5 u + c 6 u x ) I + α c 7 u x D x 1 + c 8 u x D x 1 ( u I ) .
    To move the operator D x to the right, use, for example,
    D x 3 u x D x 1 ( u I ) = D x 2 u x x D x 1 ( u I ) + u u x I = D x u 3 x D x 1 ( u I ) + 2 u u x x I + u x 2 I + u u x D x = u 4 x D x 1 ( u I ) + 4 u x u x x I + 3 u u 3 x I + 3 u u x x D x + 2 u x 2 D x + u u x D x 2 ,
    and similar formulas (expressed in a more general form in [55]). The expanded expression for F R has 73 terms. Substitute the simplified expressions for (67), (70), and (73) into (57) and require that the resulting expression (which has 36 terms) vanishes identically, i.e., all monomials in u , u x , , I , D x , D x 2 , ⋯, should be treated as independent. This yields c 2 = c 3 = c 6 = 0 and a linear system,
    2 ( 3 c 4 2 β c 1 ) = 0 , α ( 3 c 5 2 c 1 ) = 0 , α ( 3 c 7 c 1 ) = 0 , 3 c 8 2 β c 1 = 0 , 3 α ( c 7 + c 5 c 1 ) = 0 , 3 ( c 8 + 2 c 4 2 β c 1 ) = 0 ,
    involving the remaining nonzero constants. Solve the system to get
    c 4 = 2 3 β c 1 , c 5 = 2 3 c 1 , c 7 = 1 3 c 1 , c 8 = 2 3 β c 1 .
    Finally, set c 1 = 1 and substitute (76) into (66), yielding
    R = D x 2 + 2 3 α u + β u 2 I + 1 3 α u x D x 1 + 2 3 β u x D x 1 ( u I ) ,
    which matches R in (53).

8. Lax Pair

Another way to prove the complete integrability of a nonlinear PDE of type (4) is to construct a Lax pair consisting of two linear PDEs in an auxiliary function whose compatibility requires that the original PDE is satisfied.
There are two flavors of Lax pairs. One is an operator formulation where the Lax pair consists of a pair of differential operators which leads to a higher order linear equation involving an auxiliary function. Alternatively, in the matrix formulation, the Lax pair is a set of two matrices satisfying a system of equations of first order in x and t, respectively. Only Lax pairs in operator form will be covered in this section. The reader is referred to the literature [1,2,63] for the matrix formalism (also called the zero curvature representation).
Finding a Lax pair in operator form is a nontrivial task and requires an educated guess about the order of the differential operators. But once the order is selected, one can take advantage of the scaling symmetry of the given PDE to construct a candidate for each operator because they inherit that scaling symmetry. Here again, the defining equation for a Lax pair is evaluated on solutions of the given PDE which supports that claim.
Using the KdV equation as prime example, Lax [64] showed that a completely integrable nonlinear PDE has an associated system of linear PDEs involving a pair of linear differential operators ( L , M ) and an auxiliary function ψ ( x , t ) ,
L ψ = λ ψ , ψ t = M ψ ,
where L and M are expressed in powers of D x with coefficients depending on u , u x , etc., and  ψ is an eigenfunction of L corresponding to eigenvalue λ . To guarantee complete integrability of (4), at least one non-trivial Lax pair ( L , M ) should exist, and the eigenvalue should not change in time which makes the problem isospectral.
We will show below that a one-parameter family of Lax pairs of (1) is given by
L = D x 2 + 2 ϵ u D x + 1 6 ( 6 ϵ 2 + β ) u 2 + α u + ( 6 ϵ ± 6 β ) u x I ,
M = 4 D x 3 12 ϵ u D x 2 ( 12 ϵ 2 + β ) u 2 + α u + ( 12 ϵ ± 6 β ) u x D x 1 3 ϵ ( 12 ϵ 2 + 2 β ) u 3 + 1 2 α ϵ u 2 + ( 12 ϵ 2 + β ± ϵ 6 β ) u u x + 1 2 α u x + 1 2 ( 6 ϵ ± 6 β ) u x x I ,
where ϵ is any real or complex number (not necessarily small). Substituting L and M into (78) yields
D x 2 ψ = 2 ϵ u D x ψ + λ ( ϵ 2 + 1 6 β ) u 2 1 6 α u ( ϵ ± 1 6 6 β ) u x ψ ,
D t ψ = 4 D x 3 ψ 12 ϵ u D x 2 ψ ( 12 ϵ 2 + β ) u 2 + α u + ( 12 ϵ ± 6 β ) u x D x ψ 1 3 ϵ ( 12 ϵ 2 + 2 β ) u 3 + 1 2 α ϵ u 2 + ( 12 ϵ 2 + β ± ϵ 6 β ) u u x + 1 2 α u x + 1 2 ( 6 ϵ ± 6 β ) u x x ψ .
The first equation is a Schrödinger-type equation for the arbitrary eigenfunction ψ with eigenvalue λ and potential ( ϵ 2 + 1 6 β ) u 2   + 1 6 α u   + ( ϵ ± 1 6 6 β ) u x . The second equation describes the time evolution of the eigenfunction.
A lengthy computation shows that the compatibility condition for (81) and (82) can be written as
D t D x 2 ψ D x 2 D t ψ = 1 6 ( 6 ϵ ± 6 β ) D x ( u t + α u u x + β u 2 u x + u 3 x ) + ( α + 2 ( 6 ϵ 2 + β ) u ) ( u t + α u u x + β u 2 u x + u 3 x ) ψ + 2 ϵ ( u t + α u u x + β u 2 u x + u 3 x ) D x ψ ,
where (81) and (82) were used repeatedly to eliminate D x t ψ , D t ψ , D x 5 ψ , D x 4 ψ , D x 3 ψ , and D x 2 ψ , in that order. From (83), it is clear that (81) and (82) will only be compatible on solutions of (1).
The compatibility of the equations in (78) may be expressed directly in terms of the operators L and M as follows
D t L ψ = L t ψ + L ( D t ψ ) = λ D t ψ ,
where
L t ψ D t L ψ L D t ψ .
With (78), one has
L t ψ + L M ψ = λ M ψ = M λ ψ = M L ψ .
For a non-trivial Lax pair for (4) to exist, (86) should vanish only on solutions of (4). If (88) is satisfied identically, i.e., without evaluation on solutions of the PDE, then the Lax operators L and M are considered trivial). Rearranging the terms yields
L t + L M M L ψ = L t + L , M ψ = ˙ 0 ,
or, expressed in operator form by suppressing the ψ ,
L t + L , M = ˙ O ,
where = ˙ means that equality holds only on solutions of the original PDE (4). As before, [ , ] is the commutator of the operators, and O is the zero operator. Equation (88) is called the Lax equation. We now show how the Lax pair ( L , M ) can be computed using the method discussed in [65].
  • Step 1: Construct a candidate for L and M . Since the KdV and mKdV equations are special cases of (1), it makes sense to search for L of rank 2 and M of rank 3. To construct L , list all monomials in D x , u , and  α of rank 2 or less, i.e.,  { { D x 2 , u D x , α D x , u 2 I , α u I } , { D x , u I } } , where the trivial terms α I and α 2 I have been removed. As explained in Section 7, apply D x to the elements in the second sublist and, as in (62), propagate D x to the right, yielding { D x 2 , u x I , u D x } . After removing duplicates, linearly combine the monomials from both sublists with constant coefficients to get a candidate for L :
    L = D x 2 + ( c 1 u + α c 2 ) D x + ( c 3 u 2 + α c 4 u + c 5 u x ) I ,
    where the coefficient of D x 2 has been set to one (for normalization).
To make a candidate for M , list all monomials in D x , u , and  α of rank 3 or less, i.e.,  { { D x 3 , u D x 2 , α D x 2 , u 2 D x , α u D x , α 2 D x , u x D x , u 3 I , α u 2 I , α 2 u I } , { D x 2 , u D x , α D x , u 2 I , α u I } , { D x , u I } } , where the trivial terms α I , α 2 I , and  α 3 I have been removed. Apply D x to the elements in the second sublist, yielding { D x 3 , u x D x , u D x 2 , α D x 2 , 2 u u x I , u 2 D x , α u x I , α u D x } . Apply D x 2 to the elements in the third sublist and use (62) to obtain { D x 3 , u x x I , 2 u x D x , u D x 2 } . After stripping off numerical factors and removing duplicates, linearly combine the resulting monomials with constant coefficients to obtain a candidate for M :
M = c 6 D x 3 + ( c 7 u + α c 8 ) D x 2 + ( c 9 u 2 + α c 10 u + α 2 c 11 + c 12 u x ) D x + ( c 13 u 3 + α c 14 u 2 + α 2 c 15 u + c 16 u u x + α c 17 u x + c 18 u x x ) I .
  • Step 2: Compute the undetermined coefficients. First, compute
    L t = c 1 u t D x + ( 2 c 3 u u t + α c 4 u t + c 5 u x t ) I = c 1 F D x + ( 2 c 3 u F + α c 4 F + c 5 D x F ) I ,
    which, after substituting F from (49) to replace u t and u x t , yields 14 terms. Next, compute
    L M = D x 2 + ( c 1 u + α c 2 ) D x + ( c 3 u 2 + α c 4 u + c 5 u x ) I c 6 D x 3 + ( c 7 u + α c 8 ) D x 2 + ( c 9 u 2 + α c 10 u + α 2 c 11 + c 12 u x ) D x + ( c 13 u 3 + α c 14 u 2 + α 2 c 15 u + c 16 u u x + α c 17 u x + c 18 u x x ) I ,
    which upon expansion has 125 terms, and
    M L = c 6 D x 3 + ( c 7 u + α c 8 ) D x 2 + ( c 9 u 2 + α c 10 u + α 2 c 11 + c 12 u x ) D x + ( c 13 u 3 + α c 14 u 2 + α 2 c 15 u + c 16 u u x + α c 17 u x + c 18 u x x ) I D x 2 + ( c 1 u + α c 2 ) D x + ( c 3 u 2 + α c 4 u + c 5 u x ) I ,
    which after expansion has 126 terms.
Substitute (91)–(93) into (88). The resulting expression (with 105 terms) should be identically equal to zero for any function ψ ( x , t ) . Set the coefficients of ψ , ψ x , and ψ x x equal to zero, and then set the coefficients of all monomials in u and its x-derivatives separately equal to zero. This yields a nonlinear system of 24 equations for the undetermined coefficients:
2 c 7 3 c 1 c 6 = 0 , 4 c 9 6 c 3 c 6 c 1 c 7 = 0 , c 7 3 c 1 c 6 = 0 , 4 c 9 6 c 3 c 6 c 1 c 7 = 0 , c 12 + 2 c 18 c 1 c 1 c 6 3 c 5 c 6 = 0 , α ( c 17 c 4 + c 2 c 18 c 4 c 6 c 5 c 8 ) = 0 .
For brevity, we have shown only a couple of the shortest equations (coming from u x D x 3 and u u x D x 2 , respectively), and two of the longest equations (coming from u 3 x D x and u 3 x I , respectively). Since each equation has a mixture of linear and and nonlinear terms, several solution branches occur. Mathematica’s Solve function returns five non-trivial solutions. Three of these solutions lead to Lax pairs of lower order or degenerate Lax pairs which will not be discussed. Instead, we focus on the two solutions that lead to Lax pairs that are useful in, e.g., the application of the inverse scattering transform (IST) and the Riemann–Hilbert methods to solve the Gardner equation. They have coefficients
c 2 = 6 c 4 1 3 c 1 , c 3 = 1 12 ( 3 c 1 2 + 2 β ) , c 5 = 1 6 ( 3 c 1 ± 6 β ) , c 6 = 4 , c 7 = 6 c 1 , c 9 = ( 3 c 1 2 + β ) , c 10 = c 1 c 8 1 , c 11 = ( 1 6 c 4 ) ( 1 6 c 4 c 1 c 8 ) 3 c 1 2 , c 12 = ( 6 c 1 ± 6 β ) , c 13 = 1 6 c 1 ( 3 c 1 2 + 2 β ) , c 14 = 3 c 1 2 + 2 β ( 1 6 c 4 ) c 1 c 8 ( 3 c 1 2 + 2 β ) 12 c 1 , c 15 = c 4 ( c 1 c 8 + 6 c 4 1 ) c 1 , c 16 = 1 2 ( 6 c 1 2 + 2 β ± c 1 6 β ) , c 17 = 3 c 1 ( c 1 c 8 1 ) ± ( c 1 c 8 + 6 c 4 1 ) 6 β 6 c 1 , c 18 = 1 2 ( 3 c 1 ± 6 β ) ,
where c 1 , c 4 , and c 8 are arbitrary constants. To be able to obtain the Lax pair for the KdV equation where c 1 = 0 , one should require that c 4 = 1 6 and c 8 = 0 , otherwise c 2 , c 11 , c 14 , and c 17 would become infinite. Notice that both requirements allow one to clear c 1 from all denominators. Furthermore, the coefficients then simplify into
c 2 = 0 , c 3 = 1 12 ( 3 c 1 2 + 2 β ) , c 4 = 1 6 , c 5 = 1 6 ( 3 c 1 ± 6 β ) , c 6 = 4 , c 7 = 6 c 1 , c 8 = 0 , c 9 = ( 3 c 1 2 + β ) , c 10 = 1 , c 11 = 0 , c 12 = ( 6 c 1 ± 6 β ) , c 13 = 1 6 c 1 ( 3 c 1 2 + 2 β ) , c 14 = 1 4 c 1 , c 15 = 0 , c 16 = 1 2 ( 6 c 1 2 + 2 β ± c 1 6 β ) , c 17 = 1 2 , c 18 = 1 2 ( 3 c 1 ± 6 β ) .
Finally, substitute the coefficients into (89) and (90) to get
L = D x 2 + c 1 u D x + 1 6 1 2 ( 3 c 1 2 + 2 β ) u 2 + α u + ( 3 c 1 ± 6 β ) u x I ,
M = 4 D x 3 6 c 1 u D x 2 ( 3 c 1 2 + β ) u 2 + α u + ( 6 c 1 ± 6 β ) u x D x 1 6 c 1 ( 3 c 1 2 + 2 β ) u 3 + 1 4 α c 1 u 2 + 1 2 ( 6 c 1 2 + 2 β ± c 1 6 β ) u u x + 1 2 α u x + 1 2 ( 3 c 1 ± 6 β ) u x x I ,
where the constant c 1 is arbitrary. Hence, this is a one-parameter family of Lax pairs. Set c 1 = 2 ϵ to get (79) and (80). With (97) and (98)
L t + [ L , M ] = 1 6 ( 3 c 1 ± 6 β ) D x ( u t + α u u x + β u 2 u x + u 3 x ) + ( α + ( 3 c 1 2 + 2 β ) u ) ( u t + α u u x + β u 2 u x + u 3 x ) I + c 1 ( u t + α u u x + β u 2 u x + u 3 x ) D x ,
which, after setting c 1 = 2 ϵ , is equivalent to (83).

9. Bilinear Form

In this section we show how the Gardner equation (1) can be transformed into the mKdV Equation (3) and how that helps with deriving Hirota’s bilinear representation [6] and, eventually, solitary wave and soliton solutions of (1). The existence of multi-soliton solutions (i.e., soliton solutions of any order) is yet another proof that the Gardner equation is completely integrable.
A simple shift of u allows one to remove the quadratic term α u u x . Indeed, set u = u ˜ α 2 β to replace (1) by
u ˜ t α 2 4 β u ˜ x + β u ˜ 2 u ˜ x + u ˜ 3 x = 0 ,
which is still in the original independent variables x and t. As will be shown below, the linear term in u ˜ x can also be removed by a change of independent variables.
  • Step 1: Construct a bilinear form as follows: first, integrate (100) with respect to x,
    t x u ˜ d x α 2 4 β u ˜ + 1 3 β u ˜ 3 + u ˜ x x = 0 ,
    setting the integration “constant” c ( t ) equal to zero. As with the mKdV equation [5], substitute the Hirota transformation
    u ˜ = i 6 β ln f + i g f i g x = 2 6 β Arctan f g x = 2 6 β f x g f g x f 2 + g 2
    into (101). Then, divide by 2 6 β and regroup the terms, yielding
    f 2 + g 2 g f t f g t + 3 f x g x x 3 g x f x x f g 3 x + g f 3 x α 2 4 β f g x g f x 6 g f x f g x ) ( f f x x + g g x x f x 2 g x 2 = 0 .
    Next, set the factors multiplying f 2 + g 2 and g f x f g x separately equal to zero, to obtain
    g f t f g t + 3 f x g x x 3 g x f x x f g 3 x + g f 3 x α 2 4 β f g x g f x = 0 ,
    f f x x + g g x x f x 2 g x 2 = 0 .
    Based on the scaling symmetry of (1) with weights (20) and the structure of the bilinear form of the mKdV equation [5], recast the above equations in bilinear form,
    ( c 1 D t + c 2 D x 3 + α 2 c 3 D x ) ( f · g ) = 0 ,
    c 4 D x 2 ( f · f + g · g ) = 0 ,
    with undetermined coefficients c 1 , c 2 , c 3 , and c 4 . Notice that with W ( f ) = W ( g ) = 0 , all terms in (104) and (106) have weights of three, whereas the terms in (105) and (107) have weights of two.
The bilinear operators D x and D t (not to be confused with total derivatives used in earlier sections) are defined as
D x m ( f · g ) = ( x x ) m f ( x , t ) g ( x , t ) | x = x = j = 0 m ( 1 ) m j m ! j ! ( m j ) ! j f x j m j g x m j ,
D t n ( f · g ) = ( t t ) n f ( x , t ) g ( x , t ) | t = t = j = 0 n ( 1 ) n j n ! j ! ( n j ) ! j f t j n j g t n j ,
which represent the Leibniz rule for x-derivatives (and t-derivatives, respectively) of products of functions with every other sign flipped. Explicitly,
( c 1 D t + c 2 D x 3 + α 2 c 3 D x ) ( f · g ) = c 1 ( g f t f g t ) + c 2 ( 3 f x g x x 3 g x f x x + g f 3 x f g 3 x ) + α 2 c 3 ( g f x f g x ) ,
c 4 D x 2 ( f · f + g · g ) = 2 c 4 ( f f x x + g g x x f x 2 g x 2 ) .
  • Step 2: Compute the undetermined coefficients. Equate (106) and (110) and treat all monomials in f and g and their derivatives as independent to get c 1 = c 2 = 1 , and c 3 = 1 4 β . Perform the same with (107) and (111) to get c 4 = 1 2 . Finally, substitute the constants into (106) and (107) and clear common factors to get a bilinear representation of (100):
    ( D t + D x 3 + α 2 4 β D x ) ( f · g ) = 0 , D x 2 ( f · f + g · g ) = 0 ,
    which, in light of (102), will only lead to real solutions for the focusing Gardner equation ( β > 0 ) .
The above bilinear formulation is expressed in the original variables x and t. Of course, the term α 2 4 β D x in (112) can be removed at the cost of introducing a new variable X. Indeed, using the chain rule for differentiation, one can readily verify that the Galilean transformation ( x , t , u ( x , t ) ) ( X , T , U ( X , T ) ) where
X = x + α 2 4 β t , T = t , u ( x , t ) = U ( X , T ) α 2 β ,
takes (1) into
U T + β U 2 U X + U 3 X = 0 ,
with bilinear representation [5]
( D T + D X 3 ) ( F · G ) = 0 , D X 2 ( F · F + G · G ) = 0 .
Once particular solutions F ( X , T ) and G ( X , T ) of
G F t F G t + 3 F x G x x 3 G x F x x F G 3 x + G F 3 x = 0 ,
F F x x + G G x x F x 2 G x 2 = 0 ,
are computed,
U ( X , T ) = 2 6 β F X G F G X F 2 + G 2
will solve (114). Solutions u ( x , t ) of (1) then follow from
u ( x , t ) = 2 6 β F X G F G X F 2 + G 2 α 2 β ,
after using (113) to return to the original variables u, x, and t.

10. Gardner Transformation

In this section we discuss a slight generalization of the Miura transformation [31], sometimes called the Gardner transformation or Gardner transform [32,66,67,68],
u = β γ v 2 ± 6 β v x + α β v ,
which connects solutions v ( x , t ) of the Gardner equation,
v t + α v v x + β v 2 v x + v 3 x = 0 ,
with β 0 , to solutions u ( x , t ) of the KdV equation,
u t + γ u u x + u 3 x = 0 ,
with arbitrary coefficient γ 0 . The use of γ will avoid confusion with α in (121) and, more importantly, allow us to set α = 0 to get the standard Miura transformation (see Section 12).
Substituting (120) into (122), it is straightforward to verify that
u t + γ u u x + u 3 x = ˙ β γ ( 2 v + α β ) I ± 6 β D x v t + α v v x + β v 2 v x + v 3 x ,
where = ˙ means that the left hand side is evaluated on solutions of (121). As before, I is the identity operator and D x is the total derivative operator defined in (24). Clearly, the Gardner transformation will only be real for the defocusing Gardner equation.
We will show how (120) can be computed using the scaling symmetries of the KdV and Gardner equations discussed in Section 4. Recall that W ( u ) = 2 , W ( v ) = W ( α ) = 1 , and  W ( γ ) = W ( β ) = 0 . Therefore, (120) is uniform in rank of rank two.
  • Step 1: Construct a candidate for the Gardner transformation. Make the list { { v 2 , α v , α 2 } , { v , α } } of monomials in v and α of rank 2 or less. By differentiating its elements, replace the second sublist by { v x } . Linearly combine the resulting elements with undetermined coefficients,
    u = c 1 v 2 + α c 2 v + α 2 c 3 + c 4 v x ,
    to generate a candidate for the Gardner transformation.
  • Step 2: Compute the undetermined coefficients. Substitute (124) into (122) and, using (121), replace v t and v t x . The resulting expression must vanish on the jet space, leading to c 2 c 3 = c 3 c 4 = 0 and half a dozen more complicated equations. For c 4 = 0 , (124) would become an algebraic transformation. Hence, c 3 = 0 and, after simplification of the more complicated equations, one is left with the following nonlinear system
    γ c 1 β = 0 , γ c 2 1 = 0 , 6 c 1 + γ c 4 2 = 0 , 3 γ c 1 c 2 β c 2 2 c 1 = 0 .
    Substitute the solution, c 1 = β γ , c 2 = 1 γ , c 4 = ± 1 γ 6 β , and c 3 = 0 , into (124) to get (120).
Of course, the  “uniformity-in-rank” argument also applies to (123), and therefore can be used to derive that equation. Observe that the ranks of the left and right hand sides of (123) match. Since the terms of the KdV equation (left) and Gardner equation (right) have ranks five and four, respectively, the operator that connects them must have rank one. Thus, only v I , α I , and  D x can occur in that operator. Apply the candidate for the operator,
( C 1 v + α C 2 ) I + C 3 D x ,
to (121), and equate the resulting expression to (122) after substitution of (120) but without evaluation on solutions of (121). After simplification, this yields the linear system
γ C 1 2 β = 0 , γ C 2 1 = 0 , γ C 3 ± 6 β = 0 , γ C 1 + β γ C 2 3 β = 0 .
Solve the system to obtain C 1 = 2 β γ , C 2 = 1 γ , and C 3 = ± 1 γ 6 β . Substitute the solution into (126) to get (123).
Using (120), some solutions for the KdV equation could be obtained from those of the Gardner equation (see, e.g., [11]). For α = 0 , (120) reduces to the Miura transformation, allowing one to generate solutions of the KdV equation from those of the mKdV equation. Although this is worthy of investigation, it is beyond the scope of this article. For an in-depth discussion of connections between the KdV, mKdV, and Gardner equations and additional applications of the Gardner and Miura transformations, we refer to [11,66,68].

11. The Korteweg–de Vries Equation

Since the KdV equation is a special case of the Gardner equation, its conservation laws, higher-order symmetries, and recursion operator immediately follow from those of (1) by setting β = 0 in (25)–(47), (45)–(47), and (53), respectively. The conservation laws for the KdV equation have been known since the 1970s [32,69] and have played an important role in the development of the concept of complete integrability. Its symmetries and recursion operator have been studied in, e.g.,  [54].
The well-known Lax pair [64] for the KdV equation,
L = D x 2 + 1 6 α u I , M = 4 D x 3 α u D x 1 2 α u x I ,
follows from (79) and (80) by setting β = ϵ = 0 .
The bilinear form for the KdV equation is much simpler than the one for the mKdV equation, and so are its soliton solutions. The interested reader is referred to [1,5,6,30] for the bilinear formulation and explicit formulas for the two-, three-, and  N-soliton solutions.
For completeness, we only include the solitary wave and cnoidal wave solutions, which were first derived in [28],
u ( x , t ) = 12 k 2 α sech 2 ( k x 4 k 3 t + δ ) = 12 k 2 α 1 tanh 2 ( k x 4 k 3 t + δ ) ,
u ( x , t ) = 8 k 2 α ( 1 m ) + 12 k 2 α m cn 2 ( k x 4 k 3 t + δ ; m ) ,
where m ( 0 , 1 ) is the modulus of the Jacobi elliptic cosine ( cn ) function. Both solutions are depicted in Figure 1. When m approaches 1, the peaks of the cn -squared solution become a bit taller, and the valleys become lower and flatter before they spread out horizontally to become the sech -squared solution. Both solutions (129) and (130) satisfy lim | x | u ( x , t ) = 0 for all t. The more general expressions corresponding to a nonzero boundary condition can be found in, e.g., refs. [1,2,29].
The soliton solutions of (2) can be computed with our code PDESolitonSolutions.m [24]. A discussion of their properties is outside the scope of this paper. Instead, we refer to [1,2,5,30] and the references given in [29].

12. The Modified Korteweg–de Vries Equation

Without extra work, we have the first three conservation laws [32] and higher-order symmetries [54] for the mKdV equation by setting α = 0 in (25)–(27) and (45)–(47), respectively. The recursion operator [54] connecting those symmetries follows from (53) with α = 0 .
Likewise, a one parameter Lax pair for the mKdV equation follows from (79) and (80) by setting α = 0 :
L = D x 2 + 2 ϵ u D x + 1 6 ( 6 ϵ 2 + β ) u 2 + ( 6 ϵ ± 6 β ) u x I ,
M = 4 D x 3 12 ϵ u D x 2 ( 12 ϵ 2 + β ) u 2 + ( 12 ϵ ± 6 β ) u x D x 1 3 ϵ ( 12 ϵ 2 + 2 β ) u 3 + ( 12 ϵ 2 + β ± ϵ 6 β ) u u x + 1 2 ( 6 ϵ ± 6 β ) u x x I ,
which for ϵ = 0 simplify into
L = D x 2 + 1 6 β u 2 ± 6 β u x I ,
M = 4 D x 3 ( β u 2 ± 6 β u x ) D x ( β u u x ± 1 2 6 β u x x ) I .
To our knowledge, this Lax pair first appeared in [70,71]. A further discussion of other Lax pairs for the mKdV equation can be found in [63,72]. In the latter paper, a more comprehensive algorithm to compute Lax pairs in operator form is presented, with the mKdV equation among its examples.
The Miura transformation [31,32,67],
u = β γ v 2 ± 6 β v x ,
which connects solutions of the mKdV Equation (3) (with dependent variable v ( x , t ) ) to solutions u ( x , t ) of the KdV Equation (122) readily follows from (120) by setting α = 0 . Likewise, (123) reduces to
u t + γ u u x + u 3 x = ˙ β γ 2 v I ± 6 β D x v t + β v 2 v x + v 3 x .
Various types of solutions for both the focusing and defocusing mKdV equations have been reported in the literature [14,15,34], including bright and dark solitons, breathers, rational solutions, kinks, etc. In this paper we mainly focus on solitary wave and soliton solutions of the mKdV equation. The nature of the solutions of (114) depends on the sign of β . Two cases have to be considered.
  • Case I: The focusing mKdV equation ( β > 0 ).
    The mKdV Equation (114) has solutions involving a cosh-function,
    U ( X , T ) = U 0 + 3 k 2 β U 0 ± U 0 2 + 3 k 2 2 β cosh Θ ,
    with Θ = k X ( β U 0 2 + k 2 ) k T + δ , where the boundary value U 0 , wave number k, and phase δ are arbitrary constants. To prevent blow-up in finite time we will only consider the plus sign in (137). The solution satisfying lim | X | U ( X , T ) = U 0 has been computed with Hirota’s method in [73], where only the case β > 0 was considered. Solution (137) is also valid for β < 0 with a caveat (see below). For U 0 = 0 the solution reduces to the well-known sech -solution,
    U ( X , T ) = ± 6 β k sech k X k 3 T + δ ,
    where the ± sign is due to the invariance of (114) under the discrete symmetry U U .
The soliton solutions of the focusing mKdV equation have been known for a long time and can be computed with a variety of methods. Adhering to Hirota’s method [5,6], the one-soliton solution readily follows from substitution of
F = e Θ = e k X ω T + δ and G = 1
into (116), yielding ω = k 3 and (117), which is identically satisfied. From (118), one then obtains
U ( X , T ) = 2 6 β F x 1 + F 2 = 2 6 β k e Θ 1 + e 2 Θ = 6 β k sech Θ = 6 β k sech ( k X k 3 T + δ ) = 2 6 β K sech 2 K X 8 K 3 T + δ ,
where K = k 2 , which matches (138). The two-soliton solution [5] follows from
F = e Θ 1 + e Θ 2 and G = 1 a 12 e Θ 1 + Θ 2 ,
with Θ i = k i X k i 3 T + δ i and a 12 = k 1 k 2 k 1 + k 2 2 . Then, from (118)
U ( X , T ) = 2 6 β k 1 e Θ 1 + k 2 e Θ 2 + a 12 ( k 1 e Θ 2 + k 2 e Θ 1 ) e Θ 1 + Θ 2 1 + e 2 Θ 1 + e 2 Θ 2 + 8 k 1 k 2 ( k 1 + k 2 ) 2 e Θ 1 + Θ 2 + a 12 2 e 2 Θ 1 + 2 Θ 2 .
Details of the derivation are given in [5] and [6], where formulas for the three- and N-soliton solutions can also be found.
  • Case II: The defocusing mKdV equation ( β < 0 ).
    Using the Zakharov–Shabat method, the defocusing mKdV equation was solved in [74]. Solution (137) with X replaced by X corresponds to the case N = 1 in [74]. Notice that when β < 0 , solution (137) will only exist if k 2 < 2 | β | 3 U 0 2 . This will be discussed in greater detail in the next section.
A new exact two-soliton solution of (114) with β < 0 is presented in [12]:
u ( X , T ) = 6 β sinh ( Φ + 2 δ ) e 2 Ψ + sinh ( Φ 2 δ ) e 2 Ψ + 2 sinh ( Φ ) ( 1 sinh 2 ( 2 δ ) ) sech ( 2 δ ) cosh ( Φ + 2 δ ) e 2 Ψ + cosh ( Φ 2 δ ) e 2 Ψ + 2 cosh ( Φ ) cosh ( 2 δ )
with Φ = X + 2 T , and Ψ = X + 2 ( 1 + 2 sech 2 ( 2 δ ) ) T + X 0 tanh ( 2 δ ) , with δ > 0 and X 0 arbitrary real parameters. This solution is obtained as the limit for the modulus going to one of a dark breather solution (involving Jacobi elliptic functions and elliptic integrals) of the defocusing mKdV equation.
Solutions (137) and (138) will be used in the next section to find table-top and hump-shaped solutions of (1). In turn, solution (143) will lead to a two-soliton solution of the defocusing Gardner equation.

13. Solitary Wave and Periodic Solutions

As pointed out in [37], the fact that (1) is invariant under the transformation
u u + α β
makes it possible to have solutions of different polarity for that equation, most notably, “bright” as well as “dark” solitons, depending on the initial conditions. However, a solution that vanishes at x = ± will be transformed by (144) into one that goes to α β as x ± . In this section we only cover a subset of the many types of solutions that (1) admits [75].
Depending on the sign of β in (1), it is straightforward [76] to find kink- and hump-type solitary waves solutions as well as periodic solutions in terms of the Jacobi elliptic sine and cosine functions. Indeed, using our symbolic code PDESpecialSolutions.m [26] for the tanh, sech , and Jacobi elliptic function methods, at a click of a button, one obtains simple exact solutions, which we cover first.
  • Case I: The focusing Gardner equation ( β > 0 ) .
    The most well-known solutions are
    u ( x , t ) = α 2 β ± 6 β k sech θ ,
    with θ = k x ( k 2 α 2 4 β ) k t + δ , and
    u ( x , t ) = α 2 β ± 6 β k m cn ( θ ; m ) ,
    with θ = k x + ( 1 2 m ) k 2 + α 2 4 β k t + δ . Both solutions for the minus sign (in front of the square root) are shown in Figure 2. Observe that as the value of m gets closer to 1, the  cn -function starts taking the shape of a sech -profile.
    Of course, (145) also follows directly from (138) upon application of the transformation (113). Using (137) in the same way yields
    u ( x , t ) = U 0 α 2 β + 3 k 2 β U 0 + U 0 2 + 3 k 2 2 β cosh θ ,
    with θ = k x ( k 2 + β U 0 2 α 2 4 β ) k t + δ and where U 0 is arbitrary. Setting U 0 = α 2 β yields
    u ( x , t ) = 6 k 2 α ( 1 + 1 + 6 β α 2 k 2 cosh θ ) ,
    with θ = k x k 3 t + δ = k ( x V t ) + δ , where V = k 2 denotes the speed of the wave and k 1 is the effective width of the solitary wave. This special solution is frequently used in the literature [35,66,76,77,78,79] and could also be computed via a Darboux transformation (see, e.g., [78]). For β > 0 , (148) is valid for all values of V, and since V > 0 the wave is travelling to the right.
  • Case II: The defocusing Gardner equation ( β < 0 ).
    The simplest exact solutions are
    u ( x , t ) = α 2 β ± 6 β k tanh θ ,
    with θ = k x + ( 2 k 2 + α 2 4 β ) k t + δ , and
    u ( x , t ) = α 2 β ± 6 β k m sn ( θ ; m ) ,
    with θ = k x + ( 1 + m ) k 2 + α 2 4 β k t + δ . In each of these solutions, δ is an arbitrary constant phase. The shock wave solution (149) and periodic solution (150) can be found in, e.g., ref. [74]. Both solutions for the minus sign (in front of the square root) are shown in Figure 3. As the value of m draws closer to 1, the  sn -function starts taking the shape of a tanh -profile.
    With respect to (148) where V = k 2 , note that when β < 0 , the argument of the square root, 1 + 6 β α 2 k 2 , will be zero when V = V crit = α 2 6 | β | . Thus, (148) is only valid when the speed is below that critical value ( V < V crit ). Thus, when dealing with the defocusing Gardner equation, if  V V crit , there is no solution of type (148). Turning the argument around, solutions for (148) of large amplitude (which is proportional to k 2 α ) or fast traveling waves (since V = k 2 ) can only occur if | β | is relatively small in comparison with α . For example, for  k = 2 , one must require that | β | < α 2 24 .
Figure 2. Graphs of (145) (left) and (146) (right) both with the minus signs in front of the square roots, and both for α = 12 , β = 6 , δ = 0 , and k = 2 . The curves on the right correspond to m = 0.25 (dashed line) and m = 0.9 (solid line).
Figure 2. Graphs of (145) (left) and (146) (right) both with the minus signs in front of the square roots, and both for α = 12 , β = 6 , δ = 0 , and k = 2 . The curves on the right correspond to m = 0.25 (dashed line) and m = 0.9 (solid line).
Mca 29 00091 g002
Figure 3. Graphs of (149) (left) and (150) (right), both for α = 12 , β = 6 , δ = 0 , and k = 2 . The curves on the right correspond to m = 0.25 (dashed line) and m = 0.9 (solid line).
Figure 3. Graphs of (149) (left) and (150) (right), both for α = 12 , β = 6 , δ = 0 , and k = 2 . The curves on the right correspond to m = 0.25 (dashed line) and m = 0.9 (solid line).
Mca 29 00091 g003
The graphs in Figure 4 are for α = 6 , β = 6 , and  δ = 0 (i.e., V < V crit is equivalent to k 2 < 1 ), for which (148) simplifies into
u ( x , t ) = k 2 1 + 1 k 2 cosh ( k x k 3 t ) .
As the value of k increases, the solitary wave get taller and narrower. At t = 0 and values of k very close to 1, the waves become flat at the top, hence the name table-top (or flat-top) waves. For the critical value k = 1 , the wave degenerates into a horizontal line corresponding to
lim k 1 k 2 1 + 1 k 2 cosh ( k x ) = 1 .
For values of k near 1, solution (151) can be very well approximated by a kink–antikink pair [79],
u ( x , t ) = 1 2 tanh k 2 ( x k 2 t + Δ ) tanh k 2 ( x k 2 t Δ ) ,
where Δ = 2 k 1 Arctanh ( 1 1 k 2 ) serves as a measure for the width of the table-top wave. These kink and anti-kink solutions, which correspond to the left and right flanks of the table-top solutions, are clearly visible in Figure 4 where k approaches 1. Although the expressions do not match analytically, Figure 5 shows that the graphs of (151) and (153) for t = 0 nearly overlap even when k is not close to 1.
Parenthetically, Grosse [80] (Equation (24)) computed a two-soliton solution of the defocusing mKdV equation. His analytic solution supposedly describes the interaction of two kink-solutions. It does not satisfy the equation exactly, but appears to be an excellent approximation to a solution, and therefore warrants further investigation.
Notice that all the above solutions follow the scaling homogeneity (21). Functions like cosh , tanh, cn , etc., have no weights. With regard to the weights (20), W ( α 2 β ) = 1 , as it should, because W ( u ) = 1 . Furthermore, W ( k ) = 1 because W ( x ) = 1 and W ( k x ) = 0 . All terms in any θ must have weight zero, in particular, W ( δ ) = 0 , and W ( m ) = 0 , where m is the modulus of any of the Jacobi elliptic functions. From W ( t ) = 3 , it follows that W ( V ) = 2 and W ( ω ) = 3 , where ω is the angular frequency in θ = k x ω t + δ = k ( x V t ) + δ . Hence, if  ω and V are polynomials in k, then ω can only have terms proportional to k 3 , α k 2 , and α 2 k , and V can only have terms in k 2 , α k , and  α 2 . The proportionality factors could have any powers of β since W ( β ) = 0 .

14. Soliton Solutions

When considering solitons we again must make the distinction between the focusing and defocusing Gardner equations.
  • Case I: The focusing Gardner equation ( β > 0 ).
    Using (113), solutions u ( x , t ) of the focusing version of (1) are given by
    u ( x , t ) = U ( X ( x , t ) , T ( t ) ) α 2 β ,
    where U ( X , T ) is any soliton solution of the focusing mKdV equation and X ( x , t ) = x + α 2 4 β t and T = t . For example, the two-soliton solution of (1) reads
    u ( x , t ) = α 2 β + 2 6 β k 1 e θ 1 + k 2 e θ 2 + a 12 ( k 1 e θ 2 + k 2 e θ 1 ) e θ 1 + θ 2 1 + e 2 θ 1 + e 2 θ 2 + 8 k 1 k 2 ( k 1 + k 2 ) 2 e θ 1 + θ 2 + a 12 2 e 2 θ 1 + 2 θ 2 ,
    where θ i = k i x ( k i 2 α 2 4 β ) k i t + δ i and a 12 = k 1 k 2 k 1 + k 2 2 . The elastic scattering of two solitons for the focusing Gardner equation is shown in Figure 6 and Figure 7, which have 2D and 3D graphs of (155) for α = β = 6 with k 1 = 3 2 , k 2 = 1 2 , and  δ 1 = δ 2 = 0 .
Figure 6. Graph of the two-soliton solution (155) of the focusing Gardner equation at three different moments in time.
Figure 6. Graph of the two-soliton solution (155) of the focusing Gardner equation at three different moments in time.
Mca 29 00091 g006
Figure 7. Bird’s eye view of a two-soliton collision for the focusing Gardner equation. Notice the phase shift after collision: the taller (faster) soliton is shifted forward and the shorter (slower) soliton backward relative to where they would have been if they had not collided.
Figure 7. Bird’s eye view of a two-soliton collision for the focusing Gardner equation. Notice the phase shift after collision: the taller (faster) soliton is shifted forward and the shorter (slower) soliton backward relative to where they would have been if they had not collided.
Mca 29 00091 g007
The Gardner equation has solitons of all orders N, which confirms once more that it is a completely integrable PDE.
In general, the existence of a two-soliton solution is no guarantee that the given PDE is completely integrable, but the existence of three-soliton solutions is a good indicator that the PDE is completely integrable. Indeed, as shown in [5], there are PDEs that have at most two-soliton solutions but no three-soliton solutions.
  • Case II: The defocusing Gardner equation ( β < 0 ).
    Based on (143), one obtains a two-soliton solution
    . u ( x , t ) = α 2 β + 6 β sinh ( θ + 2 δ ) e 2 η + sinh ( θ 2 δ ) e 2 η + 2 sinh ( θ ) ( 1 sinh 2 ( 2 δ ) ) sech ( 2 δ ) cosh ( θ + 2 δ ) e 2 η + cosh ( θ 2 δ ) e 2 η + 2 cosh ( θ ) cosh ( 2 δ ) ,
    with θ = x + 2 + α 2 4 β t , and η = x + 2 + α 2 4 β + 4 sech 2 ( 2 δ ) t + x 0 tanh ( 2 δ ) , with δ > 0 and x 0 arbitrary real parameters. Solution (156), which describes the coalescence of two wave fronts, is pictured in 2D and 3D in Figure 8 and Figure 9 for two different values of δ .

15. Symbolic Software

Using the concept of scaling homogeneity, we have been able to create powerful algorithms to investigate the complete integrability of systems of polynomial nonlinear PDEs. In this section, we give a brief summary of the available codes.
Our Mathematica code PainleveTest.m [8] automates the Painlevé test, which allows one to verify if a nonlinear PDE has the Painlevé property [7] as discussed in Section 3.
The Mathematica code InvariantsSymmetries.m [17] computes polynomial conserved densities and higher-order symmetries of nonlinear ( 1 + 1 ) -dimensional PDEs that can be written as a polynomial system of evolution equations. If a PDE has arbitrary parameters, the code allows one to derive conditions on these parameters so that the PDE admits conserved quantities or generalized symmetries. An example of such a “classification” problem is given in [9]. A discussion of the scope and limitations of the code can be found in [9,10].
To cover conservation laws of nonlinear PDEs in more than one space variable [23,44,81], we developed ConservationLawsMD.m [18], a Mathematica package to compute polynomial conservation laws of polynomial systems of nonlinear PDEs in space variables ( x , y , z ) and time t.
In [55], the authors show details of the algorithm to compute recursion operators for systems of nonlinear PDEs of type (4), including formulas for handling integro-differential operators used in Section 7. The Mathematica package PDERecursionOperator.m [16] performs the symbolic computation of recursion operators of systems of polynomial nonlinear evolution equations.
In Appendix B of [65], Larue presents LaxpairTester.m, a Mathematica code to verify Lax pairs in operator and matrix form. In [63], an algorithm is presented to compute Lax pairs in operator form, but that algorithm has not been implemented yet.
In addition, we developed the Mathematica package PDESpecialSolutions.m [25] to compute solitary wave solutions based on the tanh-method and generalizations for the sech , sn , and  cn functions [26].
Recently, we added the Mathematica code PDESolitonSolutions.m [24] to compute soliton solutions of polynomial PDEs based on a simplified version of Hirota’s method described in [5].
Since our codes only use tools from calculus, linear algebra, the calculus of variations, and differential geometry, these algorithms are fairly straightforward to implement in the syntax of computer algebra systems such as Mathematica, Maple, and REDUCE. Our software is open source and available in the public domain. All our Mathematica packages and notebooks are available on the Internet at https://people.mines.edu/whereman/ (accessed on 1 October 2024). A summary of the codes used in this paper can also be found at https://community.wolfram.com/groups/-/m/t/3275116 (accessed on 23 September 2024).

16. Conclusions and Future Work

The approach described in this paper and related software is applicable to large classes of nonlinear PDEs, which can be expressed as polynomial systems of evolution equations. As a prototypical example, we gave a detailed integrability analysis of the Gardner equation by computing its densities (and fluxes), higher-order symmetries, recursion operator, Lax pair, Hirota’s bilinear representation, and soliton solutions. The corresponding results for the KdV and mKdV equations were obtained by setting the coefficient of the cubic and quadratic term equal to zero, respectively.
We also showed how to compute the Gardner (Miura, resp.) transformation, which connects solutions of the Gardner (mKdV, resp.) equation to those of the KdV equation. With the Gardner transformation, some solutions of the KdV equation could be obtained from those of the Gardner equation shown in this paper. Likewise, applying the Miura transformation to solutions of the mKdV equation will lead to solutions of the KdV equation. When new solutions of the Gardner and mKdV equations are discovered, it would be worthwhile to investigate which solutions of the KdV equation they correspond to and, more importantly, if new solutions of the KdV equation could be computed that way (see, e.g., ref. [11]).
The crux of our computational strategy is a skillful use of the scaling symmetry of the PDE and relies on the observation that the defining equations for conservation laws, generalized symmetries, recursion operator, Lax pair, bilinear representation, and Gardner transformation should only hold on solutions of the given PDE. Consequently, the quantities (or operators) one computes inherit the scaling symmetry of the given PDE.
Since their defining equations are similar, it would also be possible to use this approach to compute symplectic and Hamiltonian (co-symplectic) operators of PDEs. In doing so, it would be possible to verify whether or not a PDE has a bi-Hamiltonian (or tri-Hamiltonian) structure, which is yet another criterion for its complete integrability. Further exploration of this idea, as well as the design of algorithms and codes for the computation of symplectic and Hamiltonian operators, is left for future work.
The methodology discussed in this paper might also apply to the Gardner equation in ( 2 + 1 ) dimensions [82,83], which is a combination of the Kadomtsev–Petviashvili (KP) and the modified KP equations.
The algorithms used in this paper are coded in Mathematica syntax, but can be adapted for major computer algebra systems such as Maple and REDUCE.

Author Contributions

Both authors contributed equally to this work. Conceptualization, W.H.; methodology, W.H.; software, Ü.G., validation, Ü.G., formal analysis, W.H. and Ü.G.; investigation, W.H. and Ü.G.; resources, W.H. and Ü.G.; writing—original draft preparation, W.H.; writing—review and editing, W.H. and Ü.G.; visualisation, Ü.G.; supervision, W.H.; project administration, W.H.; funding acquisition W.H. All authors have read and agreed to the published version of the manuscript.

Funding

The material presented in this paper is based in part upon research supported by the National Science Foundation (NSF) of the United States of America under Grants Nos. CCR-9300978, CCR-9625421, CCR-9901929, and CCF-0830783. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of NSF.

Data Availability Statement

Our software is open source and in the public domain. All our Mathematica packages and notebooks are available on the Internet at https://people.mines.edu/whereman/ (accessed on 2 September 2024).

Acknowledgments

The authors are grateful to Frank Verheest (Sterrenkundig Observatorium, Univ. Gent, Belgium) for valuable discussions on the use of the Gardner equation in plasma physics. We thank him, Rehana Naz (Dept. Math. Stat. Sci., Lahore School of Economics, Pakistan), Wen-Xiu Ma (Dept. Math. Stat., Univ. South Florida, Tampa, FL), and the anonymous referees for their thoughtful comments, suggestions, and references which helped further improve the paper.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Ablowitz, M.J.; Clarkson, P.A. Solitons, Nonlinear Evolution Equations and Inverse Scattering; London Mathematical Society Lecture Note Series; Cambridge University Press: Cambridge, UK, 1991; Volume 149. [Google Scholar]
  2. Ablowitz, M.J.; Segur, H. Solitons and The Inverse Scattering. In SIAM Studies in Applied Mathematics; SIAM: Philadelphia, PA, USA, 1981; Volume 4. [Google Scholar]
  3. Remoissenet, M. Waves Called Solitons: Concepts and Experiments, 3rd ed.; Springer: Berlin/Heidelberg, Germany, 1999. [Google Scholar]
  4. Zakharov, V.E. (Ed.) What is Integrability? Springer Series in Nonlinear Dynamics; Springer: Berlin/Heidelberg, Germany, 1991. [Google Scholar]
  5. Hereman, W.; Göktaş, Ü. Symbolic computation of solitary wave solutions and solitons through homogenization of degree. In Nonlinear and Modern Mathematical Physics, Proceedings of the 6th International Workshop on Nonlinear and Modern Mathematical Physics (NMMP2022), Tallahassee, FL, USA, 17–19 June 2022; Springer Proceedings in Mathematics & Statistics; Manukure, S., Ma, W.-X., Eds.; Springer: New York, NY, USA, 2024; Chapter 4; pp. 101–164. [Google Scholar]
  6. Hirota, R. The Direct Method in Soliton Theory. In Cambridge Tracts in Mathematics; Cambridge University Press: Cambridge, UK, 2004; Volume 155. [Google Scholar]
  7. Baldwin, D.; Hereman, W. Symbolic software for the Painlevé test of nonlinear ordinary and partial differential equations. J. Nonl. Math. Phys. 2006, 13, 90–110. [Google Scholar] [CrossRef]
  8. Baldwin, D.; Hereman, W. PainleveTest.m: A Mathematica Package for the Painlevé Test of Systems of Nonlinear Ordinary and Partial Differential Equations; Department of Applied Mathematics and Statistics, Colorado School of Mines: Golden, CO, USA, 2001. [Google Scholar]
  9. Göktaş, Ü.; Hereman, W. Symbolic computation of conserved densities for systems of nonlinear evolution equations. J. Symb. Comp. 1997, 24, 591–621. [Google Scholar] [CrossRef]
  10. Göktaş, Ü.; Hereman, W. Algorithmic computation of higher-order symmetries for nonlinear evolution and lattice equations. Adv. Comput. Math. 1999, 11, 55–80. [Google Scholar] [CrossRef]
  11. Kirane, M.; Stalin, S.; Arun, R.; Lakshmanan, M. Soliton molecules in Fermi–Pasta–Ulam–Tsingou lattice: Gardner equation approach. Chaos Solitons Fractals 2024, 178, 114393. [Google Scholar] [CrossRef]
  12. Mucalica, A.; Pelinovsky, D.E. Dark breathers on a snoidal wave background in the defocusing mKdV equation. Lett. Math. Phys. 2024, 14, 100. [Google Scholar] [CrossRef]
  13. Slyunyaev, A.V. Dynamics of localized waves with large amplitude in a weakly dispersive medium with a quadratic and positive cubic nonlinearity. J. Exper. Theor. Phys. 2001, 92, 529–534. [Google Scholar] [CrossRef]
  14. Slyunyaev, A.V.; Pelinovsky, E.N. Role of multiple soliton interactions in the generation of rogue waves: The modified Korteweg–de Vries framework. Phys. Rev. Lett. 2016, 117, 214501. [Google Scholar] [CrossRef]
  15. Zhang, G.; Yan, Z. Focusing and defocusing mKdV equations with nonzero boundary conditions: Inverse scattering transforms and soliton interactions. Phys. D 2014, 410, 132521. [Google Scholar] [CrossRef]
  16. Baldwin, D.; Hereman, W. PDERecursionOperator.m: A Mathematica Package for the Symbolic Computation of Recursion Operators for Nonlinear Partial Differential Equations; Department of Applied Mathematics and Statistics, Colorado School of Mines: Golden, CO, USA, 2003. [Google Scholar]
  17. Göktaş, Ü.; Hereman, W. InvariantsSymmetries.m: A Mathematica Integrability Package for the Computation of Invariants and Symmetries of Nonlinear Systems of Partial Differential Equations and Differential-Difference Equations; Department of Applied Mathematics and Statistics, Colorado School of Mines: Golden, CO, USA, 1997. [Google Scholar]
  18. Poole, D.; Hereman, W. ConservationLawsMD.m: A Mathematica Package for the Symbolic Computation of Conservation Laws of Polynomial Systems of Nonlinear PDEs in Multiple Space Dimensions; Department of Applied Mathematics and Statistics, Colorado School of Mines: Golden, CO, USA, 2009. [Google Scholar]
  19. Hereman, W.; Göktaş, Ü. Integrability Tests for Nonlinear Evolution Equations. In Computer Algebra Systems: A Practical Guide; Wester, M., Ed.; Wiley: New York, NY, USA, 1999; Chapter 12; pp. 211–232. [Google Scholar]
  20. Hereman, W.; Göktaş, Ü.; Colagrosso, M.; Miller, A. Algorithmic integrability tests for nonlinear differential and lattice equations. Comp. Phys. Commun. 1998, 115, 428–446. [Google Scholar] [CrossRef]
  21. Hereman, W.; Colagrosso, M.; Sayers, R.; Ringler, A.; Deconinck, B.; Nivala, M.; Hickman, M.S. Continuous and discrete homotopy operators with applications in integrability testing. In Differential Equations with Symbolic Computation; Wang, D., Zheng, Z., Eds.; Birkhäuser Verlag: Basel, Switzerland, 2005; Chapter 15; pp. 255–290. [Google Scholar]
  22. Hereman, W.; Adams, P.J.; Eklund, H.L.; Hickman, M.S.; Herbst, B.M. Direct methods and symbolic software for conservation laws of nonlinear equations. In Advances in Nonlinear Waves and Symbolic Computation; Yan, Z., Ed.; Nova Science Publishers: New York, NY, USA, 2009; Chapter 2; pp. 19–79. [Google Scholar]
  23. Poole, D.; Hereman, W. Symbolic computation of conservation laws for nonlinear partial differential equations in multiple space dimensions. J. Symb. Comp. 2011, 46, 1355–1377. [Google Scholar] [CrossRef]
  24. Göktaş, Ü.; Hereman, W. PDESolitonSolutions.m: A Mathematica Package for the Symbolic Computation of Solitary Wave and Soliton Solutions of Polynomial Nonlinear PDEs Using a Simplified Version of Hirota’s Method; Department of Applied Mathematics and Statistics, Colorado School of Mines: Golden, CO, USA, 2023. [Google Scholar]
  25. Baldwin, D.; Hereman, W. PDESpecialSolutions.m: A Mathematica Package for the Symbolic Computation of Exact Solutions Expressible in Hyperbolic and Elliptic Functions for Systems of Nonlinear Partial Differential Equations; Department of Applied Mathematics and Statistics, Colorado School of Mines: Golden, CO, USA, 2003. [Google Scholar]
  26. Baldwin, D.; Göktaş, Ü.; Hereman, W.; Hong, L.; Martino, R.S.; Miller, J.C. Symbolic computation of exact solutions expressible in hyperbolic and elliptic functions for nonlinear PDEs. J. Sym. Comput. 2004, 37, 669–705. [Google Scholar] [CrossRef]
  27. Gardner, C.S.; Greene, J.M.; Kruskal, M.D.; Miura, R.M. Korteweg-de Vries equation and generalizations. VI. Methods for exact solution. Commun. Pure Appl. Math. 1974, 27, 97–133. [Google Scholar] [CrossRef]
  28. Korteweg, D.J.; de Vries, G. XLI. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves. Philos. Mag. 1895, 39, 422–443. [Google Scholar] [CrossRef]
  29. Hereman, W. Shallow water waves and solitary waves. In Encyclopedia of Complexity and Systems Science; Meyers, R.A., Ed.; Springer: Berlin/Heidelberg, Germany, 2009; Volume 11, pp. 8112–8125. [Google Scholar]
  30. Drazin, P.G.; Johnson, R.S. Solitons: An Introduction. In Cambridge Texts in Applied Mathematics; Cambridge University Press: Cambridge, UK, 1989. [Google Scholar]
  31. Miura, R. Korteweg-de Vries equation and generalizations. I. A remarkable explicit nonlinear transformation. J. Math. Phys. 1968, 9, 1202–1204. [Google Scholar] [CrossRef]
  32. Miura, R.M.; Gardner, C.S.; Kruskal, M.D. Korteweg-de Vries equation and generalizations. II. Existence of conservation laws and constants of motion. J. Math. Phys. 1968, 9, 1204–1209. [Google Scholar] [CrossRef]
  33. Kakutani, T.; Yamasaki, N. Solitary waves on a two-Layer fluid. J. Phys. Soc. Jpn. 1978, 45, 674–679. [Google Scholar] [CrossRef]
  34. Zhang, D.-J.; Zhao, S.-L.; Sun, Y.Y.; Zhou, J. Solutions to the modified Korteweg–de Vries equation. Rev. Math. Phys. 2014, 26, 1430006. [Google Scholar] [CrossRef]
  35. Grimshaw, R. Nonlinear wave equations for oceanic internal solitary waves. Stud. Appl. Math. 2015, 136, 214–237. [Google Scholar] [CrossRef]
  36. Miles, J.W. On internal solitary waves. Tellus A Dyn. Meteor. Ocean. 1979, 31, 456–462. [Google Scholar] [CrossRef]
  37. Kamchatnov, A.M.; Kuo, Y.-H.; Lin, T.-C.; Horng, T.-L.; Gou, S.-C.; Clift, R.; El, G.A.; Grimshaw, R.H.J. Undular bore theory for the Gardner equation. Phys. Rev. E 2012, 86, 036605. [Google Scholar] [CrossRef]
  38. Olivier, C.P.; Verheest, F. Overtaking collisions of double layers and solitons: Tripolar structures and dynamical polarity switches. Phys. Plasmas 2020, 27, 062303. [Google Scholar] [CrossRef]
  39. Olivier, C.P.; Verheest, F.; Maharaj, S.K. A small-amplitude study of solitons near critical plasma compositions. J. Plasma Phys. 2016, 82, 905820605. [Google Scholar] [CrossRef]
  40. Torvén, S. Modified Korteweg-de Vries equation for propagating double layers in plasmas. Phys. Rev. Lett. 1981, 47, 1053–1056. [Google Scholar] [CrossRef]
  41. Verheest, F.; Hereman, W.; Olivier, C.P. On the Gardner equation for nonlinear waves in multispecies plasmas. J. Plasma Phys. 2024. in preparation. [Google Scholar] [CrossRef]
  42. Hereman, W.; Deconinck, B.; Poole, L.D. Continuous and discrete homotopy operators: A theoretical approach made concrete. Math. Comput. Simul. 2007, 74, 352–360. [Google Scholar] [CrossRef]
  43. Olver, P.J. Applications of Lie Groups to Differential Equations, 2nd ed.; Graduate Texts in Mathematics;, Springer: New York, NY, USA, 1993; Volume 107. [Google Scholar]
  44. Poole, D.; Hereman, W. The homotopy operator method for symbolic integration by parts and inversion of divergences with applications. Appl. Anal. 2010, 89, 433–455. [Google Scholar] [CrossRef]
  45. Cohen, D.; Hairer, E.; Lubich, C. Conservation of energy, momentum and actions in numerical discretizations of non-linear wave equations. Numer. Math. 2008, 110, 113–143. [Google Scholar] [CrossRef]
  46. Frasca-Caccia, G.; Hydon, P.E. Locally conservative finite difference schemes for the modified KdV equation. J. Comput. Dyn. 2019, 6, 307–323. [Google Scholar] [CrossRef]
  47. Sanz-Serna, J.M. An explicit finite-difference scheme with exact conservation properties. J. Comput. Phys. 1982, 47, 199–210. [Google Scholar] [CrossRef]
  48. Ascher, U.M.; McLachlan, R.I. On symplectic and multisymplectic schemes for the KdV equation. J. Sci. Comput. 2005, 25, 83–104. [Google Scholar]
  49. Bridges, T.J.; Reich, S. Numerical methods for Hamiltonian PDEs. J. Phys. A Math. Gen. 2006, 39, 5287–5320. [Google Scholar] [CrossRef]
  50. Quispel, R.; McLachlan, R. Preface to special issue on ‘Geometric Numerical Integration of Differential Equations’. J. Phys. A Math. Gen. 2006, 39, E01–E03. [Google Scholar]
  51. Hairer, E.; Lubich, C.; Wanner, G. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd ed.; Springer Series in Computational Mathematics; Springer: Berlin/Heidelberg, Germany, 2006; Volume 31. [Google Scholar]
  52. LeVeque, R.J. Numerical Methods for Conservation Laws, 2nd ed.; Lectures in Mathematics; Birkhäuser: Basel, Switzerland, 1992; Volume 214. [Google Scholar]
  53. Sanz-Serna, J.M.; Calvo, M.P. Numerical Hamiltonian Problems. Applied Mathematics and Mathematical Computation; Courier Dover Publications: Mineola, NY, USA, 2018. [Google Scholar]
  54. Olver, P. Evolution equations possessing infinitely many symmetries. J. Math. Phys. 1977, 18, 1212–1215. [Google Scholar] [CrossRef]
  55. Baldwin, D.; Hereman, W. A symbolic algorithm for computing recursion operators of nonlinear partial differential equations. Int. J. Computer Math. 2010, 97, 1094–1119. [Google Scholar] [CrossRef]
  56. Fuchssteiner, B.; Oevel, W.; Wiwianka, W. Computer-algebra methods for investigation of hereditary operators of higher order soliton equations. Comput. Phys. Commun. 1987, 44, 47–55. [Google Scholar] [CrossRef]
  57. Wang, J.P. Symmetries and Conservation Laws of Evolution Equations. Ph.D. Thesis, Thomas Stieltjes Institute for Mathematics, Free University of Amsterdam, Amsterdam, The Netherlands, 1998. [Google Scholar]
  58. Wang, J.P. A list of (1+1) dimensional integrable equations and their properties. J. Nonl. Math. Phys. 2002, 9, 213–233. [Google Scholar] [CrossRef]
  59. Sanders, J.A.; Wang, J.P. Integrable systems and their recursion operators. Nonlinear Anal. 2001, 47, 5213–5240. [Google Scholar] [CrossRef]
  60. Sanders, J.A.; Wang, J.P. On recursion operators. Phys. D 2001, 149, 1–10. [Google Scholar] [CrossRef]
  61. Baldwin, D.; Hereman, W.; Sayers, J. Symbolic algorithms for the Painlevé test, special solutions, and recursion operators for nonlinear PDEs. In Group Theory and Numerical Analysis; Winternitz, P., Gomez-Ullate, D., Iserles, A., Levi, D., Olver, P.J., Quispel, R., Tempesta, P., Eds.; CRM Proceedings and Lecture Notes; AMS: Providence, RI, USA, 2005; Volume 39, pp. 17–32. [Google Scholar]
  62. Bilge, A.H. On the equivalence of linearization and formal symmetries as integrability tests for evolution equations. J. Phys. A: Math. Gen. 1973, 26, 7511–7519. [Google Scholar] [CrossRef]
  63. Hickman, M.; Hereman, W.; Larue, J.; Göktaş, Ü. Scaling invariant Lax pairs of nonlinear evolution equations. Appl. Anal. 2012, 91, 381–402. [Google Scholar] [CrossRef]
  64. Lax, P.D. Integrals of nonlinear equations of evolution and solitary waves. Commun. Pure Appl. Math. 1968, 21, 467–490. [Google Scholar] [CrossRef]
  65. Larue, J. Symbolic Verification of Operator and Matrix Lax Pairs for Some Completely Integrable Nonlinear Partial Differential Equations. Master’s Thesis, Department of Mathematical and Computer Sciences, Colorado School of Mines, Golden, CO, USA, 2011. [Google Scholar]
  66. Alejo, M.A.; Muñoz, C.; Vega, L. The Gardner equation and the L2-stability of the N-soliton solution of the Korteweg–de Vries equation. Trans. Amer. Math. Soc. 2012, 365, 195–212. [Google Scholar] [CrossRef]
  67. Miura, M.R. The Korteweg de Vries equation: A survey of results. SIAM Rev. 1976, 18, 412–459. [Google Scholar] [CrossRef]
  68. Muñoz, C. The Gardner equation and the stability of multi-kink solutions of the mKdV equation. Discrete Cont. Dyn. Sys. 2016, 36, 3811–3843. [Google Scholar] [CrossRef]
  69. Kruskal, M.D.; Miura, R.M.; Gardner, C.S.; Zabusky, N.J. Korteweg-de Vries equation and generalizations. V. Uniqueness and nonexistence of polynomial conservation laws. J. Math. Phys. 1970, 11, 952–960. [Google Scholar] [CrossRef]
  70. Wadati, M. The exact solution of the modified Korteweg-de Vries equation. J. Phys. Soc. Jpn. 1972, 32, 1681. [Google Scholar] [CrossRef]
  71. Wadati, M. The modified Korteweg-de Vries equation. J. Phys. Soc. Jpn. 1973, 34, 1289–1296. [Google Scholar] [CrossRef]
  72. Burde, G.I. Lax pairs for the modified KdV equation. Axioms 2024, 13, 121. [Google Scholar] [CrossRef]
  73. Ablowitz, M.J.; Satsuma, J. Solitons and rational solutions of nonlinear evolution equations. J. Math. Phys. 1978, 19, 2180–2186. [Google Scholar] [CrossRef]
  74. Ono, H. Solitons on a background and a shock wave. J. Phys. Soc. Jpn. 1976, 40, 1487–1497. [Google Scholar] [CrossRef]
  75. Pelinovsky, E.; Talipova, T.; Didenkulova, E. Rational solitons in the Gardner-like models. Fluids 2022, 7, 294. [Google Scholar] [CrossRef]
  76. Wazwaz, A.-M. New solitons and kink solutions for the Gardner equation. Commun. Nonl. Sci. Numer. Simulat. 2007, 12, 1395–1404. [Google Scholar] [CrossRef]
  77. Grimshaw, R.; Pelinovsky, D.; Pelinovsky, E.; Slunyaev, A. Generation of large-amplitude solitons in the extended Korteweg-de Vries equation. Chaos 2002, 12, 1070–1076. [Google Scholar] [CrossRef] [PubMed]
  78. Slyunyaev, A.V.; Pelinovskiĭ, E.N. Dynamics of large-amplitude solitons. J. Exper. Theor. Phys. 1999, 89, 173–181. [Google Scholar] [CrossRef]
  79. Pelinovskiĭ, E.N.; Slyunyaev, A.V. Generation and interaction of large-amplitude solitons. J. Exper. Theor. Phys. Lett. 1998, 67, 655–661. [Google Scholar] [CrossRef]
  80. Grosse, H. Solitons of the modified KdV equation. Lett. Math. Phys. 1984, 8, 313–319. [Google Scholar] [CrossRef]
  81. Hereman, W. Symbolic computation of conservation laws of nonlinear partial differential equations in multi-dimensions. Int. J. Quant. Chem. 2006, 106, 278–299. [Google Scholar] [CrossRef]
  82. Konopelchenko, B.G. Introduction to Multidimensional Integrable Equations. The Inverse Spectral Transform in (2+1) Dimensions; Plenum Press: New York, NY, USA, 1992. [Google Scholar]
  83. Zhou, R.; Ma, W.-X. Algebro-geometric solutions of the (2+1)-dimensional Gardner equation. Nuovo Cimento B 2000, 115, 1419–1431. [Google Scholar]
Figure 1. Graphs of the solitary wave (solid line) and cnoidal wave (dashed line) solutions for α = 6 , k = 2 , m = 9 10 , and δ = 0 .
Figure 1. Graphs of the solitary wave (solid line) and cnoidal wave (dashed line) solutions for α = 6 , k = 2 , m = 9 10 , and δ = 0 .
Mca 29 00091 g001
Figure 4. Graphs of (151) for k = 0.5 , k = 0.625 , and k = 0.75 (left), and k = 0.9999 , k = 0.999995 , k = 0.9999995 , and k = 1 (right). The solitary wave becomes taller and narrower as the value of k increases.
Figure 4. Graphs of (151) for k = 0.5 , k = 0.625 , and k = 0.75 (left), and k = 0.9999 , k = 0.999995 , k = 0.9999995 , and k = 1 (right). The solitary wave becomes taller and narrower as the value of k increases.
Mca 29 00091 g004
Figure 5. Graphs of (151) (full line) and (153) (dashed line) for four different values of k.
Figure 5. Graphs of (151) (full line) and (153) (dashed line) for four different values of k.
Mca 29 00091 g005
Figure 8. 2D and 3D graphs of solution (156) for α = 6 , β = 6 , δ = 0.65 and x 0 = 0 .
Figure 8. 2D and 3D graphs of solution (156) for α = 6 , β = 6 , δ = 0.65 and x 0 = 0 .
Mca 29 00091 g008
Figure 9. 2D and 3D graphs of solution (156) for α = 6 , β = 6 , δ = 3 and x 0 = 0 .
Figure 9. 2D and 3D graphs of solution (156) for α = 6 , β = 6 , δ = 3 and x 0 = 0 .
Mca 29 00091 g009
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Hereman, W.; Göktaş, Ü. Using Symmetries to Investigate the Complete Integrability, Solitary Wave Solutions and Solitons of the Gardner Equation. Math. Comput. Appl. 2024, 29, 91. https://doi.org/10.3390/mca29050091

AMA Style

Hereman W, Göktaş Ü. Using Symmetries to Investigate the Complete Integrability, Solitary Wave Solutions and Solitons of the Gardner Equation. Mathematical and Computational Applications. 2024; 29(5):91. https://doi.org/10.3390/mca29050091

Chicago/Turabian Style

Hereman, Willy, and Ünal Göktaş. 2024. "Using Symmetries to Investigate the Complete Integrability, Solitary Wave Solutions and Solitons of the Gardner Equation" Mathematical and Computational Applications 29, no. 5: 91. https://doi.org/10.3390/mca29050091

APA Style

Hereman, W., & Göktaş, Ü. (2024). Using Symmetries to Investigate the Complete Integrability, Solitary Wave Solutions and Solitons of the Gardner Equation. Mathematical and Computational Applications, 29(5), 91. https://doi.org/10.3390/mca29050091

Article Metrics

Back to TopTop