Exact Solutions and Numerical Simulation of the Discrete Sawada–Kotera Equation

We investigated an integrable five-point differential-difference equation called the discrete Sawada–Kotera equation. On the basis of the geometric series method, a new exact soliton-like solution of the equation is obtained that propagates with positive or negative phase velocity. In terms of the Jacobi elliptic function, a class of new exact periodic solutions is constructed, in particular stationary ones. Using an exponential generating function for Catalan numbers, Cauchy’s problem with the initial condition in the form of a step is solved. As a result of numerical simulation, the elasticity of the interaction of exact localized solutions is established.


Introduction
The theory of integrable differential-difference equations (DDE), based on the study of symmetries, conservation laws and recursion operators, is actively developing at present. As in the continuous case, the integrability of an equation is understood in the sense of the existence of Lax pairs or the infinite hierarchy of symmetries [1]. The theoretical results obtained in this direction are of great practical importance, since "models of a wide range of physical processes are inherently discrete" [2,3]. Most often, nonlinear models of physical processes allow analytical research only on the basis of the continuum limit of the initial discrete equation [4]. This is due to the fact that the construction of exact solutions of nonlinear DDEs other than the Toda, Volterra and Ablowitz-Ladik lattices is difficult [5].
Most integrable lattices have the Korteweg-de Vries and nonlinear Schrodinger equations as a continuum limit [6,7]. Recently, as a result of the classification of five-point DDEs [8], integrable discretizations of the Sawada-Kotera and Kaup-Kupershmidt equations have been identified.
The aim of this work is to conduct a qualitative study and numerical simulation of the integrable five-point discrete Sawada-Kotera equation The considered DDE (1) is a linear combination of two famous integrable equations, the Volterra chain and the Ito-Narita-Bogoyavlensky chain [9][10][11], each of which has the Korteweg-de Vries equation as a continuous limit. It is known [12] that as these chains belong to different hierarchies, the flows corresponding to them do not commute, so the integrability of their linear combination is a nontrivial fact, which is established by constructing the corresponding Lax pair. It was shown in [12] that the introduction of new variables, allows one to obtain the integrable Sawada-Kotera equation [13] U τ = U xxxxx + 5UU xxx + 5U x U xx + 5U 2 U x from Equation (1) in the continuous limit for ε → 0. Therefore, Equation (1) is called the discrete Sawada-Kotera equation [14,15].
In [14], the general structure of the N-soliton solution of Equation (1) is given in terms of the Pfaffians. However, in specific cases of modeling processes of a discrete nature, information on the explicit form of physically realizable exact solutions is especially important for a researcher. Therefore, to achieve this goal, we solve several independent problems: The construction of exact solitary-wave and periodic solutions, the initial problem, and a numerical simulation of the obtained solutions.
The rest of the article is structured as follows. In the next two sections, an exact solitary-wave solution and a periodic solution will be constructed. In the fourth section, some combinatorial dependencies will be established when solving an initial value problem for Equation (1). In the fifth section, exact stationary solutions are discussed, and in the last section, numerical simulations are performed that demonstrate the stability of exact localized solutions and the elasticity of their interaction.

Solitary Wave Solution
Let us find an exact solitary-wave solution to the Equation (1). In [16,17], algorithms were developed for computer mathematics systems and exact kink-like solutions were obtained for Toda, Volterra, and Ablowitz-Ladik lattices and their generalizations. To construct an exact localized solution, we will use a modification of the geometric series method [18]. In the continuous case, this approach allows one to obtain exact soliton-like solutions of integrable and non-integrable equations with a polynomial dispersion relation as the sum of the geometric series of the perturbation method in powers of exponential function.
The transition to the traveling wave variable z = dn + ωt allows us to write Equation (1) as follows After substituting in Equation (2) a functional series with unknown coefficients C k where m = ±1, ±2, δ = e d , we group the terms in powers of the exponential function. In the lower order, at e z , we obtain a "linear dispersion relation" for determining the frequency ω: Successively equating to zero the coefficients at e 2z , e 3z , ..., we calculate C 2 , C 3 , ...: , ...
All coefficients, starting from C 2 , are expressed through three arbitrary parameters C 0 , C 1 and δ. Looking at Equations (6), it is easy to write the general formula for the coefficient C n ; however, this is not necessary. The series (3) with coefficients (6) is a geometric series and its sum gives an exact solution to Equation (2). To verify this, we use the remarkable property of Padé approximants: Successive diagonal approximants [1/1], [2/2], [3/3],. . . , calculated for an arbitrary geometric power series, coincide, starting from some order [N/N].
Let us prove the presence of this property. As is known, the Padé approximant of order [M/N] for power series S(x) = ∑ k a k x k is the ratio of polynomials whose Maclaurin series expansion T M,N (x) coincides with S(x) as long as possible [19,20]. To calculate the coefficients p k , q k of the ratio (7) we use the first M + N + 2 terms of the series S(x). Formal representation, means that the Maclaurin series T M,N (x) for Equation (7) is a geometric series with the first term P M and the denominator (1 − Q N ), and the approximant (7) itself is the exact sum of its Maclaurin series.
Suppose that for geometric series T M,N (x) we calculated a higher order approximant, for example, Obviously, its own Maclaurin series T M+1,N+1 (x) is also geometric and its first M + N + 4 terms coincide with the terms of T M,N (x). The geometric series can be continued in its first terms uniquely, therefore, the series T M,N (x) and T M+1,N+1 (x) coincide and after them the approximants [M/N] and [M + 1/N + 1] coincide identically, as required.
Replacing e z = Z in Equation (3), we begin to sequentially calculate for Equation (3) the diagonal Padé approximants, starting with [1/1]. It turns out that the second and third order approximants coincide: Without calculating an approximant of higher order, we check whether Equation (9) gives an exact solution to the problem. Returning to variables z, d in Equations (9) and (5) we obtain where (10) and (11) into Equation (2) turns the latter into identity, therefore, Equation (10) is an exact solution of Equation (2) under condition (11). After substitution z = dn + ωt expression (10) becomes a solitary-wave solution of the original Equation (1). The shape of this solution depends on the relationship between the parameters A and d and varies from a flat-peak pulse to a pulse with a sharp peak on the pedestal C 0 (Figure 1).

Periodic Solution
The leading terms balance gives a simple pole for solution of Equation (2). In accordance with the truncated expansion method [21], we will use the ansatz with the Jacobi elliptic function: Using the addition theorem for Jacobi elliptic functions, we have: where The constants S 2 , C 2 , D 2 are expressed through S, C, D by means of double argument formulas: Taking into account Equation (14), we substitute Equations (12) and (13) into Equation (2). After dividing by a common factor 4B cn (z, m) dn (z, m), the left-hand side of Equation (2) reduces to a fourth-order polynomial with respect to sn (z, m): where The system of nonlinear equations {A n = 0} , n = 0...4 , supplemented by the relations between Jacobi elliptic functions cn, sn and dn : can be solved by any of the modern computer mathematics software systems, for example, Maple. The solution contains several branches corresponding to stationary waves (ω = 0) and traveling waves (ω = 0). Two combinations of parameters correspond to traveling waves. The first of them has the form where all signs can be selected arbitrarily. The second combination is obtained from the first one by substitution D → −D. Note that the solution Equation (12) at each moment of time is determined only on a uniform grid nd + z 0 , n ∈ Z. If the ratio of the grid step d to the period T of function sn(z, m) is not a rational number, then the period of the grid function u n = A + B sn(nd + z 0 , m) tends to infinity and we obtain a non-periodic solution. It can be shown that function (12) under conditions (19) gives an exact periodic solution of Equation (1). This solution is real for 0 < D < 1 2 and its period T has six nodes of the lattice for any admissible value of D. With an increase in parameter D frequency ω and amplitude B decrease monotonously ( Figure 2). If we take into account that the intrinsic amplitude of function sn(z, m) when |m| > 1 is A sn = 1 m , then the total oscillation amplitude A sn B of solution (12) turns out to be numerically the same as the frequency ω. Examples of a traveling waveform are shown in Figure 3.

Exact Solution and Catalan Numbers
An interesting property of the modified Volterra lattice, is revealed in articles [22,23]. It turns out that for a special initial condition, we can write the exact solution of Equation (20) in closed form, and this solution is neither of periodic type, nor of traveling-wave type.
We apply the methodology [22,23] to the lattice Equation (1) with the same initial condition (21). First, note that u 0 (t) ≡ 0. Due to this equality, Equation (1) with any positive number n are independent of the values of the grid function at nodes with negative numbers u −1 , u −2 , .... We will solve the initial problem Equations (1), (21) for positive nodes u 1 , u 2 , ... By sequentially differentiating both sides of Equation (1) with respect to t and replacing the first derivatives in the right-hand sides with the help of Equation (1), we can express any higher derivative u (k) n (t) in terms of u p (t), for example, Substituting to Equation (22) the initial conditions (21), we obtain the sequence u (0) n (0), u For odd nodes (n = 2p + 1) we observe that for k > 0 the initial values of the derivatives u (k) 2p+1 (0) are zero, at least for the first eight calculated values. It can be assumed that during the evolution of the system, the values at these nodes do not change: In fact, substituting Equation (24) into Equation (1) for an odd knot, we have For even nodes the situation is more complicated. For n = 2 we have a sequence of Catalan numbers c n = (2n)! n!(n+1)! , for n = 4 we get a sequence whose generating function is inverse to the generating function for Catalan numbers, and the sequence for n = 6 is no longer recognized by Sloane's online sequence encyclopedia [24]. We explicitly write Equation (1) for even nodes taking into account (24): Obviously, Equations (26) coincide with (20) up to renumbering. In [22] it is shown that the exact solution of Equation (20) with the initial condition (21) is where w −2 = w −1 = 1 and, for k ≥ 0, In Equations (28) the notation u 1 = f (t) is taken. Thus, all unknown functions u 2 (t), u 3 (t), ... are expressed in closed form through function u 1 (t) and its derivatives u 1 (t), u 1 (t), .... You can find function u 1 (t) itself through the exponential generating function for the Catalan numbers c n , which serve as the coefficients of the expansion of u 1 (t) in the Maclaurin series: where I 0 (z) is the modified Bessel function, which is a solution to the Cauchy problem Unfortunately, the solutions of both the modified Volterra lattice Equation (20) and Equation (1) demonstrate an exponential increase in time and, therefore, are limited in physical applications ( Figure 4). The behavior of the solution depends on the node number n: We have u n = 1 for n = 1, 3, 5, ...; lim t→∞ u n (t) = +∞ for n = 2, 6, 10, ...; lim t→∞ u n (t) = 0 for n = 4, 8, 12, ... .
As N increases, the set of stationary solutions expands. For example, if for N = 3 we have a unique solution of the form {0, 0, a}, then for N = 4 we find four solutions {0, a, 0, b}, {1, a, 1, b}, {−1, a, −1, b} and a, b, 1 a , 1 b . Note that the first three options for N = 4 make it possible to construct non-periodic solutions, since the quantities a and b can be specified arbitrarily in each set of four numbers. The graph of the 4th option for N = 4 at a = 2, b = −3 is shown in Figure 5. Another variant of the non-periodic stationary solution follows from Equation (10), where the constants C 0 and d should be selected from the condition that the frequency ω defined by Equation (11) vanishes.

Numerical Simulation
The value of exact solutions decreases if they are unstable to small deviations. Naturally, these deviations are introduced, for example, by the finite-difference approximation of the time derivative used in the numerical construction of solutions. The invariance of the shape of the traveling-wave solution for a sufficient period of time indicates both the convergence of the applied difference scheme and the correctness of the analytical solution of the problem. Figure 6 shows the evolution of the solution of Equation (1) with the initial condition, the graph of which is shown in Figure 1a (impulse with a sharp peak); Figure 7 corresponds to the initial condition from Figure 1b As we see in Figures 6b and 7b, dependences u n (t) for all nodes are identical and coincide in form with the solution (1.10), the graphs of which are shown for some combinations of parameters in Figure 1. It can be seen from Figure 7a that over the time T the graph is completely restored to its shape, as it should be for a traveling wave solution. Figure 8 shows the interaction of two solitary waves of different amplitudes. The collision of waves occurs according to the elastic scenario, with a complete restoring of the shape and amplitude of each wave and a certain phase shift-the catching-up wave shifts forward in the direction of travel, and the catch-up wave shifts backward. The process is completely identical to that observed for solitons of the KdV equation. Numerical modeling of the periodic solution (1.12), (1.19) for two values D = 0.01 and D = 0.49, for which Figure 3 was constructed, showed that the solution form does not visually change over more than 100 oscillation periods.

Discussion
The geometric series method is designed to search for exact DDE solutions, and the solutions found are expressed in terms of the ratio of two polynomials in powers of the exponential function. As repeatedly noted by leading experts [25], there is no universal approach to finding the exact solutions of nonlinear equations. The closest to the geometric series method should be considered the method of exponential functions (exp-function method) [25][26][27][28] and the method of hyperbolic tangent (tanh-method) [16,17,29]. The exp-function method uses an ansatz in the form of a ratio of two polynomials in powers of an exponential function with arbitrary coefficients. After substituting the ansatz in the equation, it is required that in the numerator of the obtained expression the factors at each degree of the exponential function vanish. This requirement leads to the need to solve a system of nonlinear algebraic equations for the coefficients. Unfortunately, the nonlinear system of the exp-function method is successfully solved by modern computer mathematics systems only in the simplest cases [5]. In particular, in the case of the discrete Sawada-Kottera equation under consideration, a calculation in Maple 2016.2 for one hour on a 6-core AMD Ryzen 5 4.2 GHz workstation with 16 Gb of RAM, did not lead to any results. Unlike the exp-function method, the geometric series method is based on the sequential solution of a system of linear equations and the exact solution (10) of the discrete Sawada-Kotera equation is successfully detected within a few seconds.
The tanh-method uses an ansatz in the form of a polynomial in powers of the hyperbolic tangent, the highest degree of which is consistent with the pole order of the desired solution. If tanh(z) is expressed in terms of exp(z), then the tanh-method ansatz is represented by the ratio of two polynomials in powers of exp(z), as in the geometric series method. However, it is easy to verify that solution (10) cannot be represented by any finite polynomial in powers of tanh(z) and, therefore, cannot be obtained in principle by the tanh-method.
Sometimes an explicit form of exact DDE solutions can be obtained using Darboux transforms of the corresponding Lax pairs [30], however, this is a much more time-consuming procedure compared to the proposed direct method.
A review of available publications on the problem could not reveal solutions of the discrete Sawada-Kotera equation, the structure of which coincides with Equation (10). It can be argued that Equation (10) is a new solution. As far as we know, the exact periodic solution (12), (19) and the initial problem solution (27) were also obtained for the first time.

Conclusions
In this article, a number of problems are solved for an integrable discrete Sawada-Kotera equation. For the first time using the geometric series method, its exact solitary-wave solution in explicit form was constructed and numerical modeling was carried out, revealing its soliton properties. Exact periodic solutions are obtained in terms of the Jacobi elliptic function. When solving the problem with special initial conditions, new combinatorial dependencies were identified that connect the desired solution with the generating function for Catalan numbers. So far, questions remain about the possibility of constructing rational solutions, as well as exact solutions of DDE systems by the geometric series method. In the future, it may be able to study the problems of classifying exact solutions of integrable DDEs and constructing quasi-exact and approximate solutions of some non-integrable equations.