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 dimensions, i.e., one space variable in addition to time . 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
where
is the dependent variable (or field) which is a function of space variable
x and time
t. The subscripts denote partial derivatives, i.e.,
and
The parameters
and
are real numbers for which the values
and
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
, (
1) reduces to the ubiquitous Korteweg–de Vries equation [
28],
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
, (
1) becomes the mKdV equation [
1,
27],
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
in (
1), because if
, replacing
u by
would make the coefficient of
positive again. However, no discrete symmetries of
or
u will flip the sign of the coefficient of
, so the cases for positive and negative
will have to be treated separately. Exact solutions of (
1) with
, called the
focusing Gardner equation, are quite different from these of the
defocusing case, where
.
We focus on (
1) because it is a typical example of a scalar
-dimensional evolution equation,
of order
n in
x and first order in
t, with
in (
1). More importantly, many of the results for (
1) lead to the corresponding results for the KdV and mKdV equations by setting
or
, 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
and
. 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),
in which the coefficients
are analytic functions in a neighborhood of the singular manifold
. Furthermore,
determines the poles and
is assumed to be non-characteristic (i.e.,
). The negative integer
determines the leading order term,
, 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
, substitute
into (
1). Balance the minimal power terms in
g, namely,
and
, to get
. Next, require that the leading terms (in
vanish, yielding
Step 2: Compute the resonances. In this step, one determines which functions
in (
5) will remain arbitrary. That happens at specific values of
k, called
resonances, denoted by
r. To find the resonances, substitute
, with
and
given in (
6) or (7) into (
1), and equate the coefficients of the dominant terms (in
) that are
linear in
to get the characteristic equation
Since
the resonances are
, and
Resonance
corresponds to the arbitrariness of
g.
Step 3: Check the compatibility conditions. To do so, substitute
into (
1) and verify that
and
can be
unambiguously determined. Verify also that
and
are indeed arbitrary functions, meaning that no
compatibility conditions arise. For (
1), one readily determines the real functions,
when
(and complex expressions when
). At resonance
, the compatibility condition (arising as the coefficient of
),
is
identically satisfied upon substitution of (
6) and (
10). Likewise, a long but straightforward computation shows that the compatibility condition at
(appearing at order
),
is
identically satisfied upon substitution of (
6), (
10), and (
11).
Thus, at least in the neighborhood of
, 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
. 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
,
and
in (
5).
Moreover, notice that truncating the Laurent series at the constant level term yields an auto-Bäcklund transformation,
because for an arbitrary
, both
and
must then be solutions of (
1) with
. Setting
in (
14) motivates the Hirota transformation discussed in
Section 9, where it is shown that only for
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
, the Gardner equation reduces to the KdV Equation (
2), which is scaling homogeneous under the scaling (dilation) symmetry
where
is an arbitrary scaling parameter. Indeed, replacing
in the KdV equation in terms of
yields
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),
, and
. Equivalently,
and
The
rank of a monomial is defined as the total weight of the monomial. For example,
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
Solving (
17) gives
and
Since
is arbitrary, without loss of generality one can set
, resulting in
and
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
, the Gardner equation becomes the mKdV equation. Requiring uniformity in rank readily yields
and
Hence, the mKdV equation is scaling homogeneous under the transformation
Because
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
yielding
which expresses that (
1) has the scaling symmetry
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
where the dot means that the equality should only hold on solutions
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 etc. The density and flux could also explicitly depend on x and t, but we will not cover such exceptional cases.
In (
22),
is the total derivative with respect to
defined by
where
is the Fréchet derivative of
in the direction of
, and
M is the highest order of
in
x. In practice, one simply applies the chain rule for differentiation with respect to
t, treating
, etc., as independent functions.
Likewise,
is the total derivative with respect to
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
Indeed, (
25) is (
1) written as a conservation law. Next, (26) straightforwardly follows after multiplication of (
1) by
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,
has rank 1,
has rank 2, and
is of rank 4. The corresponding fluxes are also uniform, with ranks
and
As a matter of fact, the entire conservation laws themselves are uniform in rank, with ranks
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., , , , , . For the construction of candidate densities, the constant terms 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, Obviously, if we carry out partial integration, the highest derivative term only differs from by the x-derivative of and therefore can be ignored. Likewise, terms like , , , and can be neglected because they are x-derivatives of single-term monomials (, , etc.). There is no need to put terms like , 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
and
). Finally, linearly combine the remaining terms with constant coefficients, yielding
which is of first order in
x (
).
Step 2: Compute the undetermined coefficients. Using (
23), first compute
where
is the identity operator. Replace
and
from (
1) to get
Next, find the constants
and
so that
matches
for some flux
J (to be computed in Step 3 below). Mathematically, this means that
E must be
exact. The Euler operator (variational derivative),
allows one to test exactness [
21,
22,
23,
42], where
K is the order of the expression the Euler operator is applied to. Therefore,
for
E in (
30). Consequently, (
31) will terminate after five terms.
E will be exact if
on the jet space (treating
, 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
which must vanish identically, yielding the linear system
and
with solution
. Substitute these constants into (
28) to obtain
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
Step 3: Compute the flux. Since
, to get
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
acting on a exact expression
E of order
K, is given [
21,
22,
23,
42,
44] by
with integrand
In (
35),
means that once
is computed one must replace
u by
,
by
,
by
, etc. Use (
34), to get
which already has the terms of
, but still with incorrect coefficients (and the opposite sign). Finally, using (
35),
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
Hence,
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
is called a
generalized symmetry of (
4) if and only if it leaves (
4) invariant under the replacement
within order
, that is,
must hold up to order
on the solutions of (
4). Consequently,
G must satisfy the
linearized equation
where
is the Fréchet derivative of F in the direction of
G:
where
n is the (highest) order of
F. In (
41) and (
43), one must not only replace
u by
but also
by
,
by
, etc. As before,
is the identity operator. The total derivative operators,
and
, 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
With regard to the weights in (
20), the above symmetries are of ranks
, and 6, respectively.
Except for
, each of these symmetries leads to a completely integrable
nonlinear PDE,
, where the weight of
t increases as
j increases. Notice that
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
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., , , , , , , without the constant terms .
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
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
and
can be removed. Finally, linearly combine the remaining terms with constant coefficients, yielding
which is of fifth order (
).
Step 2: Find the undetermined coefficients. Compute
and use (
4) to remove
, etc. This produces an expression with 176 terms (not listed). Next, compute (44) using
and
G from (
48), yielding 193 terms (not listed). Then,
, which has 138 terms, must vanish identically on the jet space, yielding a linear system of 20 equations for the nine nonzero coefficients
,
through
, and
:
For brevity, we have shown only a couple of the shortest equations (coming from
and
, respectively) and two of the longest equations (coming from
and
, respectively). The 18 coefficients
through
,
through
, and
through
are all zero. Solving the system yields
Set
and substitute (
51) into (
48) to obtain
which matches
in (
47).
The code InvariantsSymmetries.m can also be used to verify if an evolution equation, e.g., , with 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 , the software returns in (46), confirming that corresponds to a generalized symmetry of .
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
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]):
where
denotes the anti-derivative (or integral) operator.
connects the symmetries sequentially (without gaps in this case; for examples with gaps, we refer to [
55]):
For example,
which is
, and
which is
.
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
is a hereditary operator [
56].
If
is a recursion operator for (
4), then the Lie derivative [
19,
43,
57,
58] of
is zero, yielding
where
and ∘ denote the commutator and composition of operators, respectively.
is the Fréchet derivative operator, defined as
where
n is the order of
F in (
4).
is the Fréchet derivative of
in the direction of
F. For recursion operators of the form
where
T is the total number of terms in
, and
and
are the orders of
U and
V, respectively, one has
With regard to (
53),
is a linear integro-differential operator [
59,
60] which naturally splits into two pieces [
55,
58,
61],
where
is a local differential operator and
is a non-local integral operator.
is a linear combination of monomials involving
,
u, and parameters with weight (if applicable) so that each monomial has the correct rank. Note also that
will always be “propagated” to the right in
. For example,
Based on the theory of recursion operators,
is a linear combination with constant coefficients of terms of the form
where
is a symmetry and
is a cosymmetry (Euler operator applied to a density
, selected such that
has the correct rank [
58,
62]. To standardize the form of
, one propagates
to the left. For example,
. 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
and
(as well as
and
) differ by 2,
must have rank 2. Based on (
20),
, 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
, and
of rank 2 or less, i.e.,
,
,
, where the trivial terms
and
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
. After stripping off numerical factors and removing duplicates, linearly combine the resulting monomials with constant coefficients to obtain a candidate local operator
Using symmetry
and densities
and
(all of low rank), compute
and
. Then, with terms of type (
63), make the candidate non-local operator
so that each term has rank 2. Add both pieces to get the candidate recursion operator
Step 2: Compute the undetermined coefficients. Separately compute all the pieces in (
57), beginning with
which can also be computed using (
60). With
F given in (
49), insert
into (
67), expand, and simplify. This yields an expression of 27 terms (not listed). Some of the terms have
and
; others involve
,
,
,
, and
. Next, using (
58), compute
With (
66) and (
69), then compute
using, for example,
to consistently move the operator
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,
converts the integral
into
by moving the
operator under the
operator from the right to the left. After integration by parts,
has 58 terms (not listed). Next, compute
To move the operator
to the right, use, for example,
and similar formulas (expressed in a more general form in [
55]). The expanded expression for
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
,
, ⋯, should be treated as independent. This yields
and a linear system,
involving the remaining nonzero constants. Solve the system to get
Finally, set
and substitute (
76) into (
66), yielding
which matches
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
and an auxiliary function
,
where
and
are expressed in powers of
with coefficients depending on
, etc., and
is an eigenfunction of
corresponding to eigenvalue
. To guarantee complete integrability of (
4), at least one non-trivial Lax pair
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
where
is any real or complex number (not necessarily small). Substituting
and
into (
78) yields
The first equation is a Schrödinger-type equation for the arbitrary eigenfunction
with eigenvalue
and potential
. 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
where (
81) and (
82) were used repeatedly to eliminate
,
,
,
,
, and
, 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
and
as follows
where
With (
78), one has
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
and
are considered
trivial). Rearranging the terms yields
or, expressed in operator form by suppressing the
,
where
means that equality holds only on solutions of the original PDE (
4). As before,
is the commutator of the operators, and
is the zero operator. Equation (
88) is called the Lax equation. We now show how the Lax pair
can be computed using the method discussed in [
65].
Step 1: Construct a candidate for
and
. Since the KdV and mKdV equations are special cases of (
1), it makes sense to search for
of rank 2 and
of rank 3. To construct
, list all monomials in
, and
of rank 2 or less, i.e.,
,
,
,
,
,
,
, where the trivial terms
and
have been removed. As explained in
Section 7, apply
to the elements in the second sublist and, as in (
62), propagate
to the right, yielding
. After removing duplicates, linearly combine the monomials from both sublists with constant coefficients to get a candidate for
:
where the coefficient of
has been set to one (for normalization).
To make a candidate for
, list all monomials in
, and
of rank 3 or less, i.e.,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
, where the trivial terms
,
, and
have been removed. Apply
to the elements in the second sublist, yielding
,
,
,
,
,
,
,
. Apply
to the elements in the third sublist and use (
62) to obtain
,
,
,
. After stripping off numerical factors and removing duplicates, linearly combine the resulting monomials with constant coefficients to obtain a candidate for
:
Step 2: Compute the undetermined coefficients. First, compute
which, after substituting
F from (
49) to replace
and
, yields 14 terms. Next, compute
which upon expansion has 125 terms, and
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
. Set the coefficients of
, and
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:
For brevity, we have shown only a couple of the shortest equations (coming from
and
respectively), and two of the longest equations (coming from
and
, 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
where
and
are arbitrary constants. To be able to obtain the Lax pair for the KdV equation where
, one should require that
and
, otherwise
,
,
, and
would become infinite. Notice that both requirements allow one to clear
from all denominators. Furthermore, the coefficients then simplify into
Finally, substitute the coefficients into (
89) and (
90) to get
where the constant
is arbitrary. Hence, this is a one-parameter family of Lax pairs. Set
to get (
79) and (80). With (
97) and (98)
which, after setting
, is equivalent to (
83).
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],
which connects solutions
of the Gardner equation,
with
, to solutions
of the KdV equation,
with arbitrary coefficient
. The use of
will avoid confusion with
in (
121) and, more importantly, allow us to set
to get the standard Miura transformation (see
Section 12).
Substituting (
120) into (
122), it is straightforward to verify that
where
means that the left hand side is evaluated on solutions of (
121). As before,
is the identity operator and
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
,
, and
. Therefore, (
120) is uniform in rank of rank two.
Step 1: Construct a candidate for the Gardner transformation. Make the list
,
of monomials in
v and
of rank 2 or less. By differentiating its elements, replace the second sublist by
. Linearly combine the resulting elements with undetermined coefficients,
to generate a candidate for the Gardner transformation.
Step 2: Compute the undetermined coefficients. Substitute (
124) into (
122) and, using (
121), replace
and
. The resulting expression must vanish on the jet space, leading to
and half a dozen more complicated equations. For
, (
124) would become an algebraic transformation. Hence,
and, after simplification of the more complicated equations, one is left with the following
nonlinear system
Substitute the solution,
,
,
, and
, 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
,
, and
can occur in that operator. Apply the candidate for the operator,
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
Solve the system to obtain
,
, and
. 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
, (
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
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,
follows from (
79) and (80) by setting
.
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],
where
is the modulus of the Jacobi elliptic cosine (
) function. Both solutions are depicted in
Figure 1. When
m approaches 1, the peaks of the
-squared solution become a bit taller, and the valleys become lower and flatter before they spread out horizontally to become the
-squared solution. Both solutions (
129) and (130) satisfy
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
in (
25)–(27) and (
45)–(47), respectively. The recursion operator [
54] connecting those symmetries follows from (
53) with
.
Likewise, a one parameter Lax pair for the mKdV equation follows from (
79) and (80) by setting
:
which for
simplify into
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],
which connects solutions of the mKdV Equation (
3) (with dependent variable
) to solutions
of the KdV Equation (
122) readily follows from (
120) by setting
. Likewise, (
123) reduces to
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.
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
into (
116), yielding
and (117), which is identically satisfied. From (
118), one then obtains
where
, which matches (
138). The two-soliton solution [
5] follows from
with
and
. Then, from (
118)
Details of the derivation are given in [
5] and [
6], where formulas for the three- and
N-soliton solutions can also be found.
A new exact two-soliton solution of (
114) with
is presented in [
12]:
with
and
with
and
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
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
will be transformed by (
144) into one that goes to
as
. 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,
, and Jacobi elliptic function methods, at a click of a button, one obtains simple exact solutions, which we cover first.
Figure 2.
Graphs of (
145) (
left) and (
146) (
right) both with the minus signs in front of the square roots, and both for
,
,
, and
. The curves on the right correspond to
(dashed line) and
(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
,
,
, and
. The curves on the right correspond to
(dashed line) and
(solid line).
Figure 3.
Graphs of (
149) (
left) and (
150) (
right), both for
,
,
, and
. The curves on the right correspond to
(dashed line) and
(solid line).
Figure 3.
Graphs of (
149) (
left) and (
150) (
right), both for
,
,
, and
. The curves on the right correspond to
(dashed line) and
(solid line).
The graphs in
Figure 4 are for
,
, and
(i.e.,
is equivalent to
), for which (
148) simplifies into
As the value of
k increases, the solitary wave get taller and narrower. At
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
, the wave degenerates into a horizontal line corresponding to
For values of
k near 1, solution (
151) can be very well approximated by a kink–antikink pair [
79],
where
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
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
, tanh,
, etc., have no weights. With regard to the weights (
20),
, as it should, because
Furthermore,
because
and
. All terms in any
must have weight zero, in particular,
, and
, where
m is the modulus of any of the Jacobi elliptic functions. From
, it follows that
and
, where
is the angular frequency in
. Hence, if
and
V are polynomials in
k, then
can only have terms proportional to
and
, and
V can only have terms in
, and
. The proportionality factors could have any powers of
since
.
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
-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
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
,
, and
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
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.