Compacton solutions and (non)integrability for nonlinear evolutionary PDEs associated with a chain of prestressed granules

TL;DR

This paper studies a nonlinear PDE related to pre-stressed granular chains, exploring its solitary wave solutions, stability, and integrability, with numerical simulations comparing continuous and discrete models.

Contribution

It introduces a family of PDEs supporting compacton solutions derived from granular chains and analyzes their stability, integrability, and numerical behavior.

Findings

PDE supports bright and dark compactons.

Numerical simulations show stable propagation and interactions.

Comparison with discrete model confirms continuous approximation validity.

Abstract

We present the results of study of a nonlinear evolutionary PDE (more precisely, a one-parameter family of PDEs) associated with the chain of pre-stressed granules. The PDE in question supports solitary waves of compression and rarefaction (bright and dark compactons) and can be written in Hamiltonian form. We investigate {\em inter alia} integrability properties of this PDE and its generalized symmetries and conservation laws. For the compacton solutions we perform a stability test followed by the numerical study. In particular, we simulate the temporal evolution of a single compacton, and the interactions of compacton pairs. The results of numerical simulations performed for our model are compared with the numerical evolution of corresponding Cauchy data for the discrete model of chain of pre-stressed elastic granules.

Equations

$$K(m,n):u_t+{\left(u^m\right)}_x+{\left(u^n\right)}_{xxx}=0,m\ge 2,n\ge 2.$$

$${\ddot{Q}}_k(t)=F(Q_{k-1}-Q_k)-F(Q_k-Q_{k+1}),k\in \left\{0,\pm 1,\pm 2\dots .\right\}$$

$$F(z)=A{z}^n,n>1.$$

$$Q_k(t)=u(t,k\cdot a)\approx u(t,x),$$

$$\begin{array}[]{rcl}Q_{k\pm 1}=u(t,x\pm a)&=&\exp(\pm aD_{x})u(t,x)\\ &=&\displaystyle\sum_{j=0}^{4}\frac{(\pm a)^{j}}{j!}\frac{\partial^{j}}{\partial x^{j}}u(t,x)+O\left(a^{5}\right),\end{array}$$

$$u_{tt}=-C{\left\{{\left(-u_x\right)}^n+\beta {\left(-u_x\right)}^{\frac{n-1}{2}}{\left[{\left(-u_x\right)}^{\frac{n+1}{2}}\right]}_{xx}\right\}}_x,$$

$$C=A{a}^{n+1},\beta =\frac{n{a}^2}{6(n+1)}.$$

$$S_{tt}=C{\left\{S^n+\beta S^{\frac{n-1}{2}}{\left[S^{\frac{n+1}{2}}\right]}_{xx}\right\}}_{xx}.$$

$$F(z)=A{z}^n+Bz.$$

$$u_{tt}=-C{\left\{{\left(-u_x\right)}^n+\beta {\left(-u_x\right)}^{\frac{n-1}{2}}{\left[{\left(-u_x\right)}^{\frac{n+1}{2}}\right]}_{xx}\right\}}_x-\gamma a^{n+3}{\left(-u_x\right)}_x.$$

$$S_{tt}=C{\left\{S^n+\beta S^{\frac{n-1}{2}}{\left[S^{\frac{n+1}{2}}\right]}_{xx}\right\}}_{xx}+\gamma a^{n+3}{S}_{xx}.$$

$$S_{\tau \tau }=\frac{C}{\gamma a^{n+3}}{\left\{S^n+\beta S^{\frac{n-1}{2}}{\left[S^{\frac{n+1}{2}}\right]}_{xx}\right\}}_{xx}+S_{xx}.$$

$$W_{\bar{T}}+\frac{A}{\gamma }{\left\{W^n+\frac{n}{6(n+1)}{W}^{\frac{n-1}{2}}{\left[W^{\frac{n+1}{2}}\right]}_{\xi \xi }\right\}}_{\xi }=0.$$

$$t=\frac{A}{\gamma }L\bar{T},x=L\xi ,$$

$$W_T+{\left\{W^n+W^{\frac{n-1}{2}}{\left[W^{\frac{n+1}{2}}\right]}_{XX}\right\}}_X=0,$$

$$F(z)=-A{z}^{2k}+Bz$$

$$W_T-{\left\{W^n+W^{\frac{n-1}{2}}{\left[W^{\frac{n+1}{2}}\right]}_{XX}\right\}}_X=0,n=2k.$$

$$W_T+{\left[\mathrm{sgn}(W)\right]}^{n+1}{\left\{W^n+W^{\frac{n-1}{2}}{\left[W^{\frac{n+1}{2}}\right]}_{XX}\right\}}_X=0.$$

$$W_T=-{\left(W^n+W^{(n-1)/2}{\left[W^{(n+1)/2}\right]}_{XX}\right)}_X$$

$$W_T=D_X\delta {\mathcal{H}}_{\mathrm{Nest}}/\delta W\equiv F.$$

$$h_{\mathrm{Nest}}=\left\{\begin{array}[]{l}\displaystyle\left(\frac{1}{4}(n+1)W^{n-1}W_{X}^{2}-W^{n+1}/(n+1)\right)\quad\mbox{for $n\neq-1$},\\[14.22636pt] \displaystyle\ln|W|\quad\mbox{for $n=-1$.}\end{array}\right.$$

$$D_X=\frac{\partial }{\partial X}+\sum _{j=0}^{\infty }{W}_{j+1}\frac{\partial }{\partial W_j},D_T=\frac{\partial }{\partial T}+\sum _{j=0}^{\infty }{D}_X^j(F)\frac{\partial }{\partial W_j}.$$

$$\frac{\delta \mathcal{F}}{\delta W}=\sum _{j=0}^{\infty }(-D_X{)}^j\left(\frac{\partial f}{\partial W_j}\right).$$

$$f_{\ast }=\sum _{j=0}^k\frac{\partial f}{\partial W_j}{D}_X^j.$$

$$D_T(\rho )=D_X(\sigma ),$$

$$E_W=\sum _{j=0}^{\infty }(-D_X{)}^j\circ \frac{\partial }{\partial W_j}.$$

$${\rho }_0=W,{\rho }_1=W^2/2,{\rho }_2=h_{\mathrm{Nest}}.$$

$${\rho }_3=X^{2}W{W}_X-T{W}_X^2/W.$$

$$W_T=K(X,W,W_X,\dots ,{\partial }^{k}W/\partial X^k)$$

$$D_T(G)=K_{\ast }(G),$$

Full text

For the compacton solutions we perform a stability test followed by the numerical study. In particular, we simulate the temporal evolution of a single compacton, and the interactions of compacton pairs. The results of numerical simulations performed for our model are compared with the numerical evolution of corresponding Cauchy data for the discrete model of chain of pre-stressed elastic granules.

Keywords: chains of pre-stressed granules; compactons; integrable systems; symmetry integrability; symmetries; conservation laws; stability test; conserved quantities; Hamiltonian structures; numerical simulation

MSC 2010 35B36; 74J35; 74H15; 37K05; 37K10

1 Introduction

This paper deals with nonlinear evolutionary PDEs associated with dynamics of a one-dimensional chain of pre-stressed granules which arises in quite a number of applications. Since Nesterenko’s pioneering works [1, 2] propagation of pulses in such media has been a subject of a great number of experimental studies and numerical simulations, see [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] and references therein. We consider a nonlinear evolutionary PDE which is derived from the infinite system of ODEs describing the dynamics of one-dimensional chain of elastic bodies interacting with each other by means of a nonlinear force. The PDE in question is obtained through the passage to continuum limit followed by the formal multi-scale decomposition.

The PDE under study turns out to admit a Hamiltonian representation and possess localized traveling wave solutions manifesting some features of solitons. For this reason, it is of interest to investigate its complete integrability. We do this below along with the study of generalized symmetries and conservation laws. We show below that the compacton traveling wave (for traveling waves in general see e.g. [14, 15] and references therein) solutions satisfy the necessary condition for the extremum of a functional associated with the Hamiltonian. Using this we also perform a stability test followed by the numerical study of the compacton solutions. Somewhat surprisingly, numerical simulations show that even in a nonintegrable case the compacton solutions recover their shapes after the collisions, yet the dynamics of interaction slightly differs from that of KdV solitons. In this connection note that compactons, i.e. soliton-like solutions with compact support, see [1, 3, 16] and references therein, exist for a number of physically relevant models and possess several interesting features making them a subject of intense research, cf. e.g. [17, 18, 19, 20, 21, 22] and references therein.

The paper is organized as follows. In section 2 we introduce the continual analog of the granular pre-stressed media with the specific interaction of the adjacent blocks which allows for the description of both the waves of compression and rarefaction. In section 3 we present the Hamiltonian structure of the equation in question. In section 4 we study the conservation laws admitted by the said equation. In section 5 we perform the integrability test that singled out an exceptional integrable case, which is studied in more detail in section 6. In section 7 we show that the compacton traveling wave (TW) solutions that satisfy factorized equations also satisfy necessary conditions of extrema for the appropriate Lagrange functionals. Next we perform stability tests for compacton solutions based on the approach developed in [23, 24, 25], and show that both dark and bright compactons pass the stability test. The results of qualitative analysis are backed and partly supplemented by the numerical study performed in section  8. We also present the results of numerical simulation of the Cauchy problem for discrete chains and compare the results obtained with the analogous simulations performed for the continual analogue of these chains. The closing section 9 contains conclusions and discussion.

2 Evolutionary PDEs associated with the granular prestressed chains

Amazing features of the solitons associated with the celebrated Korteweg–de Vries (KdV) equation, as well as other completely integrable models [14], are often ascribed to the existence of higher symmetries and infinite sets of conservation laws, cf. e.g. [26, 27, 28]. However, there exist non-integrable equations possessing localized TW solutions with quite similar behavior. A well-known example of this is provided by the $K(m,n)$ equations [16]:

[TABLE]

The members of this hierarchy are not completely integrable at least for generic values of the parameters $m$, $n$, see [17, 20] and references therein, and yet possess compactly-supported TW solutions exhibiting solitonic features [16, 29].

The $K(m,n)$ family was introduced in the 1990s as a formal generalization of the KdV hierarchy without referring to its physical context. Earlier V.F. Nesterenko [1] considered the dynamics of a chain of preloaded granules described by the following ODE system:

[TABLE]

where $Q_{k}(t)$ is the displacement of the $k$th granule center-of-mass from its equilibrium position,

[TABLE]

He has described for the first time the formation of localized wave patterns and evolution within this model [1, 2, 3]. In [1, 2, 5] he presented the nonlinear evolutionary PDEs being the quasi-continual limits of the discrete models; in this connection cf. also [30].

The transition to the continual model is achieved via the substitution

[TABLE]

where $a$ is the average distance between granules. Insert this formula, together with the substitutions

[TABLE]

into (2), and observe that the term of lowest order in $a$ on the right-hand side of (2) is proportional to $a^{n+1}$. Now expanding the right-hand side of (2) divided by $a^{n+1}$ into the (formal) Taylor series and then dropping the terms of the order $O(a^{4})$ and higher in this expansion yields from (2) the equation

[TABLE]

where

[TABLE]

Differentiating the above equation with respect to $x$ and employing the new variable $S=\left(-u_{x}\right)$ corresponding to the strain field, one obtains the Nesterenko equation [5]:

[TABLE]

Eq. (6) was derived using only one small parameter corresponding to the long wave approximation. Thus it can describe the dynamics of “strongly preloaded media” with dynamic amplitude much smaller than the preload or the dynamics of “weakly preloaded media” when the dynamic amplitude in the wave is much larger than the preload or even when the preload is equal to zero, in which case the propagation of acoustic waves is impossible (the effect of “sonic vacuum” [2]). As it is shown in [5], equation (6) possesses a one-parameter family of compacton TW solutions describing the propagation of the waves of compression.

Unfortunately, the compacton solutions supported by (6) are unstable. This can be verified by a direct numerical calculation, substituting in the corresponding difference scheme as Cauchy data known compacton solutions.

A similar situation occurs in the case of the Boussinesq equation, obtained as a continuum limit of the Fermi–Pasta–Ulam system of coupled oscillators [14]. As is well known, the Boussinesq equation possesses unstable soliton-like solutions, and the KdV equation, supporting the stable uni-directional solitons, is extracted from the Boussinesq equation by means of the asymptotic multi-scale expansion [14], cf. also [31] and references therein.

In this connection it should be also noted that the instability caused by the short wavelengths can be removed using the regularization consisting in replacing the space derivatives of the force by mixed space and time derivatives. The regularized equation for the case of general power law is nothing but Eq. (1.110) from [5], and its counterpart for a general interaction law is Eq. (1.156) from [5].

It should be further noted that, at least for $n=3/2$, equation (6) has stationary compacton solutions which are close to the numerical solutions of the discrete Hertzian chain, see e.g. [8], and the numerical simulations and experiments strongly suggest that the latter solutions are stable. For example, such solutions are generated from various initial conditions on short distances from the disturbed end and propagate in experiments despite disturbances due to inevitable dissipation and violation of periodicity, see e.g. experimental results in [5].

Another interesting observation is that the conditions for existence of solitary waves in discrete chain [32] and in the continuum approximation, see Eq.(1.154) at p. 108 in [5], are identical and based on the sign of the second derivative of the force, see p. 113 in [5].

Our approach to finding a “proper” compacton-supporting equation is as follows. We start from the discrete system (2) in which the interaction force has the form

[TABLE]

In addition, we assume that $B=\gamma a^{n+3}$, where $|\gamma|=O(|A|)$.

The interaction law in (7) is a special case of general interaction law that results in long wave equation, Eq.(1.154) at p. 108 in [5] or its regularized counterpart, Eq.(1.156) in [5]. The stationary solutions of the said long wave equation are studied in [5], where, depending on the behavior of the second derivative of the interaction law, strongly nonlinear compression or rarefaction solitary waves are predicted.

In this connection also note that for small deformations the interaction law in (7) is a special case of the situation where the first derivative is nonzero and higher derivatives are zero except for the $n$-th order one which, in conjunction with the discussion in the preceding paragraph implies, cf. [5], in particular p.110–123, that the linear part of the interaction law is, to an extent, irrelevant for the study of qualitative behavior of stationary solutions.

Inserting (4), (5) into the formula (2) and assuming that the interaction is described by (7), we obtain, up to the terms of the order $O(a^{4})$ and higher in the expansion of the right-hand side of (2) divided by $a^{n+1}$, cf. the discussion after (5), the equation

[TABLE]

Differentiating the above equation with respect to $x$ and introducing the new variable $S=\left(-u_{x}\right)$, we obtain the following equation:

[TABLE]

Now we use a series of scaling transformations. Employing the scaling $\tau=\sqrt{\gamma a^{n+3}}t$ enables us to rewrite the above equation in the form

[TABLE]

Next, the transformation $\bar{T}=\frac{1}{2}a^{q}\tau$, $\xi=a^{p}(x-\tau)$, $S=a^{r}W$ is used. If, for example, we assign the following values to the parameters $q=1$, $p=-1$, $r=5/n,$ then the higher-order coefficient $O(a^{2})$ will be that of the second derivative with respect to $\bar{T}$. So, dropping the terms proportional to $O(a^{2})$, we obtain, after the integration with respect to $\xi$, the equation:

[TABLE]

Performing the rescaling and returning to the initial notation

[TABLE]

where $L=\sqrt{\frac{6(n+1)}{n}}$, we finally obtain the sought-for equation

[TABLE]

to which we shall hereinafter refer as to the Nesterenko equation. Note that Eq.(9) appears in [19] (see also [22]) as a particular case of the $C_{1}(m,a+b)$ hierarchy introduced as a generalization of the set of $K(m,n)$ equations.

The description of waves of rarefaction in the case $n=2k$ requires the following modification of the interaction force:

[TABLE]

(for $n=2k+1$ the formula (7) describes automatically both waves of compression and of rarefaction). Applying the above machinery to (2) with the interaction (10), we obtain, in the same notation, the equation

[TABLE]

Thus, the universal equation describing waves of compression and rarefaction for arbitrary $n\in\mathbb{N}$ can be written in the form

[TABLE]

In closing note that equations (9), (11) and (12) are obtained by formal application of the multiscale decomposition method which cannot be substantiated in our case because of negativity of the index $p$, cf. [33] where this problem is discussed in a more general fashion. Further study of these equations is justified by the fact that they possess a set of compacton solutions possessing interesting dynamical features. As will be shown below, these solutions describe well enough propagation of short impulses in the chain of pre-stressed blocks.

3 Hamiltonian structure for the Nesterenko equation

Now return to (9) which we now write in the manifestly evolutionary form, that is,

[TABLE]

Note that for $n=-1$ this equation boils down to a quasilinear first-order equation $W_{T}=(W^{-1})_{X}$ which is obviously integrable, and for $n=1$ equation (13) becomes linear.

Equation (13) can be written (cf. [22]) as

[TABLE]

Thus, (14) is written in Hamiltonian form with the Hamiltonian $\mathcal{H}_{\mathrm{Nest}}$ and the Hamiltonian structure $\mathfrak{P}_{0}=D_{X}$.

This implies, in particular, that to any nontrivial local conserved density of (14) there corresponds a (generalized, but not necessarily genuinely generalized (see the definition below), and possibly trivial) symmetry of (14).

Here $\delta/\delta W$ is the variational derivative (see below for details) and $\mathcal{H}_{\mathrm{Nest}}=\displaystyle\int h_{\mathrm{Nest}}dX$ with the density

[TABLE]

Here and below the integrals are understood in the sense of formal calculus of variations, see e.g. [28, 34]. Here we put, cf. [27, 28, 34], $W_{j}=\partial^{j}W/\partial X^{j}$, $j=1,2,\dots$, $W_{0}\equiv W$, and define [26, 27, 28, 34] the total derivatives

[TABLE]

The variational derivative of a functional $\mathcal{F}=\displaystyle\int f(X,T,W,W_{1},\dots,W_{k})dX$ has the form

[TABLE]

For any $f=f(X,T,W,\dots,W_{k})$ we also define, cf. e.g. [27, 28], its linearization

[TABLE]

4 Conservation laws

Recall, cf. e.g. [26, 27, 28, 34, 35, 36, 37, 38] and references therein, that a local conservation law for (13) is, roughly speaking, a relation of the form

[TABLE]

where $\rho=\rho(X,T,W,W_{1},\dots,W_{r})$ and $\sigma=\sigma(X,T,W,W_{1},\dots,W_{s})$, which holds by virtue of (13). Here $\rho$ and $\sigma$ are called a (conserved) density and the flux of our conservation law.

Also recall, cf. e.g. [36], that a conservation law (18) is called nontrivial if there exists no function $\zeta(X,T,W,W_{1},\dots,W_{q})$ such that $\rho=D_{X}\zeta$, i.e., $\rho\not\in\mathrm{Im}D_{X}$.

It is well known, see e.g. [28, 34], that a necessary and sufficient condition for a function $f=f(X,T,W,W_{1},\dots,W_{r})$ to not belong to the image of $D_{X}$ is $E_{W}f\neq 0$, where $E_{W}$ is the Euler operator

[TABLE]

Hence $\rho$ is a conserved density for (13) if and only if $E_{W}D_{T}(\rho)=0$, and this density is nontrivial if and only if $E_{W}\rho\neq 0$.

It is readily checked that we have the following

Proposition 1**.**

For any $n$ equation (13) admits the following three conserved densities:

[TABLE]

For $n=0$ we have an extra density

[TABLE]

Moreover, for $n\neq 0,1,-1,-2$ (resp. for $n=0$) the densities (19) (resp. (19 and (20)) exhaust, modulo the addition of trivial ones, the linearly independent conserved densities of order up to five, i.e., of the form $\rho=\rho(X,T,W,W_{X},\dots,W_{XXXXX})$.

It is very likely that for $n\neq 1,-1,-2$ no local conserved densities of order greater than five (of course, again modulo trivial ones) exist at all in view of nonintegrability of (13) for $n\neq 1,-1,-2$ as discussed below.

Recall that $\rho_{2}$ is the density of the Hamiltonian $\mathcal{H}_{\mathrm{Nest}}$ for (13) with respect to the Hamiltonian structure $\mathfrak{P}_{0}=D_{X}$. To the functional $\mathcal{C}=\displaystyle\int WdX$ there corresponds a trivial symmetry, i.e., a symmetry with zero characteristic, as $D_{X}\delta\mathcal{C}/\delta W=0$, so $\mathcal{C}$ is a Casimir functional for $\mathfrak{P}_{0}$. To the functional $\mathcal{P}=\frac{1}{2}\displaystyle\int W^{2}dX$ there corresponds a symmetry with the characteristic $W_{X}=D_{X}\delta\mathcal{P}/\delta W$, that is, $X$-translation, and to $\mathcal{H}_{\mathrm{Nest}}$ there corresponds a symmetry with the characteristic equal to the r.h.s. $F$ of (14), i.e., the time translation symmetry.

For $n=0$ to the conserved functional $\mathcal{H}_{3}=\displaystyle\int\rho_{3}dX$ there corresponds a scaling symmetry with the characteristic $4TF+2XW_{X}+2W=D_{X}\delta\mathcal{H}_{3}/\delta W$. Again, it is very likely that $\rho_{i}$, $i=0,\dots,3$, are the only local conserved densities (modulo trivial ones) for (13) with $n=0$ in view of nonintegrability of this special case of (13).

5 Integrability

Integrable equations of the form (14) with the Hamiltonian of general form $\mathcal{H}=\displaystyle\int dXh(W,W_{X})$ where the density $h=h(W,W_{X})$ is such that $\partial^{2}h/\partial W_{X}^{2}\neq 0$ were classified (modulo point transformations leaving $T$ invariant) in [39]. Note that in [27, 39] and references therein integrability of an evolution equation

[TABLE]

with $k\geq 2$ means existence of an infinite hierarchy of generalized symmetries of increasing orders which do not depend explicitly on $T$. In order to avoid ambiguity we shall, following the common usage, refer below to this kind of integrability as to the symmetry integrability.

Recall, cf. e.g. [26, 27, 28, 34], that a generalized symmetry of order $r$ for (21) is111For the sake of simplicity and without loss of generality we identify here a generalized symmetry with its characteristic. a function $G=G(X,T,W,W_{1},\dots,\allowbreak W_{r})$ such that $\partial G/\partial W_{r}\neq 0$ and

[TABLE]

where now $D_{T}=\displaystyle\frac{\partial}{\partial T}+\sum\limits_{j=0}^{\infty}D_{X}^{j}(K)\displaystyle\frac{\partial}{\partial W_{j}}$.

Such a symmetry $G$ is known as genuinely generalized if it cannot be written in the form $G=c(T)K+b(X,T,W,W_{X})$ for some functions $b$ and $c$, that is, it is not equivalent to a point or contact symmetry. As far as point symmetries of the equations studied in the present paper, and, more broadly, of $\mathcal{C}_{1}(m,a,b)$ equations (see e.g. [19, 22]), cf. e.g. [40] and references therein.

Thus, symmetry integrability of (21) means existence of an infinite hierarchy of generalized symmetries of the form $G_{i}(X,W,W_{1},\dots,W_{r_{i}})$ of increasing orders $r_{i}$.

Now turn to comparison of the density $h_{\mathrm{Nest}}$ of our Hamiltonian and the densities $h$ found in [39] for which the equation $W_{T}=D_{X}(\delta\mathcal{H}/\delta W)$ with the general Hamiltonian $\mathcal{H}=\displaystyle\int h(W,W_{X})dX$ is symmetry integrable.

Proposition 2**.**

The only symmetry integrable case of (13) which is genuinely nonlinear and genuinely of third order is that of $n=-2$.

Proof. It is not difficult to observe (cf. e.g. [41]) that using point transformations leaving $t$ invariant the density $h_{\mathrm{Nest}}$ of our Hamiltonian for $n\neq-1$ can, if at all, only be transformed into just one case from [39], namely, equation (2.1) in [39], that is,

[TABLE]

where $f=c_{0}+c_{1}W+c_{2}W^{2}$, $P=\sum\limits_{i=0}^{4}d_{i}W^{i}$, and $c_{i}$ and $d_{i}$ are arbitrary constants.

Moreover, it is clear that in our case $f$ should actually be a monomial: $f=cW^{\alpha}$, $\alpha=0,1,2$.

Upon comparing the coefficients at $W_{X}^{2}$ in (23) and (15) modulo an obvious rescaling of $W$, we see that all values of $n$ for which (13) could be integrable should satisfy $n-1=0,-3,-6$. The case of $n=1$ is trivially integrable, as then (13) is just a linear equation, so we are left with two possibilities $n=-2$ and $n=-5$ corresponding to $\alpha=1$ and $\alpha=2$.

Now upon inspecting the remaining terms in $h_{\mathrm{Nest}}$ and in (23) we readily conclude that the polynomial $P$ should also reduce to a single monomial: $P=dW^{\beta}$, where $\beta=0,1,2,3,4$, so we have a system $n-1=-3\alpha$ and $n+1=\beta-\alpha$, where $\alpha=1,2$ and $\beta=0,1,2,3,4$. An obvious corollary of this system is $-3\alpha+2=\beta-\alpha$, whence $\beta=2(1-\alpha)$. However, $\beta\geq 0$ by assumption, so the case of $n=-5$, when $\alpha=2$ and we should have $\beta=-2$, is not integrable.

Thus, the only integrable case of (13) which is genuinely nonlinear and genuinely of third order is that of $n=-2$, and the result follows. $\Box$.

Recall that for $n=-1$ equation (13) degenerates and becomes a first order quasilinear equation whose general solution can be found, see above, and for $n=1$ equation (13) is just linear.

In fact, the result of Proposition 2 can be further strengthened so that absence of any generalized symmetries, rather than just those that do not depend explicitly on $T$, can be established.

To this end consider, following [27], the so-called canonical density $\rho_{-1}=(\partial F/\partial W_{XXX})^{-1/3}$. It is readily checked that $E_{W}D_{T}(\rho_{-1})\neq 0$ for $n\neq-1,-2,-5,1$. Hence for $n\neq-1,-2,-5,1$ we have $D_{T}(\rho_{-1})\not\in\mathrm{Im}\ D_{X}$, and thus $\rho_{-1}$ is not a density of a local conservation law for (13).

In turn, by virtue of the results from [42] this immediately implies

Proposition 3**.**

Equation (13) for $n\neq 1,-1,-2,-5$ has no generalized symmetries of order greater than three.

In other words, Proposition 3 means that for $n\neq 1,-1,-2,-5$ any solution $G=G(X,T,W,W_{1},\dots,W_{r})$ of the equation

[TABLE]

where $D_{T}$ and $F$ are given in (16) and (14), in fact depends at most on $X,T,W,W_{X},W_{XX},W_{XXX}$.

This implies that (13) for $n\neq-1,-2,-5,1$ admits no genuinely generalized symmetries, and hence (13) for $n\neq-1,-2,-5,1$ is unlikely to be integrable in any reasonable sense, cf.  [27].

Leaving aside the degenerate cases of $n=\pm 1$, turn to the remaining two special cases: $n=-2$ and $n=-5$. We believe that using the technique similar to that of [43] (cf. also [20, 44]) it can be shown that in the case of $n=-5$ equation (13) admits no genuinely generalized symmetries, including those with explicit dependence on $T$ and not just the time-independent ones whose nonexistence follows from the above comparison of (15) with (23), so we are left with just one integrable case of $n=-2$ which we discuss below.

6 Nesterenko equation for $n=-2$: integrability and beyond

The following result is readily checked by straightforward computation:

Proposition 4**.**

For $n=-2$ equation (13) has a Lax pair of the form

[TABLE]

and admits a recursion operator

[TABLE]

The recursion operator (26) can be found e.g. using the technique from [45] (cf. also [46]). Also note that upon passing to a new dependent variable equal to a square of $W$ the first equation of (25) can be identified with a special case of the eigenvalue problem related to the extended Harry Dym systems, see e.g. [47] and references therein.

Equation (13) for $n=-2$ also admits a second local Hamiltonian operator $\mathfrak{P}_{1}=\mathfrak{R}\circ D_{x}$, that is,

[TABLE]

which is compatible with $\mathfrak{P}_{0}=D_{X}$, so the recursion operator $\mathfrak{R}$ is hereditary and equation (13) for $n=-2$ can be written, in addition to (14), in the second Hamiltonian form as

[TABLE]

where $\tilde{h}=W/2$.

Thus, we have the following

Proposition 5**.**

Equation (13) for $n=-2$ is an integrable bihamiltonian system with two local Hamiltonian operators $\mathfrak{P}_{0}$ and $\mathfrak{P}_{1}$ and two local Hamiltonian representations (14) and (27).

Using general theory of bihamiltonian systems (see e.g. [28, Ch. 7] and [48]), we also readily obtain

Corollary 1**.**

Equation (13) for $n=-2$ possesses an infinite hierarchy of commuting generalized symmetries of the form $\mathfrak{R}^{k}W_{X}$, $k=0,1,2,\dots$ and an infinite hierarchy of local conservation laws whose densities $h_{j}$ are generated recursively through the relations

[TABLE]

where $j=0,1,2,\dots$ and $h_{0}=W/2$, and of associated integrals of motion $\mathcal{H}_{j}=\int h_{j}dX$ in involution with respect to the two Poisson brackets associated with $\mathfrak{P}_{0}$ and $\mathfrak{P}_{1}$.

The fact that the generalized symmetries $\mathfrak{R}^{k}W_{X}$ and the conserved densities $h_{k}$ for $k=0,1,2,\dots$ do not involve any nonlocal terms can be established using the results of [49] or [50] (cf. also [46]).

As we have already pointed out above, up to a suitable rescaling of $T$ and obvious change of notation equation (14) for $n=-2$ is a special case of equation (2.1c) in [39], and hence can be transformed into a special case of the well-known $S$-integrable Calogero–Degasperis–Fokas [51, 52] equation in the manner described therein.

Namely, pass first to the potential form of (14) with $n=-2$,

[TABLE]

related to (13) through the differential substitution $W=V_{X}$.

The subsequent hodograph transformation interchanging $X$ and $V$ turns the above equation into a constant separant equation

[TABLE]

or, upon a suitable rescaling of $T$,

[TABLE]

Finally, putting $V=\exp(U/2)$ turns (28) into a special case of the Calogero–Degasperis–Fokas [51, 52] equation, viz.,

[TABLE]

7 Compacton solutions and stability tests

Consider the pair of equations (9), (11), which can be represented by the single expression

[TABLE]

As we are interested in the traveling wave (TW) solutions $W=W(z)\equiv W(X-cT)$, it is convenient to pass to the TW coordinates $T\to T$, $X\to z=X-cT$. This change of variables yields from (30) the equation

[TABLE]

It is easy to check that equation (31) admits a Hamiltonian formulation

[TABLE]

where now

[TABLE]

and

[TABLE]

The above formulation up to the coefficient $\epsilon$ follows directly from the Hamiltonian form (cf. (14)) of equation (9) after the change of coordinates. Recall that both functionals $\mathcal{H}_{\mathrm{Nest}}$ and $\mathcal{P}$ are conserved in time.

Now consider the following functions:

[TABLE]

where $\epsilon=\pm 1$,

[TABLE]

It is readily checked that we have the following

Proposition 6**.**

If $n=2k+1$, $k\in\mathbb{N}$, then the functions $W_{c}^{\pm}(z)$ are weak solutions to the equation

[TABLE]

If $n=2k$, $k\in\mathbb{N}$, then the functions $W_{c}^{\pm 1}(z)$ are weak solutions to the equation

[TABLE]

So, the TW solutions (33) are the critical points of either the Lagrange functional $\Lambda=\mathcal{H}_{\mathrm{Nest}}+\beta\mathcal{P}$ (the case of $n=2k+1$) or $\Lambda^{\epsilon}=\epsilon\mathcal{H}_{\mathrm{Nest}}+\beta\mathcal{P}$ (the case of $n=2k$) with the common Lagrange multiplier $\beta=c$. As is well known, necessary and sufficient condition for $\Lambda$ (resp. $\Lambda^{\epsilon}$) to attain the minimum on the compacton solutions can be stated in terms of the positivity of the second variation of the corresponding functional, which, in turn, guarantees the orbital stability of the TW solution [53]. Here we do not touch upon the problem of strict estimating of the signs of the second variations. We follow instead the approach suggested in [23, 24, 25], which enables us to test the possibility of appearance of the local minimum on selected sets of perturbations of TW solutions.

Consider the following family of perturbations

[TABLE]

Upon choosing $\alpha=1/2$ we obtain

[TABLE]

Thus, for this choice $\mathcal{P}[\lambda]$ keeps its unperturbed value. By imposing this condition we reject “fake” perturbations associated with the translational symmetry $T_{\delta}\left[W_{c}^{\epsilon}(z)\right]=W_{c}^{\epsilon}(z+\delta)$. Indeed, since equations (34), (35) are invariant under the shift $z\to z+\delta$, $T_{\delta}W_{c}^{\epsilon}(z)$ belongs to the set of solutions as well, while formally the transformation $W_{c}^{\epsilon}(z)\rightarrow W_{c}^{\epsilon}(z+\delta)$ can be treated as a perturbation. In order to exclude the perturbations of this sort, the orthogonality condition is imposed. Introducing the representation for the perturbed solution

[TABLE]

and using the condition (37), we find

[TABLE]

so if $\mathcal{P}$ is independent of $\lambda$, then, up to $O\left(||v(z,\lambda)||^{2}\right)$ the perturbation created by the scaling transformation is orthogonal to the TW solution.

For $\alpha=1/2$ and $n\in\mathbb{N}$, we arrive at the following functions to be tested:

[TABLE]

where

[TABLE]

[TABLE]

If the functional $\Lambda^{\nu}=\nu\mathcal{H}_{\mathrm{Nest}}+c\mathcal{P}$ attains the extremal value on the compacton solution, then the function $\Lambda^{\nu}[\lambda]$ has the corresponding extremum in the point $\lambda=1$. The verification of this property is employed as a test.

A necessary condition for the extremum \frac{d}{d\lambda}\Lambda^{\nu}[\lambda]\Bigl{|}_{\lambda=1}=0 gives us the equality

[TABLE]

Using (39), we obtain the estimate

[TABLE]

which is valid for both $n=2k+1$ and $n=2k$. Thus, the generalized solutions (33) pass the test for stability, and we can state the following

Conjecture. For $n>1$ weak solutions (33) provide minima of the functional $\Lambda^{\nu}$.

Further information about the properties of the compacton solutions is provided by the numerical simulations discussed below.

8 Numerical simulations for dynamics of compactons

The dynamics of solitary waves is studied by means of direct numerical simulation based on the finite-difference scheme.

To derive a finite-difference scheme, say, for the model equation (9), we modify the scheme presented in [29]. In agreement with the methodology proposed in that paper we introduce the artificial viscosity by adding the term $\varepsilon W_{4x}$, where $\varepsilon$ is a small parameter. Thus, instead of (9) we have for the case of $n=3$ the following equation:

[TABLE]

Let us approximate the spatial derivatives as follows:

[TABLE]

where $L_{j}=W_{j}\frac{W_{j-2}^{2}-2W_{j}^{2}+W_{j+2}^{2}}{h^{2}}$ .

To integrate the system (41) in time, we use the midpoint method. Then the quantities $W_{j}$ and $\dot{W}_{j}$ are represented in the form

[TABLE]

The resulting nonlinear algebraic system with respect to $W_{j}^{n+1}$ can be solved by iterative methods.

We test the scheme (41) by considering the movement of a single compacton. Assume that the model parameters $c=1$ and the scheme parameters $N=600$, $h=30/N$, $\tau=0.01$, $\varepsilon=10^{-3}$ are fixed. The application of the scheme (41) gives us fig. 1a.

The starting profile providing the initial condition for the numerical scheme is chosen according to (33) where $n=3$, $c_{1}=1$ and $c_{2}=1/4$, namely,

[TABLE]

where $z_{1}=5$, $z_{2}=13$, and $\epsilon=+1$ for fig. 1 while $\epsilon=-1$ for fig. 2 (note that $W_{i}$ corresponds to $i$-th figure).

To study the interaction of two bright compactons, we combine the compacton having the velocity $c=1$ with the slow one characterized by the velocity $c=1/4$ and being shifted to the right at the initial moment of time. The result of modelling is presented at fig. 1b. The interaction of two dark compactons has similar properties and is depicted at fig. 2.

As we have already mentioned at the end of Section 2, there is no way of selecting the scales in the model equations (9), (11) and (12), so the scaling decomposition employed there is rather formal. Nevertheless, it leads to interesting equations possessing localized solutions with solitonic features.

Now we are going to compare the evolution of the compacton solutions with corresponding solutions of the finite (but long enough) discrete system. Since the average distance $a$ between adjacent blocks does not play the role of a small parameter anymore, we assume from now on that it is equal to one. With this assumption in mind, we can write equation (12) in the initial variables $t,x$ as follows:

[TABLE]

where

[TABLE]

It is easy to verify that equation (42) possesses the following compacton solutions:

[TABLE]

where $\epsilon=\pm 1,$ $z=x-ct,$

[TABLE]

We introduce the functions $R_{k}=Q_{k-1}-Q_{k}$ being the discrete analogs to the strain field $W(t,x)$. These functions are assumed to satisfy the system

[TABLE]

We solve this system with the following initial conditions induced by the compacton solution (43) in the respective nodes:

[TABLE]

[TABLE]

[TABLE]

where $I$ is a constant phase, $k=2,3,\dots,m-1$. Note that $A$ and $\gamma$ appear in equation (42) in the form of the ratio $Q=A/\gamma$, whereas in the system (44) they appear as independent parameters. Therefore, one should not expect a one-to-one correspondence between the solutions of the discrete and continuous problems for arbitrary values of the parameters. The numerical experiments confirm this hypothesis by showing that synchronous evolution of the same compacton perturbation within two models can be observed for a unique value of the velocity $c=c_{0}$. This value depends strongly on the parameter $\gamma$ and depends on the parameter $A$ in a much weaker fashion. It has been observed that at $c<c_{0}$ the discrete compacton moves quicker than its continuous analogue while at $c>c_{0}$ the opposite effect occurs. The result of comparison for a single compacton is shown at fig. 3. One can see that at the chosen values of the parameters the main perturbations move synchronously and do not change their form. However, in the tail part of the discrete analogue small nonvanishing oscillations appear after a while.

Since for every value of the parameter $\gamma$ there is a unique value of the wave pack velocity for which the discrete and continuous compacton perturbations move synchronously, one should not expect that the collision of compactons within these two models will proceed in the same way for any set of values of parameters. However, collision processes display not much of qualitative differences for the discrete pulses which interact elastically like their continuous analogues. This is illustrated on fig. 4 showing the evolution of two initially separated discrete compactons. For convenience, the continuous compactons which coincide with the right-hand side of the initial data (45) at $t=0$ (the leftmost graph in the first row) are also shown in this figure. Continuous curves shown on the following graphs are obtained by appropriate translations. They are presented in order to emphasize the quasi-elastic nature of interaction of the discrete pulses.

9 Conclusions and discussion

In the present paper we have studied compacton solutions supported by the nonlinear evolutionary PDEs. The equations we considered, (9), (11), and (12), are obtained from the dynamical system (2) describing one-dimensional chain of prestressed elastic bodies. Equation (8) obtained in [5] from this model without resorting to the method of multi-scaled decomposition possesses the compacton solutions which fail the stability test. Numerical simulations show that the compacton solutions supported by equation (8) are destroyed in a very short time.

In contrast with the above, equations (9) (resp. (11)), which are obtained using formal multiscale decomposition, possess families of bright (resp. dark) compacton solutions which appear to be stable. This is backed both by the stability test and the results of the numerical simulations.

As we have shown in Sections 3–5, for generic values of the parameter $n$ equation (9) does not possess an infinite set of higher symmetries or other signs of complete integrability such as infinite hierarchies of conservation laws. Nevertheless the compacton solutions to this equation possess some features which are characteristic for “genuine” soliton solutions. In this connection it would be interesting to compare the traveling wave solutions for the distinguished case $n=-2$ with any other equation of the family (9) with negative $n$. Qualitative analysis of the factorized equations describing the TW solutions shows that there are no compacton solutions for the models with the negative $n$, but nevertheless all of them seem to possess periodic solutions resembling peakons. It would be interesting to find out whether there is any difference in the qualitative behavior of periodic solutions of the only integrable case ($n=-2$) in comparison with the periodic TW solutions supported by the model characterized by other values $n<0$. Perhaps the differences will be manifested in the stability properties as this is the case with the soliton solutions supported by the family of the KdV-type equations.

A characteristic feature of equations (9), (11) related to the decomposition we used is that they describe processes with “long” temporal and “short” spatial scales. Hence it is rather questionable whether these equations can adequately describe a localized pulse propagation in discrete media in the situation when the distance between the adjacent particles is comparable to the compacton width $\Delta x$. In fact, making the “reverse” transformations $X\rightarrow\xi\rightarrow x$ we get the following formula for the width of the compacton solution (33) in the initial coordinate system:

[TABLE]

this is nothing but equation (1.130) from [5]. For $n=3/2$, corresponding to the Hertzian force between spherical particles, we get $\Delta x\approx 4.96a$. It is then interesting to notice that the same results for the particles with the spherical geometry were obtained in the course of numerical simulations, and experimental studies [1, 3, 2, 54, 55]. We wish to stress that results of our analysis as well as the main conclusions are in agreement with the earlier publications by other authors. In particular, P. Rosenau notes, when considering the general models of dense chains [19], that the natural separation of scales leading to an unidirectional PDE of first order in time does not exist.

Acknowledgments

VV gratefully acknowledges support from the Polish Ministry of Science and Higher Education. The research of AS was supported in part by the RVO funding for IČ47813059, and by the Grant Agency of the Czech Republic (GA ČR) under grant P201/12/G028. AS gratefully acknowledges warm hospitality extended to him in the course of his visits to AGH in Kraków.

We are pleased to thank the anonymous referee and M.V. Pavlov for useful suggestions.