SYMBOLIC COMPUTATION OF RECURSION OPERATORS FOR NONLINEAR DIFFERENTIAL-DIFFERENCE EQUATIONS

.


INTRODUCTION
A number of interesting problems can be modeled with nonlinear differentialdifference equations (DDEs) [1]- [3], including particle vibrations in lattices, currents in electrical networks, and pulses in biological chains.Nonlinear DDEs also play a role in queuing problems and discretizations in solid state and quantum physics, and arise in the numerical solution of nonlinear PDEs.
The study of complete integrability of nonlinear DDEs largely parallels that of nonlinear partial differential equations (PDEs) [4]- [7].Indeed, as in the continuous case, the existence of large numbers of generalized (higher-order) symmetries and conserved densities is a good indicator for complete integrability.Albeit useful, such predictors do not provide proof of complete integrability.Based on the first few densities and symmetries, quite often one can explicitly construct a recursion operator which maps higherorder symmetries of the equation into new higher-order symmetries.The existence of a recursion operator, which allows one to generate an infinite set of such symmetries stepby-step, then confirms complete integrability.
There is a vast body of work on the complete integrability of DDEs.Consult, e.g., [5,8] for additional references.In this article we describe an algorithm to symboli-cally compute recursion operators for DDEs.This algorithm builds on our related work for PDEs and DDEs [9]- [11] and work by Oevel et al [12] and Zhang et al [13].
In contrast to the general symmetry approach in [5], our algorithms rely on specific assumptions.For example, we use the dilation invariance of DDEs in the construction of densities, higher-order symmetries, and recursion operators.At the cost of generality, our algorithms can be implemented in major computer algebra systems.
Our Mathematica package InvariantsSymmetries.m [14] computes densities and generalized symmetries, and therefore aids in automated testing of complete integrability of semi-discrete lattices.Our new Mathematica package DDERecursionOperator.m [15] automates the required computations for a recursion operator.
The paper is organized as follows.In Section 2, we present key definitions, necessary tools, and prototypical examples, namely the Kac-van Moerbeke (KvM) [16] and Toda [17,18] lattices.An algorithm for the computation of recursion operators is outlined in Section 3. Usage of our package is demonstrated on an example in Section 4. Section 5 covers two additional examples, namely the Ablowitz-Ladik (AL) [19] and RelativisticToda (RT) [20] lattices.Concluding remarks about the scope and limitations of the algorithm are given in Section 6.

Definition
A nonlinear DDE is an equation of the form , ,...) , , (..., where n u and F are vector-valued functions with N components.The subscript n cor- responds to the label of the discretized space variable; the dot denotes differentiation with respect to the continuous time variable .u and a finite number of its forward and backward shifts.We assume that F is polynomial with constant coefficients.No restrictions are imposed on the shifts or the degree of nonlinearity in .F

Example
The Kac-van Moerbeke (KvM) lattice [16], also known as the Volterra lattice, , ) ( arises in the study of Langmuir oscillations in plasmas, population dynamics, etc.

Example
One of the earliest and most famous examples of completely integrable DDEs is the Toda lattice [17,18], , ) exp( ) exp( where n y is the displacement from equilibrium of the nth particle with unit mass under an exponential decaying interaction force between nearest neighbors.With the change of variables, ), exp( , due to Flaschka [21], lattice (3) can be written in polynomial form [22] .) ( ,

Definition
A DDE is said to be dilation invariant if it is invariant under a scaling (dilation) symmetry.

Example
Equation ( 4) is invariant under the scaling symmetry , ) , where  is an arbitrary scaling parameter.

Definition
We define the weight, , w of a variable as the exponent in the scaling parameter  which multiplies the variable.As a result of this definition, t is always replaced by

Definition
The rank of a monomial is defined as the total weight of the monomial.An expression is uniform in rank if all of its terms have the same rank.

Example
In the first equation of (4), all the monomials have rank 2; in the second equation all the monomials have rank 3. Conversely, requiring uniformity in rank for each equation in (4) allows one to compute the weights of the dependent variables (and thus the scaling symmetry) with elementary linear algebra.Indeed, which is consistent with (5).Dilation symmetries, which are Lie-point symmetries, are common to many lattice equations.Polynomial DDEs that do not admit a dilation symmetry can be made scaling invariant by extending the set of dependent variables with auxiliary parameters with appropriate scales.

Definition
A scalar function ) ( n n u  is a conserved density of (1) if there exists a scalar function ), ( n n J u called the associated flux, such that [23] is satisfied on the solutions of (1).
In (8), we used the (forward) difference operator, where D denotes the up-shift (forward or right-shift) operator, , D and I is the identity operator.
The operator  takes the role of a spatial derivative on the shifted variables as many DDEs arise from discretization of a PDE in 1 1 variables.Most, but not all, densities are polynomial in .n u

Example
The first four density-flux pairs [22] for (4) are , ), ln( The densities in ( 13)-( 16) are uniform of ranks 0 through 3, respectively.The corresponding fluxes are also uniform in rank with ranks 1 through 4, respectively.In general, if in ( 8) The various pieces in (8) are uniform in rank.Since (8) holds on solutions of (1), the conservation law 'inherits' the dilation symmetry of (1).
Consult [22] for our algorithm to compute polynomial conserved densities and fluxes, where we use (4) to illustrate the steps.Non-polynomial densities (which are rare) can be computed by hand or with the method given in [8].

Definition
A vector function ) ( is called a generalized (higher-order) symmetry of ( 1) For the vector case with two components n u and , In ( 18) -( 21), summation is over all positive and negative shifts (including the term without shift, i.e., 0).
D is applied repeatedly.The generalization of (20) to N com- ponents should be obvious.

Example
The first two symmetries [11] of (2) are ), ( These symmetries are uniform in rank (rank 2 and 3, respectively).Symmetries of ranks 0 and 1 are both zero.

Example
The first two non-trivial symmetries [24] of (4), are uniform in rank.For example, 3 rank The symmetries of lower ranks are trivial.
An algorithm to compute polynomial generalized symmetries is described in detail in [24].

Definition
A This happens in most, but not all, cases.For  N component systems,  is an N N x matrix operator.
The defining equation for  [6,23] is where the bracket   , denotes the commutator of operators and  the composition of operators.The operator ) ( n u F was defined in (20).  F  is the Fréchet derivative of  in the direction of .F For the scalar case, the operator  is often of the form and in that case For the vector case and the examples under consideration, the elements of the N N x operator matrix  are of the form ).
Thus, for the two-component case [7]  

Algorithm for computation of recursion operators
We will now construct the recursion operator (32) for (4).In this case all the terms in (27) are 2 x 2 matrix operators.The construction uses the following steps: Step 1 (Determine the rank of the recursion The difference in rank of symmetries is used to compute the rank of the elements of the recursion operator.Using ( 7), ( 24) and (25), .4 Assuming that , ( to compute a rank matrix associated to the operator Step 2 (Determine the form of the recursion operator): The coefficients of these terms are admissible power combinations of 1 1 and , , , (which come from the terms on the right hand sides of ( 4)), so that all the terms have the correct rank.The maximum up-shift and down-shift operator that should be included can be determined by comparing two consecutive symmetries.Indeed, if the maximum up-shift in the first symmetry is , The same line of reasoning determines the minimum down-shift operator to be included.So, in our example As in the continuous case [10], 1  is a linear combination (with constant coeffi- cients jk c ~ of sums of all suitable products of symmetries and covariants (Fréchet derivatives of conserved densities) sandwiching .
where  denotes the matrix outer product, defined as Only the pair ) , ( G can be used, otherwise the ranks in (35) would be exceeded.Using (13) and ( 21), we compute Therefore, using (38) and renaming 10 c to , Adding (36) and (41), we obtain .
Step 3 (Determine the unknown coefficients): All the terms in (27) need to be computed.Referring to [7] for details, the result is: .

THE MATHEMATICA PACKAGE
To use the code, first load the Mathematica package DDERecursionOperator.m using the command ]; m" .onOperator DDERecursi [" Get := In [2] Proceeding with the KvM lattice (2) as an example, call the function DDERe-cursionOperator (which is part of the package) : The first part of the output (which we assign to R for later use) is indeed the recursion operator given in (31).Now using the first symmetry, generate the next symmetry by calling the function GenerateSymmetries (which is also part of the package): Evaluating the next five symmetries starting from the first one, can be done as follows: Due to the length of the output we do not show this result here.The Mathematica function TableForm will nicely reformat the output in a tabular form.Our package is available at [15].

Relativistic Toda (RT) Lattice
The RT lattice [20] is given as , ) ( ), ( and the recursion operator found by our package coincides with the one in [20]: . (51)

CONCLUDING REMARKS
The existence of a recursion operator is a corner stone in establishing the complete integrability of nonlinear DDEs because the recursion operators allows one to compute an infinite sequence of generalized symmetries.
Therefore, we presented an algorithm to compute recursion operators of nonlinear DDEs with polynomial terms.The algorithm uses the scaling properties, conservation laws, and generalized symmetries of the DDE, but does not require the knowledge of the bi-Hamiltonian operators.The algorithm has been implemented in Mathematica, a leading computer algebra system.The package DDERecursionOperators.m uses In-variantsSymmetries.m to compute the conservation laws and higher-order symmetries of nonlinear DDEs.
The algorithm presented in this paper works for many nonlinear DDEs, including the Kac-van Moerbeke (Volterra), modified Volterra, and Ablowitz-Ladik lattices, as well as standard and relativistic Toda lattices.However, the algorithm does not allow one to compute recursion operators for lattices due to Blaszak-Marciniak and Belov-Chaltikian (see, e.g., [20] for references).An extension of the algorithm that would cover these lattices is under investigation.