Integrability test for evolutionary lattice equations of higher order

Journal of Symbolic Computation 74 (2016) 125–139

Contents lists available at ScienceDirect

Journal of Symbolic Computation www.elsevier.com/locate/jsc

Integrability test for evolutionary lattice equations of higher order V.E. Adler L.D. Landau Institute for Theoretical Physics, Ak. Semenov str. 1-A, 142432 Chernogolovka, Russian Federation

a r t i c l e

i n f o

Article history: Received 9 February 2015 Accepted 5 May 2015 Available online 14 June 2015 Keywords: Volterra type lattice Higher symmetry Conservation law Integrability test Summation by parts Computer algebra

a b s t r a c t A generalized summation by parts algorithm is presented for solving of difference equations of the form T m ( y ) − ay = b with variable coefficients, where T denotes the shift operator. Solvability of equations of this type with respect to the coefficients of formal symmetry (or formal recursion operator) provides a convenient integrability test for evolutionary differential-difference equations u ,t = f (u −m , . . . , um ). © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Let F be the set of locally analytic functions depending on a finite number of variables u j ∈ K, K = R or C, j ∈ Z, and let the shift operator T act on elements of F according to the rule

T ( f (u i , . . . , u j )) = f (u i +1 , . . . , u j +1 ). Our goal in this article is to develop an algorithm which allows one either to construct a solution y ∈ F of the linear difference equation

T m ( y ) − ay = b,

(1)

with an integer exponent m ≥ 1 and given coefficients a, b ∈ F , or to prove that the solution does not exist. In the latter case, the algorithm should return an obstacle, that is, some non-zero expression in terms of the coefficients which would vanish if there was a solution. An important special case

E-mail address: [email protected]. http://dx.doi.org/10.1016/j.jsc.2015.05.008 0747-7171/© 2015 Elsevier Ltd. All rights reserved.

126

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

( T − 1)( y ) = b

(2)

can be reformulated as the inversion problem for the operator of total difference T − 1. It was addressed by many authors in the context of discrete calculus of variations (Kupershmidt, 1985; a parallel conHydon and Mansfield, 2004; Mansfield and Quispel, 2005), see also Olver (1993) for  T − j ∂ j is the tinuous theory. In particular, it is well known that K ⊕ im( T − 1) = ker E, where E = difference Euler operator or the variational derivative. The preimage of T − 1 can be computed by use of the so-called summation by parts algorithm or by use of a discrete homotopy operator (Hereman et al., 2006). A variational interpretation of more general operators T m − a and generalization of discrete homotopy operator remain open questions, very interesting from theoretical standpoint; however, this approach can hardly give an effective solution method for equation (1). A motivation for study of equation (1) comes from the theory of differential-difference (or lattice) evolutionary equations of the form

∂t (un ) = f (un−m , . . . , un+m ),

n ∈ Z.

(3)

Recall, that such an equation is called integrable if it is consistent with an infinite set of evolutionary flows of higher orders, or symmetries. Although this property may be not so easy to check intermediately, one can deduce that it implies solvability of certain sequence of equations of the form (1)

T m(g j ) − a j g j = b j ,

j = 0, −1, . . . ,

(4)

where a j and b j are computed explicitly if g 0 , . . . , g j +1 are already known (see Proposition 17). This gives us a test which amounts to stepwise checking of whether equation (4) is solvable with respect to g j : if not, then equation (3) is not integrable, if so, then we have to compute g j and to go to the next condition for g j −1 . In practice, this test turns out to be very useful, although checking of infinite number of conditions is formally needed in order to prove the integrability. Moreover, the sequence of necessary integrability conditions encoded in relations (4) can be used for classification of the whole set of equations under consideration (this problem is certainly much more difficult than testing of a given equation and may require additional assumptions such as the existence of an infinite set of higher order conservation laws). In the continuous setup, this approach made it possible to solve a number of classification problems for integrable partial differential equations of Korteweg–de Vries and nonlinear Schrödinger type, see e.g. Sokolov and Shabat (1984), Mikhailov et al. (1987), Mikhailov et al. (1991), Mikhailov and Shabat (1993), Meshkov and Sokolov (2012). An exhaustive classification of integrable equations (3) at m = 1 (Volterra type lattices) was obtained by Yamilov (1983), but only few examples are known at m > 1 so far, the Bogoyavlensky lattices (Bogoyavlensky, 1991) being the most well studied ones. Several classification results for other types of lattice equations were obtained by Shabat and Yamilov (1991), Adler and Shabat (1997), Adler et al. (2000), Yamilov (2006), Adler (2008). The problems of symbolic computation of the higher symmetries, conservation laws, recursion operators and Lax pairs were discussed in many papers, see e.g. Göktas¸ and Hereman (1999), Hickman and Hereman (2003), Hereman et al. (2005), Sokolov and Wolf (2001), Tsuchida and Wolf (2005); integrability tests based on these notions were developed e.g. in Gerdt et al. (1985), Gerdt (1993), Hereman et al. (1998). In the problem of solving equation (1), the main issue is related with the complexity of the coefficients a, b. In particular, the right hand side in (4) becomes more and more involved at each step, so that applying the test may turn to be non-trivial. A simple method described in Section 2.2 makes use of the expansion ( T m − a)−1 = T −m (1 + aT −m + (aT −m )2 + . . .) and allows us to obtain the general solution of (1) in explicit form, if it exists. This approach is rather straightforward, but, unfortunately, it is practically applicable only if the coefficients are not too complicated. Section 2.3 contains a more effective ‘generalized summation by parts algorithm’, the main result of the article. In short, it is based on simplifications of equation by a sequence of suitable substitutions which exist under certain relations between a and b. If all relations are fulfilled then we construct the solution after a finite number of steps and if some relation fails then this proves that the solution does not exist. Applications to the integrability problem for the lattice equations are given in Section 3.1. It contains some basic notions and theorems of the symmetry approach which are necessary in order to

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

127

derive equations (4) and to describe the procedure of testing of a given lattice equation (3). Several simple examples are discussed in Section 3.2. An implementation of the algorithm in the Mathematica programming language and some further examples are available at Adler (2014b). 2. Algorithms for solving of equation (1) 2.1. Notations The partial derivatives with respect to the dynamical variables will be denoted ∂ j = ∂/∂ u j , f ( j ) = ∂ j ( f ), f ( j ,k) = ∂ j ∂k ( f ). Notice the identity T k ∂ j = ∂ j +k T k . Let us define the orders of a function f ∈ F as follows:

ord f = min{ j : f ( j ) =  0}, ord f = +∞,

ord f = max{ j : f ( j ) =  0},

ord f = −∞,

f = const,

f = const .

(5) (6)

We will also use the notation

J ( f ) = [ord( f ), ord( f )],

J (const) = ∅,

where [ p , q] = { j ∈ Z : p ≤ j ≤ q}. Notice, that although these operations look very elementary, their definition relies actually on the possibility to compare a given expression with zero. It is well known that this can be a difficult problem, for instance for the algebraic functions, implicit functions and so on, see e.g. Takayama (1992). We will always assume that functions under consideration belong to a class which admits the algorithmic zero recognition. The following operation (extraction) will be used in the summation by parts algorithm:

g = Xk ( f ) :

g (k) = f (k) ,

g ( j ) = 0, j = k, j ∈ / J ( f (k) ).

(7)

It is defined up to addition of an arbitrary function of variables u j , j ∈ J ( f (k) ) \ {k}. Function X k ( f ) can be computed by the formula

 Xk ( f ) =

f (k) duk ,

(8)

where we assume that the integrand is cast in a form not containing variables u j at j ∈ / J ( f (k) ) explicitly and that the integration constant does not depend on these variables either (which is natural). Of course, the applicability of this definition depends on the class of functions under consideration and also on the properties of the integrator used. Alternatively, one can accept the definition

X k ( f ) = f |u j =c j , j =  k, j ∈ / J ( f (k) ) ,

(9)

where c j are any constants such that the expression in the right hand side is correctly defined. 2.2. Formal inversion of operator T m − a Both algorithms presented in the paper depend on the form of the coefficient a. First of all, we notice that it governs the uniqueness of solution. Proposition 1. If a = T m (h)/h then ker( T m − a) = Kh, otherwise ker( T m − a) = 0. Proof. Assume that equation (1) possesses two different solutions y and y˜ , then a = T m ( y − y˜ )/ ( y − y˜ ). 2

128

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

1) Case a = const = 0. Proposition 2. Let b = const, J (b) = [ p 1 , p 2 ], r = ( p 2 − p 1 )/m . If a solution y exists then it is of the form

y=c+

r  s =1



as−1 T −sm (b)u =c , j < p , 1 j j

r=

p − p  2 1 m

,

(10)

where c j are arbitrary constants such that the sum in the right hand side is well defined and c is an undetermined constant (arbitrary if a = 1). Proof. Let y exist, then ord y = p 1 and ord y = p 2 − m (notice, that this implies p 2 ≥ p 1 + m). Let us consider the relation

y − ar T −rm ( y ) = T −m (b) + · · · + ar −1 T −rm (b),

r≥1

which follows from equation (1). The inequality p 1 > p 2 − m − rm holds for the chosen value of r, therefore y and T −rm ( y ) depend on disjoint sets of variables and we can get rid of the second function by fixing suitable constant values of the arguments. 2 According to Proposition 1, different choices of constants c j in (10) may only result in different representations of the constant c in the final answer. It is clear that if b is a polynomial then the simplest choice is to set all c j = 0. If a = 1, then the alternative formula

y = 1/(1 − ar )

r  s =1



as−1 T −sm (b)u =u j j +rm , j < p 1

(11)

avoids the problem of the choice of constants. Proposition 2 does not cover b = const, but in this case the solution is obvious:

y=c

at a = 1, b = 0,

y

at a = 1, b = 0,

y = b/(1 − a) at a = 1, where c is an arbitrary constant. The existence of the solution is checked by a direct substitution of expression (10) or (11) into equation (1). Alternatively, this can be done before the actual computation of solution by checking of whether b belongs to the image of T m − a characterized as the common kernel of Euler type operators. Proposition 3. If a = const = 1 then

im( T m − a) =

m −1 

ker

j =0



as T −sm ∂ j +sm ,

s∈Z

if a = 1 then

K ⊕ im( T m − 1) =

m −1 

ker

j =0



T −sm ∂ j +sm .

s∈Z

Proof. Let us differentiate equation (1) with respect to u j :

T m ( y ( j −m) ) − ay ( j ) = b( j ) . Summation yields the relations

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139



as T −sm (b( j +sm) ) = 0,

j = 0, 1, . . . , m − 1

129

(12)

s

which, therefore, serve as necessary solvability conditions for equation (1). that J (b) = [ p 1 , p 2 ] and equations (12) are fulfilled. This means that for any r Conversely, assume  r s −sm (b) cancels except for some boundary terms, that is, large enough the sum s=0 a T

b + · · · + ar T −rm (b) = f (u q2 , . . . , u p 2 ) − g (u p 1 −rm , . . . , u q1 −rm ), where q1 , q2 are some fixed numbers such that q1 > p 1 and q2 < p 2 . From here, we obtain

b − f + aT −m ( f ) = ar +1 T −(r +1)m (b) − g + aT −m ( g ). Since the two sides of the equality depend on disjoint sets of variables, their common value is constant, that is, the equation T m ( y ) − ay = b + const is solvable. Therefore, if a = 1 then relations (12) give an exact definition of im( T m − a) and if a = 1 then these conditions characterize the set K ⊕ im( T m − 1). 2 2) Case a = α T m (h)/h, α = const. This case is brought to the previous one by the change y = h y˜ . The question, whether the function a is of such form, is equivalent, upon the substitution z = log h, λ = log α , to solvability of auxiliary equation T m ( z) − z = log a − λ with the undetermined parameter λ. 3) Case a = α T m (h)/h. Differentiation of equation (1) with respect to u j yields

T m ( y ( j −m) ) − ay ( j ) = a( j ) y + b( j ) and elimination of derivatives of y results in the relation



as T sm (a( j −sm) y + b( j −sm) ) = 0,

(13)

s

where the sum contains only finite number of non-vanishing terms corresponding to the values of s from the interval

 j − max(ord a, ord b)  m

≤s≤

 j − min(ord a,ord b)  m

.

The coefficients as are defined by the recurrence relation a s−1 = as T sm (a), a0 = 1, that is,

as =

0

T km (a), s ≤ 0,

a s = 1/

k = s +1

s

T km (a), s ≥ 0.

k =1

Equation (13) is obviously a generalization of (12) for the case of non-constant a. The equation corresponding to j = j + m differs from (13) only by a factor. Elimination of T sm ( y ) by use of equation (1) itself allows us to bring (13) to the form

A j y = B j,

j = 0, 1 , . . . , m − 1 .

Proposition 4. The factors A j are of the form A j = value of j.

(14)

 s

T sm (∂ j −sm (log a)) and do not vanish for at least one

130

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

Proof. By induction, it follows from (1) that as−1 T sm ( y ) = ay + (terms with b). Substitution into (13) implies the formula for A j . If all A j are identically zero then log a ∈ K ⊕ im( T m − 1) according to Proposition 3, but this contradicts the assumption a = α T m (h)/h. 2 Thus, if the solution y exists then it is found explicitly from one of equations (14) and we only have to substitute it into (1) in order to make the check. Remark 5. The presented analysis proves that if a solution of equation (1) exists then it is a rational function of the coefficients a, b and functions obtained from a, b by means of the operators T , ∂ j and substitutions u j = const. In particular, if a, b are rational functions of the variables u j then the solution is rational as well. If a = const and b is a polynomial then the solution is also a polynomial. 2.3. Generalized summation by parts algorithm The above approach is rather clear and can be easily implemented in the computer algebra systems. Unfortunately, computing of sums like (10) or (13) is not too effective in practice (especially if the solution does not exist, since in this case we cannot expect cancellation of intermediate terms). An alternative algorithm makes use of a sequence of suitable substitutions of the form

y = y˜ / A ,

a˜ = aT m ( A )/ A ,

b˜ = T m ( A )b

(15)

or

y = y˜ − B ,

a˜ = a,

b˜ = b + ( T m − a)( B )

(16)

which bring (1) to equivalent equations of the same type, but with coefficients depending on a reduced set of variables. As a result of a finite number of steps, we either construct the solution y explicitly or prove that it does not exist. The substitutions can be applied in different ways, but the final answer does not depend on this, because of their reversibility. We will organize the flow of control by use of inequalities involving the orders

ord a = q1 ,

ord a = q2 ,

ord b = p 1 ,

ord b = p 2 .

Some inequalities imply the existence of required substitutions (15), (16), with functions A, B obtained from a, b by the extraction operation (7). Other inequalities mean that equation (1) is unsolvable, in such a case the algorithm returns some non-zero expression which plays the role of an obstacle for the existence of the solution. The analysis of such obstacles is important in a situation when the equation coefficients contain arbitrary parameters. 1) Case a = const = 0. If a solution y exists then J ( y ) = [ p 1 , p 2 − m]. If the inequality

p1 > p2 − m holds then solution may be only constant, that is, y = b/(1 − a) if a = 1 and y = c if a = 1, and we only have to check this by inspection. Notice, that the arbitrary constant c here is the only source of possible non-uniqueness for the whole algorithm, upon all substitutions below. Proposition 6. Let p 1 ≤ p 2 − m. If the inequality

ord b( p 2 ) ≥ p 1 + m

(17)

fails then the solution does not exist. If the inequality holds then the change

B = X p 2 (b),

y = y˜ + T −m ( B ),

b˜ = b − B + aT −m ( B )

brings (1) to an equivalent equation with J (b˜ ) ⊆ [ p 1 , p 2 − 1].

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

131

Proof. Assume that solution y exists, then differentiation of (1) with respect to u p 2 yields the relation b( p 2 ) = T m ( y ( p 2 −m) ) which implies inequality (17). Let it be fulfilled, then J (b( p 2 ) ) ⊆ [ p 1 + m, p 2 ] and therefore our substitution reduces J (b) (the case J (b˜ ) = ∅, that is, b˜ = const, is possible). 2 If condition (17) fails then the non-zero expression b(r , p 2 ) , where r = ord b( p 2 ) , is returned as an obstacle for the existence of the solution. By applying Proposition 6 while possible (no more than p 2 − p 1 + 1 times), we will either construct the solution as a finite sum y = T −m ( B + B˜ + . . .) or prove that it does not exist. Example 7. Let m = 1, b = cu 10 − u 0 . Then B = cu 10 and we apply the transformation

y = y˜ + cu 9 ,

b˜ = b − B + aT −1 ( B ) = acu 9 − u 0 .

˜

On the next step, we get y˜ = y˜˜ + acu 8 , b˜ = a2 cu 8 − u 0 , and so on. In ten steps, we find that if c = a−10 then y = a−10 u 9 + a−9 u 8 + · · · + a−1 u 0 is the solution, and if c = a−10 then the solution does not exist. 2) Reduction of the coefficient a. Proposition 8. Let a = const and

q1 ≤ q2 − m

and

ord ∂q1 (log a) ≤ q2 − m

(18)

(possibly, ∂q1 (log a) = const). Then the substitution

A = exp( X q1 (log a)),

y = y˜ / A ,

a˜ = aT m ( A )/ A ,

b˜ = T m ( A )b

reduces (1) to an equivalent equation with J (˜a) ⊆ [q1 + 1, q2 ] (possibly, with a˜ = const). Proof. Conditions (18) mean that a is of the form

a = A (u q1 , . . . , u q2 −m )ˆa(u q1 +1 , . . . , u q2 ).

2

Repetition of this transformation (at most q2 − q1 + 1 times) will lead either to the case 1) or to the case when one of inequalities (18) fails, that is,

max(q1 , ord ∂q1 (log a)) > q2 − m.

(19)

We will assume that this condition is fulfilled from now on. Further substitutions will not change a. Example 9. Let a = u 1 u 2 u 3 and m = 1, then the resulting transformation will be y = y˜ /(u 1 u 22 ), a˜ = u 33 . If m = 2 then we apply the change y = y˜ /u 1 , a˜ = u 2 u 23 , and if m = 3 then no change is applied. 3) Reduction of the coefficient b. The following statement is obvious (cf. with Proposition 6). Proposition 10. If p 2 > q2 then the substitution

B = X p 2 (b),

y = y˜ + T −m ( B ),

b˜ = b − B + aT −m ( B )

brings (1) to an equivalent equation with ord b˜ < p 2 .

(20)

132

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

Clearly, repetition of this substitution (p 2 − q2 times maximum) brings the problem to the case p 2 ≤ q2 . Notice, that instead of (20) we can apply a simpler change

b˜ = aT −m (b)

y = y˜ + T −m (b),

with the same effect. The advantage of (20) is that it drops the lower order p 1 not as much (actually, the next step of the algorithm will be to increase p 1 if possible). Example 11. Let m = 1, a = u 0 and b = 1/(u 1 u 2 u 3 ) − u −3 u −2 u −1 u 0 . Then the transformation (20) reads

B=

1 u1 u2 u3

,

y = y˜ +

1 u0 u1 u2

,

b˜ =

1 u1 u2

− u −3 u −2 u −1 u 0 .

Since the condition p 2 > q2 for b˜ is still fulfilled, we repeat the change and reduce the problem to the case

y = y˜ +

1 u0 u1 u2

+

1 u0 u1

+

1 u0

,

b˜ = 1 − u −3 u −2 u −1 u 0 .

Proposition 12. Let p 2 ≤ q2 and p 1 < q1 . If the inequality

ord(b( p 1 ) /a) ≤ q2 − m

(21)

fails then the solution does not exist. If the inequality holds then the substitution

B = X p 1 (b/a),

y = y˜ − B ,

b˜ = b + ( T m − a)( B )

brings (1) to an equivalent equation with J (b˜ ) ⊆ [ p 1 + 1, q2 ]. Proof. If a solution y exists then ord y ≤ q2 − m and − y ( p 1 ) = b( p 1 ) /a from where condition (21) follows. If it holds then J (b( p 1 ) /a) ⊆ [ p 1 , q2 − m] and our substitution increases p 1 while leaving p 2 ≤ q2 . 2 If condition (21) fails then the expression ∂r (b( p 1 ) /a), where r = ord(b( p 1 ) /a), serves as an obstacle for the solution existence. Applying Proposition 12 while possible (q1 − p 1 times or less) brings to the final step of our algorithm. Example 13 (Continuation). We have brought the previous example to the case m = 1, a = u 0 , b = 1 − u −3 u −2 u −1 u 0 . Since b( p 1 ) /a = −u −2 u −1 and condition (21) is fulfilled, we apply the substitution

B = − u −3 u −2 u −1 ,

y = y˜ + u −3 u −2 u −1 ,

b˜ = 1 − u −2 u −1 u 0 .

˜ the condition (21) remains valid and we arrive to For b, y = y˜ + u −3 u −2 u −1 + u −2 u −1 + u −1 ,

b˜ = 1 − u 0

by applying the transformation two times more. 4) Solving of a linear system. The problem is reduced now to the case J (b) ⊆ J (a) = [q1 , q2 ]. Proposition 14. Assume that the solution y exists. If q1 > q2 − m then y = b/(1 − a) = const and if q1 ≤ q2 − m then y satisfies the linear system with non-vanishing determinant with respect to y, y (q1) :

a(q1 ) y

+ ay (q1 )

= −b(q1 ) ,

a(q1 ,r ) y

+

= −b(q1 ,r ) ,

a(r ) y (q1 )

where r = ord ∂q1 (log a).

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

133

Proof. The orders of solution y are J ( y ) ⊆ [q1 , q2 − m]. Therefore, if q1 > q2 − m then equation (1) may admit only a constant solution. If q1 ≤ q2 − m then condition (19) implies that r > q2 − m ≥ ord y, and the linear system follows from (1) by differentiating with respect to u q1 and then u r . 2 According to this proposition, if a solution exists then it is given by an explicit expression and we only have to check it by substitution into equation (1). Example 15 (Completion). The previous example was reduced to the equation T ( y ) − u 0 y = 1 − u 0 . Here q1 = q2 = 0 and we obtain y = b/(1 − a) = 1 which is the solution indeed. Collecting all together, we obtain the solution y = u −3 u −2 u −1 + u −2 u −1 + u −1 + 1 + 1/u 0 + 1/(u 0 u 1 ) + 1/(u 0 u 1 u 2 ) of the original equation T ( y ) − u 0 y = 1/(u 1 u 2 u 3 ) − u −3 u −2 u −1 u 0 . Summing up, the finite sequence of substitutions described in Propositions 6, 8, 10, 12 and the explicit recipe given in Proposition 14 allow us to construct the solution of equation (1) or to prove its non-existence. The whole structure of the algorithm can be represented by a binary tree with either recursive calls or breaks on the leaves and all branches governed by inequalities in terms of the orders of functions under consideration. Moreover, all substitutions use only the arithmetic operations, shift, partial derivatives and extraction operation X k . In turn, operations ord, ord and X k are defined in terms of differentiation and comparison with zero. Therefore, the presented algorithm is applicable to equations with the coefficients in any differential field which admits the algorithmic zero recognition. In practice, the most important functional class with this property consists of expressions which are rational with respect to u j and indeterminate functions of u j . In this class, operation X k can be effectively realized by the formula (8). The use of elementary, special or implicit functions is restricted by our ability to handle the necessary identities for such functions, see examples in Adler (2014b). 3. Integrability test 3.1. The symmetry approach Let us recollect some basic notions of the symmetry approach in application to the scalar evolutionary lattice equations

∂t (un ) = f (un−m , . . . , un+m ),

n∈Z

or, in a shorthand notation,

u ,t = f (u −m , . . . , um ).

(22)

A detailed exposition can be found in Mikhailov et al. (1987), Mikhailov et al. (1991), Levi and Yamilov (1997), Yamilov (2006). For any function f ∈ F , the infinite-dimensional vector field

Dt = ∇ f =



T j ( f )∂ j

j ∈Z

is called the evolutionary derivative and the difference operator

f∗ =



f ( j) T j

j ∈Z

is called the linearization operator. The differentiation D t in virtue of equation (22) is defined in two equivalent ways D t ( g ) = ∇ f ( g ) = g ∗ ( f ). A lattice equation

u ,τ = g (u −l , . . . , ul )

134

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

is called a (generalized) symmetry of equation (22) if differentiations D t , D τ commute, that is, the relation

∇ f ( g) = ∇g ( f )

(23)

(or, equivalently, g ∗ ( f ) = f ∗ ( g )) holds identically with respect to u j . Equation (22) is considered integrable if it admits symmetries of arbitrarily large order l. Equation (23) yields, upon the linearization, a more convenient operator relation

∇ f ( g ∗ ) = ∇ g ( f ∗ ) + [ f ∗ , g ∗ ].

(24)

The degree m of the operator f ∗ is fixed and this allows us to consider g ∗ as an approximate solution of the Lax equation

D t (G ) = [ f ∗ , G ].

(25)

Let us consider the division rings

F (( T −1 )) = {



a j T j | a j ∈ F },

F (( T )) = {

j <+∞



a jT j | a j ∈ F}

j >−∞

of the formal Laurent series with respect to the negative or positive powers of T . It can be proven that the existence of a sequence of symmetries of arbitrarily large orders implies that (25) admits a solution from F (( T −1 ))

G = g k T k + · · · + g 1 T + g 0 + g −1 T −1 + . . . ,

g j ∈ F,

k>0

which is called the formal symmetry, or formal recursion operator, of lattice equation (22). Conditions of solvability of equation (25) with respect to the coefficients g j serve therefore as necessary integrability conditions for (22). A weak point here is that the degree k of G is not known in advance. It should be mentioned that in the continuous setting (for KdV type equations) the shift operator T is replaced by the total derivative D and the formal symmetry becomes a pseudo-differential operator G = gk D k + · · · + g 0 + g −1 D −1 + · · · . For such series, the operation of root extraction G → G 1/k is correctly defined and acts on solutions of (25), so that we can take k = 1 without loss of generality. In the difference situation this argument does not work: it is clear already from the fact that the existence of a root requires that the leading coefficient have the special form gk = hT (h) · · · T k−1 (h). Nevertheless, the following theorem from (Adler, 2014a) states that the degree k can always be chosen equal to the order m of equation (22) itself (although this degree may be not minimal). The proof is based on the observation that if G is a solution of (25) of degree k then another solution H of degree m can be computed explicitly by use of two relations D t ( H ) = [ f ∗ , H ] and H k = G m simultaneously, after which we can prove that each relation is satisfied separately. Theorem 16. If lattice equation (22) admits symmetries of arbitrarily large order then the Lax equation (25) admits a solution of the form

G = f (m) T m + · · · + f (1) T + g 0 + g −1 T −1 + . . . ∈ F (( T −1 )).

(26)

Now, equation (25) turns into an explicit and convenient test, since the resulting necessary integrability conditions do not depend on the actual orders of higher symmetries and can be written down intermediately from the right hand side of equation (22). Proposition 17. If equation (22) is integrable then the following recurrence relations are solvable with respect to g j ∈ F :

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

T m(g j) − a j g j = b j, bj =

T j ( f (m) )

aj =

1

D t ( g j +m ) − f (m)

f (m)

m −1 



,

135

j = 0, −1, −2, . . . ,

(27)

 ,

f (s) T s ( g j +m−s ) − g j +m−s T j +m−s ( f (s) )

(28)

s=−m

where it is assumed that gk = 0 for k > m and gm = f (m) , . . . , g 1 = f (1) . The proof follows immediately from collecting terms with T j +m in (25). The expression for b j involves only coefficients g 0 , . . . , g j +1 which are already computed. Thus, the integrability test amounts to the step-by-step checking of solvability of equations (27) which are of the form (1). Remark 18. A bit more detailed analysis of the limiting transition from equation (24) to (25) shows that the existence of a symmetry of order l ≥ m + r implies that first r equations (27) can be resolved with respect to g 0 , . . . , g −r +1 . Symmetries of orders l ≤ m give no conditions in this approach, being lost on the background of the trivial symmetry u ,τ = f . Concerning the sufficiency, the fulfilment of first r conditions (27) does not formally guarantee the existence of even one generalized symmetry; however, if r is large enough then it is a very strong evidence of integrability. In the theory of lattice equations, the shifts T and T −1 are on equal footing. This means that, in addition to (27), (28), there is a complementary sequence of conditions corresponding to the formal symmetry in F (( T )). Moreover, these conditions can be considered also in a more general situation for the lattice equations

u ,t = f (u −m¯ , . . . , um ),

¯ >0 m, m

¯ In such a case, the existence of symmetries of lower order −¯l, where ¯l is arbitrarily large, with m = m. implies that the Lax equation (25) is solvable with respect to the formal symmetry of the form G¯ = f (−m¯ ) T −m¯ + · · · + f (−1) T −1 + g¯ 0 + g¯ 1 T + . . . ∈ F (( T )). If the lattice equation admits conservation laws D t (ρ ) = ( T − 1)(σ ) of arbitrarily large orders then ¯ = m and the two formal symmetries are related by the formula G¯ † = − R G R −1 . Here, necessarily m † denotes the conjugation (aT j )† = T − j a and the series

R = rl T l + rl−1 T l−1 + . . . ∈ F (( T −1 )),

0≤l
satisfies the relation †

D t ( R ) + f ∗ R + R f ∗ = 0.

(29)

It is clear that equations for g¯ j and r j are also of the type (1) and can be checked analogously to (27). How many conditions we have to check in order to guarantee the integrability of a given equation? In general, the answer is not known. At m = 1, the integrable Volterra type lattice equations

u ,t = f (u −1 , u , u 1 )

(30)

were classified by Yamilov (1983, 2006). He proved that existence of few higher symmetries and conservation laws implies

D t (log f (1) ) ∈ im( T − 1), D t (a) + 2 f (0) ∈ im( T − 1),

af (1) + T (a) f (−1) = 0, (31)

where a = a(u −1 , u ), and obtained a finite list of equations (30) by resolving these conditions. Technically, (31) can be easily derived from equations (25) and (29). Moreover, Yamilov proved that each equation in the list admits an infinite hierarchy of higher symmetries, that is, (31) are, in fact, the necessary and sufficient integrability conditions for (30).

136

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

In view of this result, checking of conditions (31) is all we need at m = 1. However, it is worth noticing that in this case the whole sequence of conditions (27) can be brought to the form ( T − 1)( y j ) = b˜ j by substitution g j = T j −1 ( f (1) ) · · · f (1) y j . Moreover, there exists a more complicated, but still invertible substitution which allows us to rewrite these conditions in the form of conservation laws (possibly, trivial)

D t (ρ j ) = ( T − 1)(σ j ),

j ≥ 0,

(32)

where the so-called canonical density ρ j is equivalent modulo im( T − 1) to j −1 coef T 0 G j , j > 0. Although this form gives no essential advantage when testing a particular equation, it clarifies a general structure of the integrability conditions. Proposition 19. At m = 1, the solvability of equations (27), (28) is equivalent to solvability with respect to

σ j ∈ F of conservation laws (32), where densities ρ j are defined by the recurrence relations

ρ0 = log f (1) , ρ1 = f (0) + σ0 , P j +1 (−ρ1 , . . . , −ρ j +1 ) + f (−1) T −1 ( f (1) P j −1 (ρ1 , . . . , ρ j −1 )) + σ j = 0,

j > 0,

with polynomials P j = P j (ρ1 , . . . , ρ j ) defined by the generating function

P 0 + P 1 λ + P 2 λ2 + . . . = exp(ρ1 λ + ρ2 λ2 + ρ3 λ3 + . . .). The proof can be found in Adler (2014a). The first few polynomials P j are

P 0 = 1,

P 1 = ρ1 ,

P 2 = ρ2 +

ρ12 2

,

P 3 = ρ3 + ρ1 ρ2 +

ρ13 6

, ...

and the corresponding conserved densities are

1

ρ2 = f (−1) T −1 ( f (1) ) + ρ12 + σ1 , 2

1

ρ3 = f (−1) T −1 ( f (1) ρ1 ) + ρ1 ρ2 − ρ13 + σ2 , 6

1

1

1

1

2

2

2

24

ρ4 = f (−1) T −1 ( f (1) (ρ2 + ρ12 )) + ρ1 ρ3 + ρ22 − ρ12 ρ2 +

ρ14 + σ3 .

In the general case m > 1, only part of conditions (27) can be rewritten as conservation laws. For instance, it is easy to prove directly from equation (25) that its solvability implies the existence of functions σ , σ1 ∈ F such that

D t (log f (m) ) = ( T m − 1)(σ ),

D t ( f (0) + σ ) = ( T − 1)(σ1 ).

3.2. Examples In this section we present several simple examples, in order to clarify various computational aspects rather than to obtain new results. The Matematica programs available at Adler (2014b) include an extended set of examples, the full Yamilov’s list of first order integrable equations, and several higher order equations. Example 20. Let us solve equations (27), (28) for the Volterra lattice

u ,t = u (u 1 − u −1 ). In this case computation of formal symmetry encounters no obstacles and setting all integration constants to zero yields

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

g1 = u ,

g0 = u + u 1 ,

g −1 =

It is easy to see that the series G =



uu 1 u −1

,

gj =

u ( u 1 − u −1 ) uj

,

137

j < −1.

g j T j can be rewritten in a closed form

G = uT + u + u 1 + uT −1 + u (u 1 − u −1 )( T − 1)−1

1 u

which is the well known recursion operator for the Volterra lattice. However, in most cases expressions for g j are much more complicated and search of corresponding recursion operators is a very nontrivial problem. For instance, in the case of the second order Bogoyavlensky lattice

u ,t = u (u 2 + u 1 − u −1 − u −2 ), we get

g2 = u , g −2 = g −4 =

g1 = u , 1

u −2

g0 = u + u 1 + u 2 ,

(u −1 u 1 + uu 1 + uu 2 ),

1 u −4 u −2

g −1 = 0,

g −3 = −

1 u −3

(u −2 u + u −1 u + u −1 u 1 ),

((u −3 + u −2 )u −1 u 1 + u −2 u (u 1 + u 2 )),

...

which gives little hint on the factored form of G

G = u (1 + T −1 + T −2 )( T 2 u − uT −1 )( T u − uT −1 )−1 ( T u − uT −2 )(u − uT −2 )−1 . This is a particular example of recursion operators found by Wang (2012) for the Bogoyavlensky lattices of any order m. Notice, that in these operators all inverse factors are binomial and therefore computation of G ( f ) for a given function f amounts to solving of a sequence of relations of the type (1). Example 21. As a sample ‘classification’ problem, consider an equation of the Bogoyavlensky type

u ,t = u (u 2 + k1 u 1 + k2 u + k3 u −1 + k4 u −2 ) with undetermined coefficients. Application of the formal symmetry test yields on the first step the obstacle

(1 + k2 + k4 )u + (k1 + k3 )u 1 = 0

⇒

k2 = −1 − k4 ,

k3 = −k1 .

After the substitutions, g 0 is successfully found, but computing of g −1 encounters the next obstacle

k1 (k1 + k4 )(u − u 2 ) + k1 (k1 − 1)(u −1 − u 1 ) = 0. If k1 = 0 then k1 = 1, k4 = −1 and we arrive at the Bogoyavlensky lattice. If k1 = 0 then computation of g −2 brings to the obstacle

k4 (1 + k4 )/u −2 = 0 and we get two more integrable (albeit disappointing) cases: a linearizable equation u ,t = u (u 2 − u ) and the stretched Volterra lattice u ,t = u (u 2 − u −2 ). Of course, this example is over-simplified. In fact, solving of any realistic classification problem, like the analysis of conditions (31), requires quite different techniques and our algorithm should be viewed just as an auxiliary tool for verifying the answers obtained during the classification.

138

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

Example 22. According to Yamilov (2006), the lattice equation

u ,t = h(u 1 − u ) + h(u − u −1 ) is integrable if h satisfies equation

h = αh2 + β h + γ

(33)

with arbitrary constant coefficients. Equation (33) can be solved in elementary functions, but this leads to consideration of several cases corresponding to different parameter sets and special solutions. In order to handle the whole family in a uniform manner we only have to compute the coefficients b j (28) modulo a rule which replaces the derivatives of h in virtue of (33). This allows us to perform computations in the field of rational functions on h(u j +1 − u j ) where the basic operations of orders (5) and extraction (7) are correctly defined. The summation by parts algorithm confirms the integrability of the lattice equation for arbitrary α , β , γ . Example 23. In the above examples, equations pass the test for any choice of integration constants, but this is not always the case. Consider the modified Bogoyavlensky lattice

u ,t = u (u 2 u 1 − u −1 u −2 )

(34)

= uu 1 . It is easy to see that the operator T − for any j, so that the general solution of equation (27) for g j contains an arbitrary constant c j on each step. However, it turns out that c −2 j +1 becomes an obstacle when we proceed to computing of g −2 j . As a result, the test passes only if we set to zero all integration constants with odd numbers. This indicates that the minimal degree of the formal symmetry G is equal to 2, so that (34) cannot be a symmetry of an equation of order 1. The same is true for equation with f (2)

2

j

T ( f (2) )/ f (2) possesses a nontrivial kernel

u ,t = (u 2 + 1)((u 21 + 1)(u 2 − u ) + (u 2−1 + 1)(u − u −2 )) which is related by the non-autonomous change un = (−1)n un to the second order symmetry of the modified Volterra lattice

u ,τ = (u 2 + 1)(u 1 − u −1 ). In contrast, this symmetry itself passes the test for arbitrary integration constants, as expected. 4. Conclusion The presented algorithm is designed for straightforward computation of the formal symmetry for a scalar evolutionary lattice equation of arbitrary order. It is suitable mainly for testing integrability of a single equation or a family depending on several parameters. Further generalizations may include equations with two or more components such as Toda or Ablowitz–Ladik type lattices and their ‘hungry’ analogs. In this case the coefficients of formal symmetry are matrices, and inversion of difference operators becomes a more difficult problem. The vectorial case (see e.g. Adler, 2008) can also be handled by enlarging the set of dynamical variables. Concerning the classification problem in general, many results were obtained in the continuous case for scalar evolutionary equations of orders 3, 5, 7, see Meshkov and Sokolov (2012) and the references therein. Moreover, there is a conjecture that all integrable higher order equations are symmetries of equations of orders 3 and 5, so that there is just a finite set of integrable hierarchies (this was proven in the homogeneous polynomial case by Sanders and Wang, 1998). In the difference case, the classification is much more difficult and the Yamilov (1983) list of first order equations remains the only rigorous result obtained so far. The known examples show that there are primitive integrable lattice equations of any orders, and their description is a challenging problem, even in the polynomial case.

V.E. Adler / Journal of Symbolic Computation 74 (2016) 125–139

139

Acknowledgements I am grateful to Evelyne Hubert and the anonymous referees for their valuable and helpful comments and suggestions for improvements. This work was supported by the Russian Foundation for Basic Research under grant 13-01-00402a and the Leading Scientific Schools Program under grant 3139.2014.2. References Adler, V.E., 2008. Classification of integrable Volterra-type lattices on the sphere: isotropic case. J. Phys. A 41 (14), 145201. Adler, V.E., 2014a. Necessary integrability conditions for evolutionary lattice equations. Theor. Math. Phys. 181 (2), 1367–1382. Adler, V.E., 2014b. A Mathematica implementation of the generalized summation by parts algorithm. Available at http:// adler.itp.ac.ru/programs/. Adler, V.E., Shabat, A.B., 1997. Generalized Legendre transformations. Theor. Math. Phys. 112 (2), 935–948. Adler, V.E., Shabat, A.B., Yamilov, R.I., 2000. Symmetry approach to the integrability problem. Theor. Math. Phys. 125 (3), 1603–1661. Bogoyavlensky, O.I., 1991. Algebraic constructions of integrable dynamical systems—extensions of the Volterra system. Russ. Math. Surv. 46 (3), 1–64. Gerdt, V.P., 1993. Computer algebra, symmetry analysis and integrability of nonlinear evolution equations. Int. J. Mod. Phys. C 4, 279–286. Gerdt, V.P., Shvachka, A.B., Zharkov, A.Yu., 1985. Computer algebra application for classification of integrable non-linear evolution equations. J. Symb. Comput. 1 (1), 101–107. Göktas, ¸ Ü., Hereman, W., 1999. Algorithmic computation of generalized symmetries of nonlinear evolution and lattice equations. Adv. Comput. Math. 11, 55–80. Hereman, W., Deconinck, B., Poole, L.D., 2006. Continuous and discrete homotopy operators: a theoretical approach made concrete. Math. Comput. Simul. 74 (45), 352–360. Hereman, W., Göktas, ¸ Ü., Colagrosso, M.D., Miller, A.J., 1998. Algorithmic integrability tests for nonlinear differential and lattice equations. Comput. Phys. Commun. 115 (2–3), 428–446. Hereman, W., Sanders, J.A., Sayers, J., Wang, J.P., 2005. Symbolic computation of conserved densities, generalized symmetries, and recursion operators for nonlinear differential-difference equations. In: Winternitz, P., et al. (Eds.), Group Theory and Numerical Methods. In: CRM Proc. Lecture Notes, vol. 39. AMS, Providence, Rhode Island, pp. 267–282. Hickman, M., Hereman, W., 2003. Computation of densities and fluxes of nonlinear differential-difference equations. Proc. R. Soc., Math. Phys. Eng. Sci. 459, 2705–2729. Hydon, P.E., Mansfield, E.L., 2004. On the variational complex for difference equations. Found. Comput. Math. 4, 187–217. Kupershmidt, B.A., 1985. Discrete Lax Equations and Differential-Difference Calculus. Asterisque, Paris. Levi, D., Yamilov, R.I., 1997. Conditions for the existence of higher symmetries of evolutionary equations on the lattice. J. Math. Phys. 38, 6648–6674. Mansfield, E.L., Quispel, G.R.W., 2005. Towards a variational complex for the finite element method. In: Winternitz, P., et al. (Eds.), Group Theory and Numerical Methods. In: CRM Proc. Lecture Notes, vol. 39. AMS, Providence, Rhode Island, pp. 207–231. Meshkov, A.G., Sokolov, V.V., 2012. Integrable evolution equations with the constant separant. Ufa Math. J. 4 (3), 104–153. Mikhailov, A.V., Shabat, A.B., 1993. Symmetries—test of integrability. In: Fokas, A., Zakharov, V. (Eds.), Important Developments in Soliton Theory. Springer-Verlag, Berlin, pp. 355–372. Mikhailov, A.V., Shabat, A.B., Sokolov, V.V., 1991. The symmetry approach to classification of integrable equations. In: Zakharov, V.E. (Ed.), What Is Integrability?. Springer-Verlag, pp. 115–184. Mikhailov, A.V., Shabat, A.B., Yamilov, R.I., 1987. The symmetry approach to classification of nonlinear equations. Complete lists of integrable systems. Russ. Math. Surv. 42 (4), 1–63. Olver, P.J., 1993. Applications of Lie Groups to Differential Equations, 2nd ed. Grad. Texts Math., vol. 107. Springer-Verlag, New York. Sanders, J.A., Wang, J.P., 1998. On the integrability of homogeneous scalar evolution equations. J. Differ. Equ. 147 (2), 410–434. Shabat, A.B., Yamilov, R.I., 1991. Symmetries of nonlinear chains. Leningr. Math. J. 2 (2), 377–399. Sokolov, V.V., Shabat, A.B., 1984. Classification of integrable evolution equations. Sov. Sci. Rev., C, Math. Phys. Rev. 4, 221–280. Sokolov, V.V., Wolf, T., 2001. Classification of integrable polynomial vector evolution equations. J. Phys. A 34 (49), 11139–11148. Takayama, N., 1992. An approach to the zero recognition problem by Buchberger algorithm. J. Symb. Comput. 14 (2–3), 265–282. Tsuchida, T., Wolf, T., 2005. Classification of polynomial integrable systems of mixed scalar and vector evolution equations. I. J. Phys. A 38 (35), 7691–7733. Wang, J.P., 2012. Recursion operator of the Narita–Itoh–Bogoyavlensky lattice. Stud. Appl. Math. 129 (3), 309–327. Yamilov, R.I., 1983. On classification of discrete evolution equations. Usp. Mat. Nauk 38 (6), 155–156 (in Russian). Yamilov, R.I., 2006. Symmetries as integrability criteria for differential difference equations. J. Phys. A 39 (45), R541–R623.