The fifth-order partial differential equation for the description of the {\alpha}+{\beta} Fermi-Pasta-Ulam model

TL;DR

This paper derives and analyzes a new fifth-order nonlinear PDE for the Fermi-Pasta-Ulam model, extending the Gardner equation, and investigates its solutions and integrability properties.

Contribution

It introduces a higher-order PDE for the Fermi-Pasta-Ulam model, incorporating quadratic and cubic interactions, and analyzes its integrability and solutions.

Findings

The equation does not pass the Painlevé test.

It cannot be solved by inverse scattering transform.

Traveling wave solutions are obtained using analytical methods.

Abstract

We study a new nonlinear partial differential equation of the fifth order for the description of perturbations in the Fermi-Pasta-Ulam mass chain. This fifth-order equation is an expansion of the Gardner equation for the description of the Fermi--Pasta--Ulam model. We use the potential of interaction between neighbouring masses with both quadratic and cubic terms. The equation is derived using the continuous limit. Unlike the previous works, we take into account higher order terms in the Taylor series expansions. We investigate the equation using the Painlev\'e approach. We show that the equation does not pass the Painlev\'e test and can not be integrated by the inverse scattering transform. We use the logistic function method and the Laurent expansion method to find travelling wave solutions of the fifth--order equation. We use the pseudospectral method for the numerical simulation of wave processes, described by the equation.

Tables

Table 1. Table 1: Periods of maximal similarity for the numerical solutions at different time moments.

t1	t2	T	diff
5	27.25	12.25	7.2705
10	32.25	12.25	8.6293
15	37.25	12.25	6.3845
20	42	12	6.4039
25	47.25	12.25	6.1192
27.25	49.5	12.25	7.1418
17.5	39.75	12.25	6.1301
19.5	41.5	12	6.6140
21	43	12	6.5962

Equations

$$\begin{array}{c}m\frac{d^2{y}_i}{d{\tau }^2}=\gamma (y_{i+1}-y_i)+\alpha (y_{i+1}-y_i{)}^2-\beta (y_{i+1}-y_i{)}^3-\\ -\gamma (y_i-y_{i-1})-\alpha (y_i-y_{i-1}{)}^2+\beta (y_i-y_{i-1}{)}^3=\\ =(y_{i+1}-2y_i+y_{i-1})(\gamma +\alpha (y_{i+1}-y_{i-1})-\\ -\beta (y_{i+1}^2+y_i^2+y_{i-1}^2-y_i{y}_{i+1}-y_i{y}_{i-1}-y_{i+1}{y}_{i-1})),\end{array}$$

$$\begin{array}{c}{y}_{i\pm 1}=y\pm h{y}_{\xi }+\frac{h^2}{2}{y}_{\xi \xi }\pm \frac{h^3}{6}{y}_{\xi \xi \xi }+\frac{h^4}{24}{y}_{\xi \xi \xi \xi }\pm \\ \pm \frac{h^5}{120}{y}_{\xi \xi \xi \xi \xi }+\frac{h^6}{720}{y}_{\xi \xi \xi \xi \xi \xi }+\dots ,\end{array}$$

$$\begin{array}{c}m\frac{d^{2}y}{d{\tau }^2}=h^2\gamma y_{\xi \xi }+2\alpha h^3{y}_{\xi }{y}_{\xi \xi }-3\beta h^4{y}_{\xi }^2{y}_{\xi \xi }+\frac{\gamma h^4}{12}{y}_{\xi \xi \xi \xi }++\frac{\alpha h^5}{3}{y}_{\xi \xi }{y}_{\xi \xi \xi }+\\ +\frac{\alpha h^5}{6}{y}_{\xi }{y}_{\xi \xi \xi \xi }-\beta h^6{y}_{\xi }{y}_{\xi \xi }{y}_{\xi \xi \xi }-\frac{\beta h^6}{4}{y}_{\xi \xi }^3-\frac{\beta h^6}{4}{y}_{\xi }^2{y}_{\xi \xi \xi \xi }+\frac{\gamma h^6}{360}{y}_{\xi \xi \xi \xi \xi \xi }.\end{array}$$

$$c^2=\frac{\gamma h^2}{m},$$

$$\begin{array}{c}y=\frac{\gamma \sqrt{\gamma }}{2c\alpha \sqrt{m}}{y}^{\prime },T=\frac{c}{2}\tau ,\end{array}$$

$$\begin{array}{c}\frac{1}{4}\frac{d^{2}y}{d{T}^2}=y_{\xi \xi }+y_{\xi }{y}_{\xi \xi }-\mu y_{\xi }^2{y}_{\xi \xi }+\delta y_{\xi \xi \xi \xi }+2\delta y_{\xi \xi }{y}_{\xi \xi \xi }+\delta y_{\xi }{y}_{\xi \xi \xi \xi }-\\ -4\mu \delta y_{\xi }{y}_{\xi \xi }{y}_{\xi \xi \xi }-\mu \delta y_{\xi \xi }^3-\mu \delta y_{\xi }^2{y}_{\xi \xi \xi \xi }+\frac{2}{5}{\delta }^2{y}_{\xi \xi \xi \xi \xi \xi },\end{array}$$

$$y(\xi ,T)=f(x,t)+\epsilon y_1(\xi ,T),x=\xi -2T,t=\epsilon T.$$

$$\begin{array}{c}{f}_{xt}+f_x{f}_{xx}-\mu f_x^2{f}_{xx}+\delta f_{xxxx}+2\delta f_{xx}{f}_{xxx}+\delta f_x{f}_{xxxx}-\\ -4\mu \delta f_x{f}_{xx}{f}_{xxx}-\mu \delta f_{xx}^3-\mu \delta f_x^2{f}_{xxxx}+\frac{2}{5}{\delta }^2{f}_{xxxxxx}=0.\end{array}$$

$$u(x,t)=f_x(x,t),$$

$$\begin{array}{c}{u}_t+u{u}_x-\mu u^2{u}_x+{\delta }^2{u}_{xxx}+2{\delta }^2{u}_x{u}_{xx}+{\delta }^{2}u{u}_{xxx}-\\ -4{\delta }^2\mu u{u}_x{u}_{xx}-{\delta }^2\mu u_x^3-{\delta }^2\mu u^2{u}_{xxx}+\frac{2}{5}{\delta }^4{u}_{xxxxx}=0,\end{array}$$

$$u_t+u{u}_x-\mu u^2{u}_x+{\delta }^2{u}_{xxx}=0.$$

$$u(x,t)=v(z),z=x-C_{0}t.$$

$$\begin{array}{c}-C_0{v}^{\prime }+v{v}^{\prime }-\mu v^2{v}^{\prime }+{\delta }^2{v}^{\prime \prime \prime }+2{\delta }^2{v}^{\prime }{v}^{\prime \prime }+{\delta }^{2}v{v}^{\prime \prime \prime }-\\ -4\mu {\delta }^{2}v{v}^{\prime }{v}^{\prime \prime }-\mu {\delta }^2{v}^{\prime 3}-\mu {\delta }^2{v}^2{v}^{\prime \prime \prime }+\frac{2}{5}{\delta }^4{v}^{(V)}=0.\end{array}$$

$$\begin{array}{c}{C}_1-C_{0}v+\frac{1}{2}{v}^2-\mu \frac{1}{3}{v}^3+{\delta }^2{v}^{\prime \prime }+\frac{1}{2}{\delta }^2{v}^{\prime 2}+{\delta }^{2}v{v}^{\prime \prime }-\\ -\mu {\delta }^{2}v{v}^{\prime 2}-\mu {\delta }^2{v}^2{v}^{\prime \prime }+\frac{2}{5}{\delta }^4{v}^{(IV)}=0.\end{array}$$

$$\begin{array}{c}{C}_2+C_{1}v-\frac{C_0}{2}{v}^2+\frac{1}{6}{v}^3-\mu \frac{1}{12}{v}^4+\frac{{\delta }^2}{2}{v}^{\prime 2}+\frac{{\delta }^2}{2}v{v}^{\prime 2}-\\ -\frac{1}{2}\mu {\delta }^2{v}^2{v}^{\prime 2}+\frac{2}{5}{\delta }^4{v}^{\prime }{v}^{\prime \prime \prime }-\frac{1}{5}{\delta }^4{v}^{\prime \prime 2}=0.\end{array}$$

$$-\frac{1}{2}\mu {\delta }^2{v}^2{v}^{\prime 2}+\frac{2}{5}{\delta }^4{v}^{\prime }{v}^{\prime \prime \prime }-\frac{1}{5}{\delta }^4{v}^{\prime \prime 2}=0.$$

$$v\simeq \pm \frac{4\sqrt{5}\delta }{5\sqrt{\mu }}\frac{1}{z}.$$

$$v=\pm \frac{4\sqrt{5}\delta }{5\sqrt{\mu }}\frac{1}{z}+v_j{z}^{j-p},$$

$$j_1=-1,j_{2,3}=\frac{5\pm i\sqrt{7}}{2}.$$

$$v=A_0+A_1\theta ,\theta =\frac{1}{1-\exp (-k(z-z_0))},$$

$$\begin{array}{c}{C}_0=\frac{15-56\mu }{180\mu },C_1=-\frac{15+56\mu }{360{\mu }^2},\\ C_2=-\frac{3136{\mu }^2+1680\mu +225}{43200{\mu }^3},\end{array}$$

$$\begin{array}{c}v(z)=\pm \sqrt{\frac{2\left(28\mu +15\right)}{15{\mu }^2}}\left(\frac{1}{1-\exp (-k(z-z_0))}-\frac{1}{2}\right)+\frac{1}{2\mu },\\ k=\frac{\sqrt{6\mu \left(28\mu +15\right)}}{12\mu \delta }.\end{array}$$

$$\begin{array}{c}v(z)=\pm \sqrt{\frac{2\left(28\mu +15\right)}{15{\mu }^2}}\left(\frac{1}{1+\exp (-k(z-z_0))}-\frac{1}{2}\right)+\frac{1}{2\mu },\\ k=\frac{\sqrt{6\mu \left(28\mu +15\right)}}{12\mu \delta }.\end{array}$$

$$\frac{1}{1+\exp (-k(z-z_0))}=\frac{1}{2}\left(1+\tanh \left(\frac{k(z-z_0)}{2}\right)\right),$$

$$\begin{array}{c}v(z)=\pm \sqrt{\frac{\left(28\mu +15\right)}{30{\mu }^2}}\tanh \left(\frac{k(z-z_0)}{2}\right)+\frac{1}{2\mu },\\ k=\frac{\sqrt{6\mu \left(28\mu +15\right)}}{12\mu \delta }.\end{array}$$

$$v(z)=H+\frac{A+B{\wp }_z(z-z_0,g_2,g_3)}{C+\wp (z-z_0,g_2,g_3)},$$

$$\begin{array}{c}H=\frac{1}{2\mu }A=0,B=\frac{2\delta }{\sqrt{5\mu }},C=-\frac{28\mu +15}{288\mu {\delta }^2},\\ C_0=-\frac{2985984{\delta }^6{g}_3{\mu }^3+78400{\mu }^3+50400{\mu }^2+10800\mu +3375}{6480{\mu }^2(15+28\mu )},\\ C_1=-\frac{2985984{\delta }^6{g}_3{\mu }^3+78400{\mu }^3+80640{\mu }^2+27000\mu +3375}{12960{\mu }^3(15+28\mu )},\\ C_2=-\frac{336{\delta }^6{g}_3}{25\mu }+\frac{1552{\delta }^6{g}_3}{25(15+28\mu )}-\frac{343}{150(15+28\mu )}-\frac{833}{180\mu (15+28\mu )}-\\ -\frac{7}{2{\mu }^2(15+28\mu )}-\frac{75}{64{\mu }^3(15+28\mu )}+\frac{1296{\delta }^6{g}_3}{5\mu (15+28\mu )}-\\ -\frac{75}{512{\mu }^4(15+28\mu )}+\frac{3577}{48600\mu }+\frac{721}{4320{\mu }^2}+\frac{103}{1152{\mu }^3}+\frac{65}{4608{\mu }^4},\\ g_2=-\frac{5971968{\delta }^6{g}_3{\mu }^3-21952{\mu }^3-35280{\mu }^2-18900\mu -3375}{20736{\delta }^4{\mu }^2(15+28\mu )}.\end{array}$$

$$y^{\prime 2}=4y^3-g_{2}y-g_3.$$

$$y^{\prime }=\pm \sqrt{4\left(y-g_3^{1/3}\right){\left(y+\frac{1}{2}{g}_3^{1/3}\right)}^2}.$$

$$z-z_0=\pm \int \frac{dy}{\sqrt{4\left(y-g_3^{1/3}\right){\left(y+\frac{1}{2}{g}_3^{1/3}\right)}^2}}.$$

Full text

Keywords: Fermi–Pasta–Ulam model; continuous limit; Painlevé test; exact solution; pseudospectral method

1 Introduction

The Fermi–Pasta–Ulam (FPU) model was first studied in work [1]. It describes perturbations in the chain of nonlinear coupling among masses. It is shown [1] that energy remains in a very few modes and the long-time dynamic of the system is recurrent. This fact was called the FPU paradox. A set of approaches was used to explain it (see for example [2, 3, 4, 5, 6, 7]). The main ideas of these approaches and their results can be found in work [8].

Usually $\alpha$ or $\beta$ FPU model is considered in literature, but in paper [9] the $\alpha+\beta$ FPU model is investigated. In this paper, we also use the $\alpha+\beta$ model but we use the continuous limit approximation for the investigation.

This approach was used for the first time by N. Zabusky and M. Kruskal in work [10], where the authors derived the Korteweg–de Vries (KdV) equation for the description of the $\alpha$ FPU model. This equation is integrable by the inverse scattering transform and has been widely studied by now ([11, 12]). Zabusky and Kruskal explained the FPU paradox (recurrence of initial conditions) using the numerical simulation of the KdV equation. Another result of work [10] was an introduction of soliton. In article [13] the author took into account high order (in comparison with [10]) terms in the Taylor series expansions for the continuous limit approximation of the $\alpha$ FPU model. The generalized KdV equation was obtained for a more accurate description of wave processes in the FPU model. Some exact solutions of the derived equation were found . It is shown [13] that if one takes into account more terms in the Taylor series and acts in analogy to paper [10] then he can not derive an integrable equation for the description of wave processes in the FPU mass chain. Recurrence of initial statement is also unrepresentative for wave processes, described by the fifth-order equation with arbitrary initial conditions. In paper [14] wave processes for the generalized KdV equation are simulated numerically.

This brings up the question, whether the situation is similar for the $\beta$ and the $\alpha+\beta$ FPU models? It is known [11] that one can obtain the Gardner equation for the description of the $\alpha+\beta$ FPU model using the continuous limit. The Gardner equation is integrable by the inverse scatterring transform and can be reduced to the modified KdV equation [15]. If we take more terms in the Taylor series expansion, we can derive the fifth-order partial differential equation for the description of the $\alpha+\beta$ FPU model. To the best of our knowledge, this equation has not been obtained and investigated before. This is the main aim of the current manuscript.

The rest of the work is organized as follows. In section 2 we derive the fifth-order partial differential equation for the description of the $\alpha+\beta$ FPU model. In section 3 we use the Painlevé test to analyze the derived equation. The logistic function method is used in section 4 to find exact solutions of the derived equation. In section 5 we obtain an elliptic solution for this equation. In section 6 we present the results for the numerical simulation of wave processes, described by the derived equation. In section 7 we briefly discuss our results.

2 The fifth-order partial differential equation for the description of the $\alpha+\beta$ Fermi–Pasta–Ulam model

Let us consider the $\alpha+\beta$ Fermi–Pasta–Ulam model. It appears as the following system of equations:

[TABLE]

where $y_{i}$ is a displacement of the particle number $i$ from its equilibrium position, $i$ runs through values from $1$ to $N$, where $N$ is the number of particles, $m$ is the mass of a single particle, $\alpha,\;\beta,\;\gamma$ are some positive constants that characterize the potential of interaction and $\tau$ is the time. Let us find the continuous limit approximation of the model. In order to do this we use the following Taylor series expansions:

[TABLE]

where $y=y(\xi,\tau)$ is a displacement of the infinitely small part of the chord with coordinate $\xi$ in time moment $\tau$, $h$ is a small parameter. Substituting expansion (2) into the system of equations (1) we derive a partial differential equation for the description of the system dynamics in the continuous limit approximation. Let us take into account only terms up to $h^{6}$.

[TABLE]

Using a new parameter:

[TABLE]

and scaling transformations:

[TABLE]

we reduce equation (3) to the form:

[TABLE]

where $\mu=\frac{3\beta\gamma}{4\alpha^{2}},\;\delta=\frac{mc^{2}}{12\gamma}$ and primes are omitted. According to [11], the solution of this equation can be found as two oppositely directed running waves. Following [10], we only study the wave propagating to the right. For this purpose we search for the solution of equation (6) in the following form:

[TABLE]

Here $f(x,t)$ corresponds to the wave profile on long distance. We omit members that have high order in epsilon. Then from equation (6), we derive an equation for $f(x,t)$:

[TABLE]

Let us use a new variable:

[TABLE]

then equation (8) takes the form:

[TABLE]

Equation (10) is the logical extension of the equation from the article [13], because in case of quadratic potential of interaction ($\mu=0$), one can see that equation (10) takes the form of equation, studied in [13]. Parameter $\delta$ has also been introduced in the same form as in earlier studies of N. Zabusky [10] and N. Kudryashov [13]. Another type of the extended fifth-ordered KdV equations is considered in [16]. Cases of the integrability for a lot of fifth-order PDEs can be found in [17]. However, equation (10) is new and worth being investigated.

If we use only terms up to the $h^{4}$ in (3), we obtain the Gardner equation for the description of wave processes in the $\alpha+\beta$ FPU model:

[TABLE]

As it was mentioned above, equation (11) can be reduced to the modified KdV equation. It is integrable and has solutions in the form of breather waves and solitons.

3 Painlevé test for equation (10)

Let us use the Painlevé test for equation (10) using the traveling wave variables. Below we will show that it is not necessary to analyse equation (10) in general case. The traveling wave variables for equation (10) have the form:

[TABLE]

Using (12) we get from (10):

[TABLE]

Integrating equation (13) with respect to $z$ we find that the first integral of equation (13) has the form:

[TABLE]

We multiply both sides of equation (14) by $v^{\prime}$ and integrate the result one more time with respect to z. Then we derive the second integral for equation (13):

[TABLE]

Hence we have reduced the fifth–order differential equation to the third-order differential equation. It simplifies the process of further investigation.

We apply the Painlevé test for equation (15) using three steps (see for example [18]). On the first step we have to determine the order of the pole and the first term in expansion of the solution in the Laurent series. We use the equation with the leading terms of equation (15):

[TABLE]

Substituting $v=\frac{a_{0}}{z^{p}}$ into equation (16) we find that the general solution of equation (15) has the pole of the first order ($p=1$) and two branches:

[TABLE]

On the second step we determine Fuchs indices. When an equation passes the Painlevé test, Fuchs indices indicate the numbers of arbitrary constants in the Laurent expansion of the general solution. To find the Fuchs indices, we substitute expression:

[TABLE]

into equation (16). Equating the linear term with respect to $v_{j}$ in the expression derived to zero, we get Fuchs indices for equation (15):

[TABLE]

Fuchs indices are the same on both branches of the solution. Two of three indices are complex numbers. It means that equation (15) does not pass the Painlevé test. Thus it does not have a general solution. However, it is possible to find some exact solutions of equation (15). This fact is demonstrated below.

The third step of the test includes examination, whether constants in the Laurent expansion of the solution can be taken as an arbitrary according to the numbers of the Fuchs indices. In our case, this step is not necessary. We also do not need to analyze equation (10) using the Painlevé approach because it is known that if a travelling wave reduction of an equation does not pass the Painlevé test, then the equation does not pass it in general case.

We see the analogy to the equations describing the $\alpha$ FPU model. If the high order terms in the Taylor series expansions are not taken into account, then it is possible to describe the model with an integrable equation. This equation admits recurrence of an initial statement and excludes any chaotic behavior. However, if higher order terms are taken into account then it is impossible to derive an integrable equation.

4 The logistic function method for finding the exact solution of equation (15).

There are various methods to find a solution of a partial differential equation (see for example [19, 20, 21, 22, 23]). We use the logistic function method. This method often allows to use the results of the Painlevé analysis for a partial differential equation for the construction of its exact solutions [24, 25, 26]. Let us search for the solution of equation (15) in the following form:

[TABLE]

where $\theta(z)$ is a so-called logistic function. During the Painlevé test, we found the pole order of the general solution. This result is taken into account in expression (20), where the pole order corresponds to the order found above. We expand (20) in the Laurent series, substitute it into equation (15), equate coefficients at the same powers of $z$ to zero, and solve the derived system of algebraic equations. As a result, we find constants:

[TABLE]

and exact solution of equation (10):

[TABLE]

The velocity of the wave propagation $C_{0}$ can be both positive and negative It depends on the correlation of coefficients in the potential of interaction between neighbouring particles. Thus the wave can propagate either to the right or to the left. The velocity as the function of the parameter $\mu$ is presented on figure 1.

Choosing an arbitrary constant $z_{0}=z_{0}+i\pi/k$ in solution (22), we obtain the solution without poles on the real line:

[TABLE]

Using the following relation:

[TABLE]

we get the solution of equation (15) in the form:

[TABLE]

Solution (25) has the form of the kink and is illustrated on figure 2 at $t=0$ for different values of parameters. It can be seen from the picture, that the parameter $\delta$ corresponds to the slope of the wave front, while the parameter $\mu$ is inverse to the amplitude of the perturbation and shifts the solution alongside the ordinate axis.

5 Elliptic solution for equation (15)

Let us use the Laurent expansion method [13, 14] to obtain the elliptic solution for equation (15). The Laurent expansion method uses the results of the Painlevé test as well as the logistic function method. To correlate the pole’s order with the order found above, we look for the solution of equation (15) in the following form [27]:

[TABLE]

where $H,\;A,\;B,\;C$ are some constants, $\wp$ is the elliptic Weierstrass $\wp$-function with invariants $g_{2}$ and $g_{3}$. We expand expression (26) in the Laurent series, substitute this expansion for $v$ in equation (15), equate coefficients at the same powers of $z$ to zero, and solve the derived system of algebraic equations. As a result we find the following relations:

[TABLE]

Note that $\delta,\mu$, and $g_{3}$ are parameters in the solution. One can consider solution (26) with constants (27) and exclude poles from the real line. This solution is presented on figure 3 at $\delta=1,\;\mu=0.5,\;C_{0}=5$.

It is useful to investigate the case of the degeneration of the elliptic solution. We take an equation for the $\wp$ function:

[TABLE]

This equation has the solution $y=\wp(z-z_{0},g_{2},g_{3})$. We explore the right side of equation (28) as an algebraic equation. As $g_{2}=3g_{3}^{2/3}$ it has roots $y_{1}=g_{3}^{1/3},\;y_{2,3}=\frac{1}{2}g_{3}^{1/3}$. In this case, equation (28) takes the form:

[TABLE]

We split the variables and obtain the solution for equation (29):

[TABLE]

Considering $y$ as the function of $z$ we derive the solution of equation (28) in quadratures:

[TABLE]

Hence we have the connection between the Weierstrass function in the case of degradation and trigonometric functions.

[TABLE]

It means, that in case of $g_{2}=3g_{3}^{2/3}$, an elliptic solution of the differential equation can be expressed in terms of the trigonometric functions. The period of the elliptic solution converges to infinity. If the poles can be excluded from the real line, then this approach allows to obtain solitary wave solutions from the periodic ones (see for example [11]). However for equation (10) the poles stay on the real line. Thus, solution (26) with coefficients (27) can not be reduced to the form of the solitary wave solution.

6 Numerical simulation of wave processes, described by equation (10) with periodic boundary conditions

We use the pseudospectral method [28, 29] for the numerical simulation of the boundary value problem for equation (10) with periodic boundary conditions. This method was successfully used for the numerical simulation of different problems of mathematical physics [26]. We chose the pseudospectral method because it allows us to solve boundary problems for the equations of high order. This approach is not as demanding to the time step as methods, based on finite difference schemes.

The main idea of the algorithm is the following. We distinguish the linear($L$) and the nonlinear ($N$) operators for equation (10):

[TABLE]

and introduce the spatial-temporal grid:

[TABLE]

where $h$ and $\tau$ are spatial and time steps correspondingly, and $N_{x},\;n$ are integer numbers. In each grid point we apply the Fourier space transform:

[TABLE]

where $\hat{u}$ is the Fourier transform of the function $u$. The derived differential equation is solved with one of the appropriate methods. We use the fourth-order Runge–Kutta method.

The smoothing error appears while using the discrete Fourier transform. To correct this the filter is used. The Fourier transform is calculated according to the rule: $\widehat{\frac{\partial^{n}f}{\partial{x^{n}}}}=(ik_{x})^{n}p\hat{f}$, where $f$ is the transformed function and $p$ is the filter. There are various filters, we use the one that is constructed by the rule of $\frac{2}{3}$:

[TABLE]

After that, the original function $u$ is calculated with the inverse Fourier transform.

The Fast Fourier transform is used to accelerate the calculation. It requires spatial grid with the number of grid points equal to the power of 2. Another method to reduce the computational time is the integrating factor method. According to this method, we substitute the unknown function in (35) in the following way:

[TABLE]

Final formulae for the solution of ODE with the Runge-Kutta method are:

[TABLE]

To check our numerical strategy, we use exact solution (25) as initial conditions for the simulation. To keep periodic boundary conditions, we add the reflective part to the exact solution. Hence only a part of the initial conditions is an exact solution for equation (10). The result is presented on figure 4 at $\delta=0.6$ and $\mu=2$. The solid line is for the numerical solution and the dotted line is for the exact one. It can be seen from the picture, that the solution propagates to the right and does not change its shape. Part of the numerical solution that is obtained from the exact solution coincides with the exact solution. Part of the numerical solution that does not appear to be an exact solution travels without any changes in its shape and velocity or an emission. The computational error was calculated as a ratio $err=\frac{\max|u_{analytic}-u_{calc}|}{\max|u_{calc}|}$, where $u_{analytic}$ and $u_{calc}$ are the exact and the numerical solutions of equation (10) correspondingly. To calculate $err$ we use only grid points, where the initial conditions for the numerical solution coincides with the exact solution. Evolution of the computational error with time is presented on figure 5. We see from the figure that the error is below $6\times 10^{-3}$, the accuracy of the simulation is quite acceptable.

At $\mu=0$, equation (10) describes waves in the $\alpha$ FPU model. Solitary wave solutions that look like solitons are found for this case in [13]. We know [30] that soliton-like waves described by non-integrable equations usually propagate with radiation. At $\mu\neq 0$, equation (10) does not have such exact solutions. We can investigate numerically the influence of the additional terms at $\mu\neq 0$ on the stability of such solitary wave solution. We use the solution of equation (10) at $\mu=0$ from [13]:

[TABLE]

as initial conditions. It is reasonably to perform the simulation with small values of the parameter $\mu$. In this case the influence of the quadratic members in the potential of interaction on wave propagation is significantly bigger than the influence of the cubic terms. Figure 6 illustrates the comparison between numerical solutions of boundary value problems for equation (10) at $\mu=0$ (the dotted line) and $\mu=0.05$ (the solid line) with $k=1,\;\delta=2$. It can be seen from the reference time that solitary wave slows down and the emission appears at $\mu\neq 0$. Waves created by the emission grow at $t=9$, and the hole destruction of the solitary wave is being observed at $t=17$. At the same time, the numerical solution for boundary value problem at $\mu=0$ changes neither its speed nor its shape.

During the simulation one can notice the similar wave profiles at different time moments. The example is presented on figure 7. We can assume that the recurrence of the perturbations appears as it was observed by Zabusky and Kruskal in work [10]. To study this, we use the following approach.

We fix the numerical solution for equation (10) in an arbitrary time moment. 2. 2.

We skip a short period of time. Then on each time layer we use the circular convolution method to compare the fixed and current solutions. 3. 3.

We plot the result of calculations as evolution of the absolute value of difference. The difference is calculated in the same way as the computational error.

The result of using this approach for $t=5$ is presented on figure 8. From this figure we can see that minimum of norm of the difference between the two numerical solutions has an obvious period. We can use the program to find this period. The result is presented in table 1 with an accuracy up to $0.25$. The time moment when the numerical solution was fixed is in the first column of the table. Time moment for the minimum of norm of the difference is in the second column. The third column presents the difference between the first and the second columns. Absolute value of the difference between the solutions at these time moments is in the fourth column. We see from the table that the period of the recurrence is approximately equal to 22.25. Thus that there is no absolute recurrence of the perturbations but we can observe quasi recurrent behaviour. One of the possible reasons for this quasi recurrent dynamics is the short simulation time.

As we have shown above, equation (10) extends Gardner equation (11). The Gardner equation has solutions in the form of the soliton:

[TABLE]

Let us investigate the wave processes, described by equation (10) with the solitary wave (40) as initial conditions. The result is presented on figure 9 at $\delta=1,\;\mu=0.1,\;C_{0}=1$. From the picture we can see that at $t=8.75$ the wave profile is being deformed, creating the burst of waves. This process is accompanied with additional perturbations. The waves are slowing down. Situation is developing at $t=16.25$ and at $t=50$ we can see almost the same picture.

In paper [10] the authors provided numerical modelling of wave processes in the FPU mass chain, described by the Korteweg-de Vries equation. They have conducted the modelling at $\delta=0.022$ in the spatial area $0\leq x\leq 2$. The following initial conditions were used: $u(x,0)=\cos(\pi x)$. We have compared wave processes, described by equation (10) at $\mu=1$ with the results for the KdV equation using the same initial conditions as in [10]. The result is illustrated on figure 10. The dotted line is for the KdV equation and the solid line is for equation (10). We can see from the figure that the numerical solutions of the boundary value problem for the KdV equation consists of solitons. The solitons interact without emission. The solution at $t=10.6$ has the same shape as at $t=1.14$. However, the numerical solution of the boundary value problem for equation (10) demonstrates a chaotic behaviour. Initial perturbations do not divide into the separated solitons with the elastic interaction. As a result the emission appears and energy leaves initial modes. The recurrence of initial statement can not be observed for wave processes described by equation (10).

7 Conclusion

We have derived the fifth-order partial differential equation for the description of wave processes in the $\alpha+\beta$ Fermi–Pasta–Ulam model, taking into account high order terms in the Taylor series expansions. We have shown that the equation does not pass the Painlevé test but it has some exact solutions. We have used the logistic function method to find a kink-type solution for the travelling wave reduction of the derived equation. The elliptic solution for the derived equation is constructed with the Laurent expansion method. We have numerically investigated wave processes described by the equation. We have shown that the recurrence of the initial conditions is not specific for an arbitrary initial state. Although we have found a quasi recurrent behaviour of the numerical solution.

8 Acknowledgments

This research was supported by Russian Science Foundation Grant no. 14-11-00258.