arXiv:1109.3745 [nlin.PS]

A KdV-like advection-dispersion equation with some remarkable properties
Abhijit Sen1 , Dilip P. Ahalpara1 , Anantanarayanan Thyagaraja2 and
Govind S. Krishnaswami3

arXiv:1109.3745v3 [nlin.PS] 21 Sep 2014

1

2

Institute for Plasma Research, Bhat, Gandhinagar 382 428, India,
H.H. Wills Physics Laboratory, University of Bristol, Tyndall Avenue, BS8 1TL, UK,
3
Chennai Mathematical Institute, SIPCOT IT Park, Siruseri 603103, India.
Email: senabhijit@gmail.com, dpa1951@gmail.com,
a.thyagaraja@bristol.ac.uk & govind@cmi.ac.in
February 25, 2012

Abstract
We discuss a new non-linear PDE, ut + (2u xx /u)u x = u xxx , invariant under scaling of dependent variable and referred to here as SIdV. It is one of the simplest such translation and space-time
reflection-symmetric first order advection-dispersion equations. This PDE (with dispersion coefficient unity) was discovered in a genetic programming search for equations sharing the KdV
solitary wave solution. It provides a bridge between non-linear advection, diffusion and dispersion. Special cases include the mKdV and linear dispersive equations. We identify two conservation laws, though initial investigations indicate that SIdV does not follow from a polynomial
Lagrangian of the KdV sort. Nevertheless, it possesses solitary and periodic travelling waves.
Moreover, numerical simulations reveal recurrence properties usually associated with integrable
systems. KdV and SIdV are the simplest in an infinite dimensional family of equations sharing
the KdV solitary wave. SIdV and its generalizations may serve as a testing ground for numerical
and analytical techniques and be a rich source for further explorations.

Keywords: genetic programming, advection dispersion equation, travelling waves, recurrence

Communications in Nonlinear Science and Numerical Simulation 17 (2012), pp. 4115-4124

Contents
1

Introduction

2

2

General properties of the SIdV equation

4

2.1

Symmetries and Conservation laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2

Preservation of positivity of u(x) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.3

Behaviour of SIdV at some special values of dispersion coefficient . . . . . . . . . . . . . . . . . . . . . .

5

2.3.1

Reduction to a linear dispersive wave equation when  = −2/3 . . . . . . . . . . . . . . . . . . . .

5

2.3.2

The dispersionless limit  = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.3.3

SIdV to mKdV when  = 2/3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

3

Some similarity and travelling wave solutions of SIdV

6

3.1

Similarity solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

3.2

Travelling waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

4

Search for a variational principle for SIdV

9

5

Numerical evolution of SIdV solitary waves: recurrent behavior

10

6

Advection-dispersion equations sharing KdV solitary waves

11

7

Discussion and conclusion

13

A Computer-aided discovery of SIdV equation

1

14

Introduction

Several equations of physics have been discovered by searching for one that admits a particular type
of solution that was known to exist on physical grounds. Perhaps the most celebrated example is
the discovery of Schrödinger’s equation of quantum mechanics via de Broglie’s hypothesis that free
particles are described by plane matter waves. The KdV equation of fluid flow in a canal [1],
ut + 6uu x + u xxx = 0

(1)

was motivated (and in part discovered) by the search for an equation possessing Russell’s [2] ‘wave
of translation’ as a solution. This is the solitary wave
!
√
c
c
2
u(x, t) = sech
(x − ct − x0 )
for c > 0, x0 ∈ R.
(2)
2
2
The KdV equation is the simplest conservative 1-dimensional wave equation with weak advective
non-linearity and dispersion. Thus it is widely applicable and has been used to model acoustic solitons in plasmas [3, 4], internal gravity waves [5] in the oceans and even blood pressure pulses [6].
Its significance was greatly amplified by Zabusky & Kruskal’s discovery [3] that KdV displays the
Fermi-Pasta-Ulam recurrent behavior and lack of thermalization [7] in a spatially periodic domain
(see also [8, 9]). This was subsequently attributed to the even more remarkable asymptotic superposition principle for scattering of KdV solitary waves and existence of infinitely many local conserved
2

quantities in involution [10]. Moreover, KdV is the prototype for an integrable non-linear PDE, its
initial value problem can be solved by the inverse scattering transform [11].
We were therefore excited to find that the KdV equation is not the only one with the sech 2 solitary
wave solution. This serendipitous discovery happened in the course of an investigation undertaken
by two of the present authors [12] to improve the efficiency and accuracy of a genetic programming
(GP) based method (see A) to deduce model equations from a known analytic solution [13]. As a
benchmark exercise to test out the method for application to nonlinear PDEs, the travelling wave (2)
was given to the program, expecting it to find the KdV equation. But, surprisingly, before finding the
KdV equation, it found a different equation
!
2u xx
ut +
u x = u xxx
(3)
u
which had the same solitary wave solution. Subsequently, we found that (3) is the simplest in a vast
family of equations sharing the KdV solitary wave. We think of (3) as a non-linear wave equation
for the dispersive advection of the real wave amplitude u. Unlike the KdV equation, where the
advecting velocity V = 6u is linear, here it is a quotient V = (2u xx /u). As a consequence (3) is
‘scale-invariant’ under dilation u → λu of dependent variable. Scale-invariant advective velocities
have appeared before e.g., the E×B
velocity that is invariant under a rescaling of fields, charges and
B2
currents in a plasma. We refer to (3) as the ( = 1 case of the) SIdV equation1 . Like KdV, faster SIdV
solitary waves are narrower, but due to scale invariance, height and speed are generally unrelated, as
in classical linear wave equations, but unlike KdV and the non-linear Schrödinger equation (NLSE).
In (3) the dimension- L3 /T coefficients of the dispersive and advective terms are both equal to
one. This is of course very special. More generally we consider the SIdV equation
!
2au xx
ut +
u x =  a u xxx .
(4)
u
√
x → 3 ax eliminates the L3 /T -dimensional constant a, leaving one dimensionless parameter 
measuring the strength of dispersion relative to advection. Unlike KdV, where they scale differently
in x, here both scale as L−3 . Though we have not found an experimental system modelled by (4),
it is just as simple and universal among scale-invariant advection-dispersion equations as KdV is
among all such equations. For generic  , we have identified 2 conserved densities. This is similar to
inviscid Eulerian hydrodynamics, but unlike KdV and NLSE which possess an infinite number. While
SIdV shares solitary waves with KdV at  = 1, non-linear diffusion and solvability emerge elsewhere.
Despite being non-linear, scale-invariance ensures SIdV has exact plane wave solutions. Furthermore,
it possesses bounded spatially periodic travelling waves and similarity solutions. However, SIdV
cannot arise from a polynomial Lagrangian in the sort of variable that works for KdV. We evade this
obstruction at some special values of  . Remarkably, numerical evolution of (3) shows Fermi-PastaUlam-like Birkhoff recurrence [9] despite no sign of soliton scattering!
1

The acronym SIdV highlighting scale-invariance is related to KdV as Sine-Gordon is related to Klein-Gordon.

3

2

General properties of the SIdV equation

2.1

Symmetries and Conservation laws

SIdV (4) is a non-linear advection-dispersion equation for the real wave amplitude u(x, t) that is being
advected by the flow V = 2u xx /u. It is scale-invariant2 under u → λu. In fact, (4) is the simplest
non-linear translation and scale-invariant advection-dispersion equation that is first order in time. The
lowest order dispersive term is u xxx , so any such equation can be written as ut + Vu x = au xxx for
some scale-invariant advective velocity V . Without requiring scale-invariance, the simplest choice
V ∝ u leads to KdV. KdV is symmetric under space-time (PT) reflection (x, t) → (−x, −t). Now
if we also require scale-invariance, the simplest advecting velocities that preserve PT symmetry are
V ∝ uuxx and uuxx . The former leads to (4) in units where a = 1. The latter choice too has some notable
properties (see §7).
Non-zero constants are the simplest solutions of (4). Linearization about a constant yields a plane
wave with cubic dispersion ω = k3 and the characteristic property that the ratio of phase to group
velocity is 1/3, as in the linear KdV equation.
Remarkably,

 despite being non-linear, (4) also admits
exact plane wave solutions u = A sin kx − ( − 2)k3 t + ϕ .
KdV is invariant under Galilean boosts x → x + ct if u transforms as u → u(x + ct, t) − (c/6).
However, SIdV is not invariant under δu = b+ctu x for any constants b and c , 0. However, Galilean
boosts could be implemented in a more intricate manner that we have not identified.


Multiplying by u, and using uu3x = uu xx − 21 u2x , (4) is written in conservation form
x


1 h 2 i 

u + 1 + u2x − uu xx = 0.
t
2
2
x

(5)

R
It follows that I = u2 dx is conserved. This suggests u2 is the concentration of some substance
whose total amount is conserved. Similarly, multiplying (4) by u−2/ we get




u1−2/ + (2 − ) u−2/ u xx = 0.
(6)
t

x

R
So J = u1−2/ dx is also conserved. For e.g., when  = 1, J = u1 dx, so on a bounded domain, J
is finite for any strictly positive/negative initial condition. These integrals of motion, travelling waves
(§3.2) and numerical evolution (§5) indicate that SIdV is generically non-dissipative.
R

2.2

Preservation of positivity of u(x)

At first sight, it appears that u = 0 is a singular point of (4). But u can vanish at points where u x
or u xx also vanish, provided u xx u x /u is finite, e.g. the above plane wave vanishes at isolated points.
Among travelling waves, this is generic, near one where u and u xx have common zeros, there is
another solution with the same property. However, in the numerical and analytical examples studied,
R
if u(x, 0) > 0, it remains positive for t > 0. Let us use the second conserved quantity J = 1u dx
to sketch why this is the case for  = 1 on a bounded domain3 . Suppose u(x, 0) > 0 and J is finite
at t = 0. At t1 > 0, let u(x) develop its first zero, this cannot be a first-order zero as u was strictly
2
3

SIdV (4) is not invariant under dilations of x and has a dimensional scale a which we set to 1 .
A similar argument can be given for any  for which a conserved density diverges at u = 0 sufficiently fast.

4

positive. Then J(t1 ) = ∞, contradicting the constancy of J ! So strictly positive initial data cannot
develop a zero and therefore must remain positive4 .

2.3
2.3.1

Behaviour of SIdV at some special values of dispersion coefficient
Reduction to a linear dispersive wave equation when  = −2/3

For  = −2/3, SIdV may be reduced to a linear dispersive equation for u2 . If we write uu xx =
(uu x ) x − u2x = 12 (u2 ) x − u2x in the conservation law form (5) of SIdV, we get
"
!
#
1 2
3 2   2 
∂t u + ∂ x 1 +
u − u
= 0.
(7)
xx
2
2 x 2
If  = −2/3, we see that ρ = u2 satisfies a linear KdV equation
2
ρt + ρ xxx = 0.
3

(8)

So for  = − 32 , ρ = u2 is a sort of Cole-Hopf transformation that linearizes the equation. The general
solution is a linear combination of plane waves
Z
Z
2
i(kx+ 32 k3 t) dk
ρ(x, t) = u (x, t) =
ρ̃(k) e
where ρ̃(k) =
u2 (x, 0) e−ikx dx.
(9)
2π
R
R
This solution could serve as the 0th order of a perturbative solution of SIdV for nearby  .
2.3.2

The dispersionless limit  = 0

The special case  = 0 of (4) gives a non-dispersive non-linear advection equation
2u xx
u x = 0.
(10)
u
 
It may be written as a conservation law (u2 )t + 2u2x = 0 for the ‘charge’ density u2 with flux u2x .
x
Being a second order parabolic PDE, (10) may also be viewed as an unusual non-linear diffusion
equation ut = αu xx for the ‘temperature’ u. The effective thermal diffusivity α = −2u x /u could be of
either sign. Thus SIdV is a remarkable bridge connecting dispersion, non-linear advection and nonlinear diffusion. We expect (10) to have instabilities if α becomes negative since time-reversed heat
equations are ill-posed. So we may think of the dispersive term in SIdV as a regularization of (10),
just as KdV is a regularization of the kinematic wave equation (KWE) ut + 6uu x = 0. Remarkably,
even without a dispersive regularization, (10) has smooth solutions. Indeed, unlike the KWE, (10)
admits waves that preserve their shape. The general travelling wave of (10) is u = A cos(kx+2k3 t+φ).
These plane waves may however be unstable, as the diffusivity oscillates in sign!
ut +

4
We assume that u(x) cannot develop an ‘integrable’ zero for which J is finite. An ‘integrable’ zero where u ∼ |x|a
for a < 1 would mean u forms a cusp and ceases to be thrice continuously differentiable. Assuming the solution remains
sufficiently smooth, such possibilities are eliminated. Our numerical simulations did not indicate cusp formation, though it
is an open question whether SIdV preserves regularity of initial data.

5

2.3.3

SIdV to mKdV when  = 2/3

We did not find any value of  at which a transformation reduces the SIdV to the KdV equation. But
there are special values of  at which it comes close. First, the rational non-linearity of SIdV can be
written as a polynomial while retaining the advection-dispersion structure of the equation. To do so,
we use the invariance of SIdV under rescaling to choose u dimensionless and put u = ew . Then5
h
i
wt + (2 − )w2x + (2 − 3)w xx w x = w xxx .
(11)
This polynomial form of SIdV was convenient for numerical evolution and also in our search for a
Lagrangian. It also indicates that  = 2, 32 are somewhat special. At these values, we get KdV-like
dispersive wave equations with advecting velocities ∝ w2x and w xx . These are among the simplest PT
symmetric advecting velocities beyond KdV. Moreover, at  = 2/3, the sign of the ‘local diffusivity’
is reversed. We see qualitative effects of this reversal in the stability of our numerical simulations as
 is decreased below 2/3. What is more, at  = 2/3 SIdV is reducible to the modified KdV (mKdV)
equation. Differentiating in x at  = 2/3, putting v = w x and letting x → −x and t → 23 t , we get the
defocusing mKdV equation, which is integrable and related to KdV via the Miura transform [14]
vt − 6v2 v x + v3x = 0.

3

(12)

Some similarity and travelling wave solutions of SIdV

3.1

Similarity solutions

SIdV is invariant under6 x → λx and t → λ3 t . So we seek solutions u(x, t) = f (z) in the similarity
variable z = x3 /54t . We get a third order non-linear ODE for f (z)


z2 000
2z
f f − z2 f 00 f 0 + z f 00 f − f 02 +
+ z f 0 f = 0.
2
3
9

(13)

By the substitution g = f 0 / f we reduce this to a cubic 2nd order ODE with variable coefficients
!
!




z2 00
3
2
2 0
0
2 3
g +
− 1 z gg + zg +
− 1 z g +  − zg2 +
+ z g = 0.
(14)
2
2
2
3
9
We haven’t solved the similarity ODE in general, but in the dispersionless limit  = 0, it becomes a
linear ODE with a regular singularity at z = 0
z f 00 +

2 0
f − f = 0.
3

(15)

The two linearly independent solutions may be expressed in terms of the confluent hypergeometric
function 0 F1 (a, z) or the modified Bessel function of the 1st kind In (z):
!
!
 √
2
3z
2 1
f1 (z) = 0 F1 , z = 1 +
+ ... = Γ
z 6 I− 1 2 z
and
3
3
2
3
5

Of course, if we restrict to real w , this equation will apply to solutions where u remains everywhere positive.
More generally, SIdV is invariant under x → λx, t → λ3 t, u → λγ u for any γ . So t−γ/3 u(x, t) = f (z) is a scale invariant
combination. Since γ is arbitrary, we restrict here to the simplest case γ = 0 .
6

6

f2 (z) = z

1
3

4
,z
3

0 F1

!
= z

1
3

3z
1+
+ ...
4

!

!
4 1  √
= Γ
z6 I1 2 z .
3
3

(16)

√

Both f1 and f2 are monotonic and grow ∝ e2 z as z → ∞. So they are bounded at late times (z → 0
or t  x3 ) but are unbounded at early times (z → ∞ or x3  t → 0). Interestingly, there is a
unique (up to scale) linear combination f (z) = Γ(4/3) f1 (z) − Γ(2/3) f2 (z) that is bounded for all z ≥ 0
( x, t ≥ 0). It begins at f (0) = Γ(4/3) and monotonically decays to zero as z → ∞.

3.2

Travelling waves

Here we discuss travelling waves (u = f (ξ) with ξ = x − ct ) for SIdV (4) on the unbounded domain
−∞ < ξ < ∞. Travelling waves must satisfy the third order non-linear ODE −c f f 0 +2 f 0 f 00 − f f 000 =
0. We may write this in ‘conservation law’ form
−

c  2 0 2 +   02 0
f +
f
− ( f f 00 )0 = 0
2
2

(17)

and integrate once to get
2 f f 00 − ( + 2) f 02 + c f 2 + 3B = 0.

(18)

The substitutions p = f 0 and F = p2 give us a first order linear ODE for F( f ):
 f F 0 ( f ) − ( + 2)F + c f 2 + 3B = 0

or

F 0 (w) − ( + 2)F + ce2w + 3B = 0,

(19)

where f = ew . This inhomogeneous 1st order ODE is reduced to quadrature using the integrating
factor e−(1+2/)w . For  , 0, in terms of r(w) = e−(1+2/)w F(w), we get
r0 + 3Be−(1+2/)w + ce(1−2/)w = 0.

(20)

For  , 0, ±2, we integrate to get
2
3B
c e2w
F(w) = Ae(1+  )w +
+
.
2−
+2

(21)

In other words, the equation for travelling waves has been reduced to quadrature

Z

dξ = ±

Z

df
p
F( f )

where

 c 2 3B



2f + 2



1
3B


 2 A f 2 − c f 2 log f + 4
F=

c 2
3B

f − 1A


4 f + 2 log

  2

2

A f (1+  ) + c f 2 + 3B
2−
+2

if  = 0.
if  = 2,
if  = −2, and
otherwise.

(22)

This travelling wave integral can be understood by a mechanical analogy [14]. It is the zero ‘energy’
condition E = f 02 + V( f ) = 0 for the ‘coordinate’ f (ξ) at ‘time’ ξ of a non-relativistic particle of
mass 2 moving in the 1-dimensional potential V( f ) = −F( f ). Bounded travelling waves correspond
to bound trajectories of this particle. Plotting V( f ) shows that, for appropriate ranges of A, B and
c, there are bounded spatially periodic/solitary travelling waves with heights
R between successive real
simple/double zeros of V( f ). For generic  , the travelling wave integral F −1/2 d f cannot be evaluated using elementary functions. But for  = ∞, 0, −2/3, the integral is trigonometric/exponential
and for  = 2/3, 1 it is elliptic.
7

To illustrate, we consider the case  =R 1 where SIdV shares solitary wave solutions with KdV.
Here we must evaluate an elliptic integral (A f 3 + c f 2 + B)−1/2 d f , and bounded travelling waves are
(limits of) cnoidal waves. Suppose A > 0 and F( f ) = A f 3 + c f 2 + B = A( f − f1 )( f − f2 )( f − f3 ) has
three simple real zeros 0 ≥ f1 < f2 < f3 ≥ 0. Then we have a periodic cnoidal wave with trough at
f1 and crest at f2 , determined by
Z f
dg
± (ξ − ξ1 ) =
where f (ξ1 ) = f1 .
(23)
p
A(g − f1 )(g − f2 )(g − f3 )
f1
Transforming to g = f1 + ( f2 − f1 ) sin2 θ and defining the shape parameter 0 ≤ m = ff32 −− ff11 ≤ 1, (23)
becomes a standard incomplete elliptic integral of the first kind. It is inverted in terms of a Jacobi
elliptic function cn(u; m) = cos φ, with modulus k (m = k2 )
Z φ
dθ
2u
where u =
.
(24)
ξ = ξ1 ± p
p
A( f3 − f1 )
0
1 − m sin2 θ
At the upper limit g = f = f2 − ( f2 − f1 ) cos2 φ, so the cnoidal wave for A > 0 is
!
p
2 1
f = f2 − ( f2 − f1 ) cn
A( f3 − f1 )(ξ − ξ1 ) ; m .
2

(25)

Its shape depends on m while its wavelength and speed are
λ= p
Here K(m) =

R π/2
0

√

dθ

1−m sin2 θ

4K(m)
A( f3 − f1 )

and

c = −A( f1 + f2 + f3 ).

(26)

is the complete elliptic integral of the 1st kind.
00

The advecting velocity field V = 2 ff for cnoidal waves is finite since f and f 00 have common
zeros. Moreover, by modifying the parameters A, fi , we get nearby waves with the same feature.
Unlike for KdV, where the shape, speed and wavelength of cnoidal waves are non-trivially modified
upon a rescaling of amplitude, here f 7→ λ f produces a new cnoidal wave with the same m, c, λ and
phase ξ1 , since the constants transform as (A, B, fi ) 7→ (A/λ, λ2 B, λ fi ).
If A < 0, the cnoidal wave extends between f2 and f3 and is given by
!
p
2 1
f (ξ) = f2 + ( f3 − f2 ) cn
A( f1 − f3 ) (ξ − ξ3 ) ; m̃
2
p
where m̃ = ( f3 − f2 )/( f3 − f1 ) and λ = 4K(m̃) A( f1 − f3 ).

(27)

Solitary waves are cnoidal waves of infinite wavelength. They occur when a pair of simple zeros
of F coalesce to form a double zero. For example, if f3 → f2 in (25) holding f1,2 and A > 0 fixed,
then m → 1− , K(m) → ∞ and we get a left-moving solitary wave of depression
!
p
2 1
f (ξ) → f2 − ( f2 − f1 ) sech
A( f2 − f1 )(ξ − ξ1 )
for A > 0.
(28)
2
If the zeros f1 → f2 coalesce in (27) we get a left-moving solitary wave of elevation
!
1p
f (ξ) → f2 + ( f3 − f2 ) sech 2
A( f2 − f3 )(ξ − ξ3 )
for A < 0.
2
Finally, if A → 0, cnoidal waves reduce to sinusoidal waves with cubic dispersion,
√

f (ξ) = N sin −c (ξ − ξ0 ) , c < 0, N arbitrary.
8

(29)

(30)

4

Search for a variational principle for SIdV

To relate symmetries to integrals of motion, it is interesting to find a Lagrangian or Hamiltonian
formulation for SIdV. The existence of a Lagrangian for a given equation depends on the field
variables used and the sort of Lagrangian allowed. For instance, the dispersive term u xxx in KdV
ut + 6uu x + u xxx = 0 isn’t the variation of any polynomial in u and its derivatives. This is because
every quadratic differential polynomial in u involving three x-derivatives is a total derivative. Yet, as
is well-known, if we put u = χ x , KdV follows from the Lagrangian density 21 χt χ x + χ3x − 21 χ2xx .
In looking for a Lagrangian for SIdV, we choose to work with w = log u, which satisfies the
KdV-like equation (11) with polynomial non-linearity. By analogy with KdV, we put w = φ x and
seek a polynomial action in φ and its derivatives, whose Euler-Lagrange (EL) equations are
h
i
φ xt + (2 − )φ2xx + (2 − 3)φ xxx φ xx = φ xxxx .
(31)
φ is a natural variable since the linear part of the equation φ xt = φ4x admits the polynomial Lagrangian L0 = 21 φt φ x + 2 φ2xx . We wish to add potentials V, W to L0 to reproduce the quadratic
φ xx φ3x and cubic terms φ3xx in (31). To do so, we note a couple of general features. If V is a monomial of degree n in φ and its x-derivatives, then (1) the resulting terms in the EL equation form a
differential polynomial of degree n − 1, and (2) the total number of x-derivatives in each term of
the differential polynomial are the same as in V . Therefore, to produce φ3xx in the EL equation, W
must be a quartic differential polynomial with 6 x-derivatives, and to give φ xx φ3x , V must be a cubic
differential polynomial with 5 derivatives. In other words,
X
V = a1 φφφ5x + a2 φφ x φ4x + a3 φφ2x φ3x + a4 φ x φ x φ3x + a5 φ x φ2x φ2x =
ai Vi .
(32)
i

Can ai be chosen so that the variation of V dx gives (2 − 3)φ xx φ3x ? Unfortunately not, as7
Z
δ
V(φ(y)) dy = − (10a1 − 5a2 + a3 + 4a4 − 2a5 ) (2φ xx φ xxx + φ x φ4x ) .
δφ(x)
R

(33)

For no choice of ai can we produce just a quadratic monomial ∝ φ xx φ xxx in the EL equation. Similarly, we showed that there is no quartic differential polynomial W that gives φ3xx upon variation. We
conclude that there is no polynomial Lagrangian in φ and its derivatives leading to SIdV.
However, the polynomiality assumption is quite strong. There may be a non-polynomial Lagrangian in φ or one in a variable non-locally related to φ. An interesting example of such a possibility occurs when  = − 32 and SIdV becomes the linear dispersive equation ρt + (2/3)ρ xxx = 0
uponR substituting ρ = u2 . This equation follows from L = 12 ψt ψ x − 31 ψ2xx , where ψ x = ρ. However,
x 2φy
ψ=
e dy is non-locally related to φ, so the Lagrangian is non-local in φ.
Another way around this negative result is that there may be a Hamiltonian that is a differential
polynomial in φ, but with non-canonical Poisson brackets. Such a possibility is realized if  = 2/3,
when SIdV can be transformed into the mKdV equation vt −6v2 v x +v3x = 0 by the substitution v = uux
and a rescaling (§2.3.3). mKdV admits a Hamiltonian formulation vt = {H, v} with
!
Z
1
v4
H=
v2x +
dx and {v(x), v(y)} = −∂ x δ(x − y).
(34)
2
6
So for  = 2/3, SIdV admits a polynomial Hamiltonian in the variable v = (log u) x = w x = φ xx . It
would be interesting to find a Hamiltonian/Lagrangian formulation of SIdV for other values of  .
7

In particular, this means all the Vi differ from one another by total x -derivatives.

9

5

Numerical evolution of SIdV solitary waves: recurrent behavior

We numerically solved8 the SIdV initial value problem for  = 1 on the interval [−π, π] with periodic
boundary conditions. Numerical evolution of one sech 2 wave produced a right-moving travelling
wave, as expected from the exact solution on (−∞, ∞). We also considered two solitary waves9
!
!
√
√
X
c1
c2
2
2
(x − x1 + 2π j) + A2 sech
(x − x2 + 2π j) .
u(x, 0) =
A1 sech
(35)
2
2
j
When the two solitary waves were initially separated by some distance but had the same heights and
speed, they were observed to travel without much interaction, just like individual solitary waves.
Next, we gave the two waves the same initial amplitudes A1,2 = 1 but different speeds (or widths)
(c1 , c2 ) = (4, 2) and locations (x1 , x2 ) = (−π/2, π/4) at t = 0 (fig. 1(a)). So it would take each
wave by itself a time of π2 and π to traverse the 2π-interval. Solitary waves that decay at ∞ must be
right-moving, so we couldn’t give them opposing velocities. But c1 > c2 , so wave-1 caught up with
wave-2 due to periodic boundary conditions and collided with it from the rear. Then they separated
into a small leading wave and a larger trailing wave. The original solitary waves did not retain their
shapes. The smaller wave that emerged from the collision moved faster and caught up with the bigger
one by going round the circle. During the next collision, the smaller wave rear-ended the larger one.
The large one in front morphed into a fast small wave, leaving behind a slow large wave. This is
illustrated in the last three plots of fig. 1(a). Qualitatively, this pattern seemed to repeat as the IVP
was solved up to t = 40, allowing more than two dozen collisions to be observed10 .
Though not KdV solitons, the solitary waves that were involved in these collisions displayed a
certain coherence, they did not dissipate nor degenerate into ripples. Despite not being periodic, the
evolution seemed to approximately revisit earlier configurations. Interestingly, there was no equipartitioning of wave intensity. We illustrate this in fig. 2(a) by plotting the absolute squares of the first
few Fourier coefficients of u(x) as a function of time11 . There is some exchange of intensity among
the first 3 or so Fourier modes c0,1,2 , with an approximate periodicity of T ≈ 3. But there is no
appreciable leakage to higher Fourier modes |c3,4,5 |2 , which are uniformly three to six orders of magnitude smaller than |c0 |2 . There is some growth in |c4,5 |2 . But this can’t be distinguished
from an
R
2
accumulation of numerical errors, which also caused the integral of motion I = u dx to increase
by 0.3% over a time 0 ≤ t ≤ 40 (fig. 1(b)).
What is more, though the Raleigh quotient (‘gradient energy’ or mean square mode number)
R
P 2
|u x (x, t)|2 dx
n |cn (t)|2
Q(t) = R
= Pn
(36)
2
|u(x, t)|2 dx
n |cn (t)|
is not conserved, it seems
to oscillate between bounded limits (fig 2(b)). The rms mode number
√
hovers around ν = Q ≈ 12 . A calibration of Q using the linear dispersive equation ut = u3x
8
The evolution was done using NDSolve on Mathematica. Stability of the numerical evolution was slightly enhanced
by working with w = log u which satisfies wt = w3x + w x w xx − w3x . Positivity of u(x, 0) was preserved at all times.
9
The sum over j ∈ Z ensures that the initial condition satisfies periodic boundary conditions u(−π) = u(π) . In practice,
the sum was restricted to | j| ≤ 3 since the remaining terms are exponentially small for −π ≤ x ≤ π .
10
A maximum of 1200 grid points were placed at an average spacing of 0.1% the domain width (2π) . The qualitative
features reported here were unchanged by adding 200 more points. Over the time interval 0 ≤ t ≤ 40 , 1800 nearly equal
steps were taken
R with an average time step of 1.4% of the time it took wave- 1 in isolation to traverse the domain.
1
11
cn (t) = 2π
u(x, t) e−inx dx , the negative coefficients |c−n |2 = |cn |2 contain no new information for real u(x, t) .

10

1
0

0.617517

1.4
1.2
1.0
0.8
0.6
0.4
-3

-2

-1

1

2

-3

3

-2

-1

-3

-2

-1

0.6695
1

2

3

0.6690

1.45

1.4
1.2
1.0
0.8
0.6
0.4

1.4
1.2
1.0
0.8
0.6
0.4
0

1

2

-3

3

-2

-1

-Π

0.6700

0

1.03671

2
à u âx Versus time

0.6705

1.4
1.2
1.0
0.8
0.6
0.4
0

Π

2Π

0.6685

0.6680

0

1

2

3

10

(a) u(x) at t = 0 , before, during and after collision of
right-moving waves on 2π -periodic interval.

20

30

40

R
(b) u2 dx is conserved upto numerical errors, which
cause it to grow by about .3% for 0 ≤ t ≤ 40 .

Figure 1: SIdV evolution of two solitary wave initial state. Times 0 ≤ t ≤ 40 are indicated above the
plots in fig. 1(a).
RMS mode number Ν =

Log10 Ècn 2 versus time for Fourier modes n=0,1,2,3,4,5
10

20

30

40

æ

-1

0.65

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 12 1 1 1 1
2
222222 222 222 222222 222 222222 222 22 2 2 22
-2 2 2 2 2
2
2
2
2
2
2
2
2
2

0.60

3

3

3

3

3

3

3

3

3

3

33

3

3

3

5

33

3

æ
æ
æ
æ
æ

3

æ

æ

æ

æ

3

33
3
-3 3 3
3
3
33
33
3
4 34
3
4
4
44
33
4
4
4
4 4 4
4
4
3
5
5
4
3
4
4 5 5
3443 3
3
5 543
4
-4
44
4
4 44
3
5
5
5
444
54 4
4 4 4
44
35 553 5
4 4 4
5 4
44
4
5
4
4
4
3
5
5
4 5
4
5
5 55 53
5 55 5
-5
55
5
5
5 5
5
5 5
45
555
5
545
5
5
3
3
5
5
-6
3

Q vs time Hfrom truncated Fourier seriesL

0.70

æ

æ

æ

ææ

æ

æ
æ

æ
æ

æ
æ

æ

0.55

æ

æ
æ

æ

æ

0.50

æ

æ

æ
æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

0.45

æ
æ

5
5

-7

2

(a) log10 |cn (t)| for modes 0 ≤ n ≤ 5 (top to bottom).

0

10

20

30

40

th

(b) Rayleigh quotient from 5 order Fourier series

Figure 2: Time dependence of Fourier coefficients of u(x) and gradient energy for a two-solitarywave initial state.
√
indicates that there are about 2N + 1 = 12Q + 1 ≈ 2 active degrees of freedom present, c0 being
the dominant one, with some contribution from c±1 . These numerical simulations indicate that in a
periodic domain, the SIdV equation displays recurrent behaviour (for  = 1) despite possessing only
two (known) constants of motion.

6

Advection-dispersion equations sharing KdV solitary waves

One of the questions that puzzled us after the discovery of (3) by genetic programming, was whether
there are other such non-linear advection-dispersion equations, sharing the sech 2 wave with KdV
[16, 17, 18]. To explore this question let us consider the following generalized form of an advection
dispersion equation,
ut + Vu x = δ a u xxx ,
(37)

11

where δ is an arbitrary dimensionless parameter and V(u, u x . . .) is an arbitrary function. We work
in units where a = 1 and assume the advecting velocity V(u, u x . . .) to be translation-invariant, so
that constants are generically solutions. We look for all V and δ for which (37) admits every asymptotically decaying KdV travelling wave as a solution. We proceed by supposing that every decaying
travelling wave u(x, t) = f (x − ct) ≡ f (ξ) that solves (1) also solves (37) for the same speed c. Then
− c f 0 + 6 f f 0 + f 000 = 0

− c f 0 + V f 0 − δ f 000 = 0.

and

(38)

Eliminating f 000 we get
{V − (δ + 1)c + 6δ f } f 0 = 0.

(39)

This equation for V and δ must be satisfied for arbitrary decaying KdV travelling waves. Since KdV
admits non-constant travelling waves, f 0 . 0. So we must have V = (δ + 1)c − 6δ f . c is eliminated
using the KdV equation −c f 0 + 6 f f 0 + f 000 = 0. However, simply substituting c = 6 f + f 000 / f 0 gives
V1 = 6u + (δ + 1) uu3xx leading back to the KdV equation. On the other hand, we could integrate either
once or twice while omitting integration constants for asymptotically decaying waves and find
− c f + 3 f 2 + f 00 = 0

−

or

1
c 2
f + f 3 + ( f 0 )2 = 0.
2
2

(40)

Eliminating c from (39), we find δ is arbitrary and V can take the functional forms
V2 = 3(1 − δ)u + (1 + δ)

u xx
u

or

V3 = 2(1 − 2δ)u + (1 + δ)

u2x
.
u2

(41)

However, these are not the only possibilities. Instead of integrating −c f 0 + 6 f f 0 + f 000 = 0, we could
differentiate this expression (any number of times!) and then use it to eliminate c. In this manner we
get an infinite sequence of higher order advecting velocities for which (37) admits (2) as a solution:
!
!
6u2x u4x
18u xx u x u5x
V4 = 6u + (1 + δ)
+
,
V5 = 6u + (1 + δ)
+
,
u xx u2x
u3x
u3x
!
24u3x u x 18u2xx u6x
V6 = 6u + (1 + δ)
+
+
, ...
(42)
u4x
u4x
u4x
V2 and V3 are distinguished in that they are linear combinations of 6u (KdV) and new scale-invariant
u2
advecting velocities uuxx or u2x . Vn≥4 differ from KdV by non-scale-invariant advecting velocities.
They also involve higher order derivatives than in the dispersive term u xxx .
Thus, KdV is at the center of an infinite dimensional space of advection-dispersion equations with
sech 2 solitary wave solutions. In the space of equations, we can go out a ‘distance’ 1 + δ in any of
the directions defined by Vn while retaining this solution. In the simplest (least non-linear) case V2 ,
we get the KdV-SIdV family interpolating between KdV (δ = −1) and (3) (δ = 1)

u xx 
ut + 3(1 − δ)u + (δ + 1)
u x = δ u xxx .
(43)
u
It may be noted that δ and  are distinctly different parameters and we have therefore kept a separate
notation.  is a dispersion coefficient that is exclusive to SIdV advection and δ is a different parameter
that allows the KdV and the SIdV(when  = 1) equations
R to be regarded as part of the SIdV-KdV one
parameter family. Like KdV and SIdV, (43) preserves u2 dx since it can be written in conservation
law form
!
!
1  2
1
3
2
u + (1 − δ)u +
+ δ u x − δuu xx = 0.
(44)
t
2
2
x
12

As we omitted a constant of integration in (40), the above argument does not guarantee that every
other KdV travelling wave is a solution of (43). For example, there are non-decaying KdV travelling
waves (e.g. 2sech 2 (x + 2t) − 1) which are not solutions of (43) when δ = 1. Nevertheless, (3) has
cnoidal and non-decaying sech 2 solitary waves of its own. Had we kept
the constant of integration
f 00
2
00
−c f + 3 f + f = A, we would have found V = 3(1 − δ) f + (δ + 1) f − Af . By varying A we get
equations that share particular asymptotic classes of travelling wave solutions with KdV.

7

Discussion and conclusion

In this paper we studied a scale-invariant analogue of the KdV equation, ut + 2u xx u x /u = u xxx ,
which we named the SIdV equation. SIdV is one of the two simplest translation, scale and spacetime parity-invariant non-linear advection-dispersion equations12 . The dimensionless parameter 
measures the strength of dispersion relative to advection. For  = 1, SIdV shares the sech 2 solitary
wave with KdV, as originally discovered by genetic programming. When  = ±2/3, 0, the equation
reduces to the integrable mKdV equation, to a linear dispersive wave equation and to a non-linear
diffusion equation. Like KdV, SIdV admits similarity solutions, solitary and cnoidal travelling waves,
and remarkably, even plane waves. SIdV may be written in conservation law form in two ways (5, 6),
leading to two integrals of motions. In general, we have not been able to find a Lagrangian for SIdV,
indeed we could show that there is no Lagrangian polynomial in φ, where φ x = log u. However, we
found a Lagrangian in a different variable at  = −2/3 and a Hamiltonian when  = 2/3.
For  = 1, numerical evolution of a pair of solitary waves did not show KdV-like soliton scattering, nor did the wave intensity tend towards equipartition among all Fourier modes, even after three
dozen collisions. The wave intensity appeared to circulate among the modes present in the initial
state. Effectively, a finite number of degrees of freedom appeared to take part in the dynamics. This
was also manifested in the boundedness of the gradient energy Q in our simulations. So despite
the apparent absence of soliton scattering and presence of only two known conserved quantities, (3)
seems to display ‘Birkhoff recurrence’ like the Fermi-Pasta-Ulam or KdV systems. Thus, it may
provide a counter example to the idea that integrability is necessary for recurrence. This was suggested in [8], where the number of effective degrees of freedom of a conservative nonlinear wave
equation was identified as the possible origin of a Birkhoff recurrence, a concept generalizing the
well-known Poincaré recurrence of finite-dimensional Hamiltonian systems. The above-mentioned
SIdV-like equation with reciprocal advecting
R velocity V = u/u xx is also interesting from this viewpoint, as it admits the conserved quantity u2x dx. It follows by Poincare’s inequality that the mean
square mode number Q is bounded, and we expect it too to display recurrence!
There are of course several other interesting directions such as perturbation theory for (43) around
the KdV limit, and for SIdV around the solvable  = ±2/3 limits. Additional conserved quantities
and solutions on bounded domains with homogeneous boundary conditions or on a half-line are also
of interest. More generally, we would like to study structural issues like regularity and stability. We
are also examining multi-dimensional analogues of SIdV as well as a complex version (iψt + ψψxx∗ ψ x =
iψ xxx ) which has a locally conserved ‘probability’ density |ψ|2 .
Motivated by the computer-aided discovery of (3) we searched for other advection-dispersion
equations ut +Vu x = δu xxx sharing the KdV solitary wave. We found that there is an infinite sequence
of 1-parameter families of advective velocities Vn (u; δ) (41, 42) that generalize the KdV equation
12

The other one involves the reciprocal of SIdV’s advecting velocity: ut +

13

u
u
u xx x

= ±u xxx .

while retaining its sech 2 solitary wave solution. They involve a linear combination of KdV-like
V ∝ u and new rational advective velocities. The first two of these have the feature of being scaleinvariant: V ∝ u xx /u (SIdV) and u2x /u2 . Among all these equations sharing KdV solitary waves, the
KdV-SIdV family (43) appears to be unique in not involving higher order derivatives and admitting a
positive definite conserved density (u2 ). It would be interesting to study this generalization in greater
detail.
Acknowledgements: The work of GSK was supported by a Ramanujan grant of the Dept. of Science
and Technology, Govt. of India.

A

Computer-aided discovery of SIdV equation

The basic idea behind Genetic Programming (GP) is to simulate a stochastic process by which genetic
traits evolve in offspring, through a random combination of the genes of the parents. Following the
seminal work by Koza [15], the GP framework provides a very useful stochastic engine to discover
various solution regimes in a complex search terrain of a given problem. It is known that an evolutionary method is especially useful when direct methods are not available. In order to set up a GP
engine, a non-linear chromosome structure representing a candidate solution is set up that can potentially grow to a true solution by successive applications of GP operators of selection, copy, crossover
and mutation. The quality of a given chromosome is defined and scaled down typically to a fitness
range [0, 1] with fitness 1 signifying a true solution. Stochastically generated chromosomes fill an
initial pool that is evolved through successive generations in which potentially strong candidate chromosomes are selected based on their fitness values. They undergo possible refinements through GP
operators and hopefully march towards a solution with fitness=1. In using GP to deduce a PDE like
the KdV (in symbolic form) from an analytic travelling wave solution, one begins by considering a
general expression for a third order ODE,

f 000 (ξ) = C ξ, f (ξ), f 0 (ξ), f 00 (ξ)
(45)
where f (ξ) is a function of the travelling wave phase ξ = x − ct . For example, a chromosome during
a GP iteration could be C = 1.1 f 00 + 2 f 0 ( f 0 − 3 f ) + ξ . The fitness parameter is then estimated at
each stage by examining the mean-squared difference between the ‘chromosomal value’ of f 000 and
its value at the given analytic solution. GP follows a ‘fitness driven evolution path’ by minimizing the
error in admitting the given function as a solution. We carried out a number of GP experiments for
the sech2 KdV solitary wave (2). Due to a stochastic search procedure adopted by GP, it was found to
be too slow. We improved and accelerated it by introducing a sniffer technique [12] that carried out a
local search at regular intervals to enhance the minimization procedure. Our improved GP approach
was quite successful in inferring PDEs. Starting with the KdV solitary wave and its derivatives, the
method not only reproduces the KdV equation, but also gives the SIdV equation (3).

References
[1] Korteweg, D. J. and de Vries, G. 1895 On the change of form of long waves advancing in a rectangular
canal, and on a new type of long stationary waves. Philosophical Magazine 39, 422–443.
[2] Scott Russell, J. 1844 Report on waves. Fourteenth meeting of the British Association for the Advancement
of Science.

14

[3] Zabusky, N. J. and Kruskal, M. D. 1965 Interaction of solitons in a collisionless plasma and the recurrence
of initial states. Phys. Rev. Lett. 15, 240–243.
[4] Washimi, H. and Taniuti, T. 1966 Propagation of ion acoustic solitary waves of small amplitude. Phys.
Rev. Lett. 17, 996–998.
[5] Benney, D. J. 1966 Long nonlinear waves in fluid flows. J. Math. Phys., 45, 52–69.
[6] Yomosa, S. 1987 Solitary Waves in Large Blood Vessels. Journal of the Physical society of Japan 56,
506–520.
[7] Fermi, E., Pasta, J., and Ulam, S. 1955 Studies of nonlinear problems. Reprinted in Lectures in applied
mathematics 15, American Math. Society, New York.
[8] Thyagaraja, A. 1979 Recurrent motions in certain continuum dynamical systems. Phys. Fluids 22 (11),
2093–2096.
[9] Thyagaraja, A. 1983 Recurrence phenomena and the number of effective degrees of freedom in nonlinear
wave motion. Nonlinear waves Ch. 17, Ed. L. Debnath, Camb. Univ. Press, 308–325.
[10] Miura, R. M., Gardner, C. S., Kruskal, M. D. 1968 Korteweg-de Vries equation and generalizations. II.
Existence of conservation laws and constants of motion. J. Math. Phys. 9, 1204–1209.
[11] Miura, R. M. 1976 The Korteweg-de Vries equation: a survey of results. SIAM Rev. 18, 412–59.
[12] Ahalpara, D.P. and Sen, A. 2011 A Sniffer Technique for an Efficient Deduction of Model Dynamical
Equations using Genetic Programming. Procs. 14th European Conference, EuroGP 2011(eds: S. Silva et
al.), Torino, Italy, April 2011, Lecture Notes in Computer Science (LNCS 6621), 1–12.
[13] Kudryashov, N.A. 2004 Nonlinear differential equations with exact solutions expressed via the Weierstrass function. Z. Naturforsch. 59a, 443 - 454.
[14] Drazin, P. G. and Johnson, R.S. 1996 Solitons: an introduction. Cambridge Texts in Applied Mathematics.
[15] Koza J. R. 1992 Genetic Programming: On the Programming of Computers by Means of Natural selection
and Genetics, Cambridge, MA: MIT Press.
[16] Kudryashov, N.A., Sinelshchikov, D.I. (2010) Nonlinear waves in bubbly liquids with consideration for
viscosity and heat transfer. Physics LettersA 374 2011-2016.
[17] Randrüüt, M. 2011 On the Kudryashov-Sinelshchikov Equation for waves in bubbly liquids. Physics
Letters A 375 3687-3692.
[18] Ryabov, P.N. 2010 Exact solutions of the Kudryashov-Sinelshchikov equation. Applied Mathematics and
Computation 217 3585-3590.

15

