Symbolic Computation of Polynomial Conserved Densities,

Generalized Symmetries, and Recursion Operators for Nonlinear Differential-Difference Equations

Willy Hereman, Jan A. Sanders, Jack Sayers, and Jing Ping Wang

This paper is dedicated to Ryan Sayers (1982-2003)

Abstract. Algorithms for the symbolic computation of polynomial conserved densities, fluxes, generalized symmetries, and recursion operators for systems of nonlinear differential-difference equations are presented. In the algorithms we use discrete versions of the Fréchet and variational derivatives, as well as discrete Euler and homotopy operators.

The algorithms are illustrated for prototypical nonlinear polynomial lattices, including the Kac-van Moerbeke (Volterra) and Toda lattices. Results are shown for the modified Volterra and Ablowitz-Ladik lattices.

1. Introduction

The study of complete integrability of nonlinear differential-difference equations (DDEs) largely parallels that of PDEs [AC91, ASY00, W98]. Indeed, as in the continuous case, the existence of large numbers of generalized 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 symmetries of the equation into new symmetries. The existence of a recursion operator, which allows one to generate commuting symmetries, then confirms complete integrability.

There is a vast body of work on the complete integrability of DDEs. Consult e.g. [ASY00, HH03] for various approaches and additional references. In this article we describe algorithms to symbolically compute polynomial conservation laws, generalized symmetries, and recursion operators for DDEs. The design of these algorithms relies on related work for PDEs [GH97a, HG99, HGCM98] and work by Oevel et al [OZF93] and Zhang et al [ZTOF93].

1991 Mathematics Subject Classification. Primary 37J35, 37K10; Secondary 35Q58, 37K05.

Key words and phrases. Homotopy operators, conserved densities, fluxes, recursion operators, complete integrability, differential-difference, DDEs, semi-discrete lattices.

This material is based upon work supported by the National Science Foundation under Grant No. CCR-9901929.

In our algorithms we use discrete versions of the Fréchet and variational derivatives, as well as discrete Euler and homotopy operators. We draw on the close analogy between the continuous and discrete (in space) cases. As shown in [HM04, MH02], this analogy can be completely formalized and both theories can be formulated in terms of complexes. The same applies for the formulation in terms of Lie algebra complexes in [W98]. This allows one to translate by analogy the existing algorithms immediately (although complications arise when there is explicit dependence on the lattice point, n, in the discrete case). One of the more useful tools following from the abstract theory is the homotopy operator, which is based on scaling vector fields, and goes back to Poincaré in the continuous case. This operator allows one to directly integrate differential forms and can be straightforwardly implemented in computer algebra systems (CAS), since it reduces to integration over one scaling parameter. The discrete analogue as described in [HM04, MH02] does the corresponding job in the discrete case and we use it to compute fluxes. In this paper we do not explicitly use the abstract framework, yet, it has been a motivating force for the development of our algorithms.

In contrast to the symmetry approach in [ASY00], our algorithms rely on specific assumptions. For example, we will use the dilation invariance of DDEs in the construction of densities, symmetries, and recursion operators. At the cost of generality, our algorithms can be implemented in major CAS. Our Mathematica package InvariantsSymmetries.m [GH97b] computes densities and symmetries, and therefore aids in automated testing of complete integrability of semi-discrete lattices. Mathematica code to compute recursion operators of DDEs is under development.

The paper is organized as follows. In Section 2, we present key definitions, necessary tools, and prototypical examples which will be used throughout the paper. The examples include the Kac-van Moerbeke (KvM) [KvM75], Toda [T81] and Ablowitz-Ladik (AL) lattices [AL75]. The discrete higher Euler operators (variational derivatives) and the discrete homotopy operator are introduced in Section 3. These operators are applied in the construction of densities and fluxes in Section 5. The algorithm for higher-order symmetries is outlined in Section 6. The algorithms for scalar and matrix recursion operators are given in Section 7 and 8. The paper concludes with two more examples in Section 9.

2. Key Definitions and Prototypical Examples

Definition 2.1. A nonlinear (autonomous) DDE is an equation of the form

(2.1)  u̇_{n} = F(..., u_{n−1}, u_{n}, u_{n+1}, ...),

where u_{n} and F are vector-valued functions with N components. The integer n corresponds to discretization in space; the dot denotes one derivative with respect to the continuous time variable (t).

For simplicity, we denote the components of u_{n} by (u_{n}, v_{n}, w_{n}, ⋯) and write F(u_{n}), although F typically depends on u_{n} and a finite number of its forward and backward shifts. We assume that F is polynomial with constant coefficients. If present, parameters are denoted by lower-case Greek letters. No restrictions are imposed on the forward or backward shifts or the degree of nonlinearity in F.

Example 2.2. The AL lattice [AL75],

u̇_{n} = (u_{n+1} − 2u_{n} + u_{n−1}) + u_{n}v_{n}(u_{n+1} + u_{n−1}),
(2.2)  v̇_{n} = −(v_{n+1} − 2v_{n} + v_{n−1}) − u_{n}v_{n}(v_{n+1} + v_{n−1}),

DENSITIES, SYMMETRIES, AND RECURSION OPERATORS FOR NONLINEAR DDES 3

is a completely integrable discretization of the nonlinear Schrödinger equation.

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

Example 2.4. The KvM lattice [KvM75],

(2.3)  u̇_{n} = u_{n}(u_{n+1} − u_{n−1}),

is invariant under (t, u_{n}) → (λ^{−1}t, λu_{n}), where λ is an arbitrary scaling parameter.

Example 2.5. The Toda lattice [T81] in polynomial form [GH98],

(2.4)  u̇_{n} = v_{n−1} − v_{n}, v̇_{n} = v_{n}(u_{n} − u_{n+1}),

is invariant under the scaling symmetry

(2.5)  (t, u_{n}, v_{n}) → (λ^{−1}t, λ^{1}u_{n}, λ^{2}v_{n}).

Definition 2.6. We define the weight, w, of a variable as the exponent in the scaling parameter λ which multiplies the variable.

Since t is replaced by \frac{t}{λ}, we set w(\frac{d}{dt}) = w(D_{t}) = 1. In view of (2.5), we have w(u_{n}) = 1, and w(v_{n}) = 2 for the Toda lattice. Weights of dependent variables are nonnegative, rational, and independent of n. For example, w(u_{n−1}) = w(u_{n}) = w(u_{n+1}), etc.

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

Example 2.8. In the first equation of (2.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 (2.4) allows one to compute the weights of the dependent variables (and thus the scaling symmetry) with simple linear algebra. Indeed,

(2.6)  w(u_{n}) + 1 = w(v_{n}), w(v_{n}) + 1 = w(u_{n}) + w(v_{n}),

yields

(2.7)  w(u_{n}) = 1, w(v_{n}) = 2,

which is consistent with (2.5).

Dilation symmetries, which are special Lie-point symmetries, are common to many lattice equations. Lattices described by 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.

Example 2.9. The AL lattice is not dilation invariant. We therefore introduce an auxiliary parameter α. Replacing (2.2) by

u̇_{n} = α(u_{n+1} − 2u_{n} + u_{n−1}) + u_{n}v_{n}(u_{n+1} + u_{n−1}),
(2.8)  v̇_{n} = −α(v_{n+1} − 2v_{n} + v_{n−1}) − u_{n}v_{n}(v_{n+1} + v_{n−1}),

and requiring uniformity in rank, then gives

w(u_{n}) + 1 = w(α) + w(u_{n}) = 2w(u_{n}) + w(v_{n}),
(2.9)  w(v_{n}) + 1 = w(α) + w(v_{n}) = 2w(v_{n}) + w(u_{n}).

Obviously, w(u_{n}) + w(v_{n}) = w(α) = 1. Several scales are possible. The choice w(u_{n}) = w(v_{n}) = \frac{1}{2}, w(α) = 1, corresponds to the scaling symmetry

(2.10)  (t, u_{n}, v_{n}, α) → (λ^{−1}t, λ^{\frac{1}{2}} u_{n}, λ^{\frac{1}{2}} v_{n}, λα).

Definition 2.10. A scalar function ρ_{n}(u_{n}) is a conserved density of (2.1) if there exists a scalar function J_{n}(u_{n}) called the associated flux, such that [O93]

(2.11)  D_{t}ρ_{n} + Δ J_{n} = 0

is satisfied on the solutions of (2.1).

In (2.11), we used the (forward) difference operator, Δ = D − I, defined by

(2.12)  Δ J_{n} = (D − I) J_{n} = J_{n+1} − J_{n},

where D denotes the up-shift (forward or right-shift) operator, DJ_{n} = J_{n+1}, and I is the identity operator. 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 u_{n}.

Example 2.11. The first three density-flux pairs [HGCM98] for (2.3) are

(2.13) ρ_{n}^{(0)} = ln(u_{n}),  J_{n}^{(0)} = u_{n} + u_{n−1},
(2.14) ρ_{n}^{(1)} = u_{n},  J_{n}^{(1)} = u_{n}u_{n−1},
(2.15) ρ_{n}^{(2)} = \frac{1}{2}u_{n}^{2} + u_{n}u_{n+1},  J_{n}^{(2)} = u_{n−1}u_{n}(u_{n} + u_{n+1}).

Example 2.12. The first four density-flux pairs [GH98] for (2.4) are

(2.16) ρ_{n}^{(0)} = ln(v_{n}),  J_{n}^{(0)} = u_{n},
(2.17) ρ_{n}^{(1)} = u_{n},  J_{n}^{(1)} = v_{n−1},
(2.18) ρ_{n}^{(2)} = \frac{1}{2}u_{n}^{2} + v_{n},  J_{n}^{(2)} = u_{n}v_{n−1},
(2.19) ρ_{n}^{(3)} = \frac{1}{3}u_{n}^{3} + u_{n}(v_{n−1} + v_{n}), J_{n}^{(3)} = u_{n−1}u_{n}v_{n−1} + v_{n−1}^{2}.

The densities in (2.16)-(2.19) are uniform of ranks 0 through 3. The corresponding fluxes are also uniform in rank with ranks 1 through 4. In general, if in (2.11) rank ρ_{n} = R then rank J_{n} = R + 1, since w(D_{t}) = 1. All the pieces in (2.11) are uniform in rank. This comes as no surprise since the conservation law (2.11) holds on solutions of (2.1), hence it ‘inherits’ the dilation symmetry of (2.1).

In Section 5 we will give an algorithm to compute polynomial conserved densities and fluxes and use (2.4) to illustrate the steps. Non-polynomial densities (which are easy to find by hand) can be computed with the method given in [HH03].

Definition 2.13. Compositions of D and D^{−1} define an equivalence relation (≡) on monomials in the components of u_{n}. All shifted monomials are equivalent.

Example 2.14. For example, u_{n−1}v_{n+1} ≡ u_{n}v_{n+2} ≡ u_{n+1}v_{n+3} ≡ u_{n+2}v_{n+4}.

Factors in a monomial in u_{n} and its shifts are ordered by u_{n+j} ≺ u_{n+k} if j < k.

Definition 2.15. The main representative of an equivalence class of monomials in u_{n} and its shifts is the monomial with u_{n} in the first position.

Example 2.16. For example, u_{n}u_{n+2} (not u_{n−2}u_{n}) is the main representative of the class {⋯, u_{n−2}u_{n}, u_{n−1}u_{n+1}, u_{n}u_{n+2}, u_{n+1}u_{n+3}, ⋯ }.

For monomials involving u_{n}, v_{n}, w_{n}, etc. and their shifts, we lexicographically order the variables, i.e. u_{n} ≺ v_{n} ≺ w_{n}, etc.

Example 2.17. For example, u_{n}v_{n+2} (not u_{n−2}v_{n}) is the main representative of the class {⋯, u_{n−2}v_{n}, u_{n−1}v_{n+1}, u_{n}v_{n+2}, ⋯ }.

Definition 2.18. A vector function G(u_{n}) is called a generalized symmetry of (2.1) if the infinitesimal transformation u_{n} → u_{n} + ϵG leaves (2.1) invariant up to order ϵ. Consequently, G must satisfy [O93]

(2.20)  D_{t}G = F′(u_{n})[G]

on solutions of (2.1). F′(u_{n})[G] is the Fréchet derivative of F in the direction of G.

For the scalar case (N = 1), the Fréchet derivative in the direction of G is

(2.21)  F′(u_{n})[G] = \frac{∂}{∂ϵ}F (u_{n} + ϵG)|_{ϵ=0} = ∑_{k} \frac{∂F}{∂u_{n+k}} D^{k}G,

which defines the Fréchet derivative operator

(2.22)  F′(u_{n}) = ∑_{k} \frac{∂F}{∂u_{n+k}} D^{k}.

In the vector case with components u_{n} and v_{n}, the Fréchet derivative operator is

(2.23)  F′(u_{n}) = ( \begin{matrix}∑_{k} \frac{∂F_{1}}{∂u_{n+k}} D^{k} & ∑_{k} \frac{∂F_{1}}{∂v_{n+k}} D^{k} \\ ∑_{k} \frac{∂F_{2}}{∂u_{n+k}} D^{k} & ∑_{k} \frac{∂F_{2}}{∂v_{n+k}} D^{k} \end{matrix} ) .

Applied to G = (G_{1} G_{2})^{T}, where T is transpose, one obtains

(2.24)  F_{i}′(u_{n})[G] = ∑_{k} \frac{∂F_{i}}{∂u_{n+k}} D^{k}G_{1} + ∑_{k} \frac{∂F_{i}}{∂v_{n+k}} D^{k}G_{2}, i = 1, 2.

In (2.21) and (2.24) summation is over all positive and negative shifts (including k = 0). For k > 0, D^{k} = D∘D∘⋯∘D (k times). Similarly, for k < 0 the down-shift operator D^{−1} is applied repeatedly. The generalization of (2.23) to a system with N components should be obvious.

Example 2.19. The first two symmetries [HGCM98] of (2.3) are

(2.25) G^{(1)} = u_{n}(u_{n+1} − u_{n−1}),
(2.26) G^{(2)} = u_{n}u_{n+1}(u_{n} + u_{n+1} + u_{n+2}) − u_{n−1}u_{n}(u_{n−2} + u_{n−1} + u_{n}).

These symmetries are uniform in rank (rank 2 and 3, respectively). The symmetries of ranks 0 and 1 are both zero.

Example 2.20. The first two non-trivial symmetries [HGCM98] of (2.4) are

(2.27)  G^{(1)} = ( \begin{matrix}v_{n−1} − v_{n} \\ v_{n}(u_{n} − u_{n+1}) \end{matrix} ) ,
(2.28)  G^{(2)} = ( \begin{matrix}v_{n}(u_{n} + u_{n+1}) − v_{n−1}(u_{n−1} + u_{n}) \\ v_{n}(u_{n+1}^{2} − u_{n}^{2} + v_{n+1} − v_{n−1}) \end{matrix} ) .

The above symmetries are uniform in rank. For example, rank G_{1}^{(2)} = 3 and rank G_{2}^{(2)} = 4. Symmetries of lower ranks are trivial.

An algorithm to compute polynomial generalized symmetries will be outlined in Section 6 and applied to (2.4).

Definition 2.21. A recursion operator R connects symmetries

(2.29)  G^{(j+s)} = RG^{(j)},

where j = 1, 2, ..., and s is the gap length. The symmetries are linked consecutively if s = 1. This happens in most, but not all, cases. For N-component systems, R is an N × N matrix operator.

The defining equation for R [O93, W98] is

(2.30)  D_{t}R + [R, F′(u_{n})] = \frac{∂R}{∂t} + R′[F] + R ∘ F′(u_{n}) − F′(u_{n}) ∘ R = 0,

where [ , ] denotes the commutator and ∘ the composition of operators. The operator F′(u_{n}) was defined in (2.23). R′[F] is the Fréchet derivative of R in the direction of F. For the scalar case, the operator R is often of the form

(2.31)  R = U(u_{n}) O((D − I)^{−1}, D^{−1}, I, D) V (u_{n}),

and then

(2.32)  R′[F] = ∑_{k} (D^{k}F) \frac{∂U}{∂u_{n+k}} O V + ∑_{k} UO (D^{k}F) \frac{∂V}{∂u_{n+k}} .

For the vector case, the elements of the N × N operator matrix R are often of the form R_{ij} = U_{ij}(u_{n}) O_{ij}((D − I)^{−1}, D^{−1}, I, D) V_{ij}(u_{n}). For the two-component case

R′[F]_{ij} = ∑_{k} (D^{k}F_{1}) \frac{∂U_{ij}}{∂u_{n+k}} O_{ij} V_{ij} + ∑_{k} (D^{k}F_{2}) \frac{∂U_{ij}}{∂v_{n+k}} O_{ij} V_{ij}

+∑_{k} U_{ij}O_{ij} (D^{k}F_{1}) \frac{∂V_{ij}}{∂u_{n+k}} + ∑_{k} U_{ij}O_{ij} (D^{k}F_{2}) \frac{∂V_{ij}}{∂v_{n+k}} . (2.33)

Example 2.22. The KvM lattice (2.3) has recursion operator

R = u_{n}(I + D)(u_{n}D − D^{−1}u_{n})(D − I)^{−1} \frac{1}{u_{n}} I
(2.34) = u_{n}D^{−1} + (u_{n} + u_{n+1})I + u_{n}D + u_{n}(u_{n+1} − u_{n−1})(D − I)^{−1} \frac{1}{u_{n}} I.

Example 2.23. The Toda lattice (2.4) has recursion operator

(2.35)  R = ( \begin{matrix}−u_{n}I & −D^{−1} − I + (v_{n−1} − v_{n})(D − I)^{−1} \frac{1}{v_{n}} I \\ −v_{n}I − v_{n}D & −u_{n+1}I + v_{n}(u_{n} − u_{n+1})(D − I)^{−1} \frac{1}{v_{n}} I \end{matrix} ) .

In Section 7 we will give an algorithm for the computation of scalar recursion operators like (2.34). In Section 8 we cover the matrix case and show how (2.35) is computed. The algorithms complement those for recursion operators of PDEs presented in [HG99] and elsewhere in these proceedings [BHS04]. We now introduce two powerful tools which will be used in the computation of densities and fluxes.

3. The Discrete Variational Derivative (Euler Operator)

Definition 3.1. A function E_{n}(u_{n}) is a total difference if there exists a function J_{n}(u_{n}) such that E_{n} = Δ J_{n} = (D − I) J_{n}.

Theorem 3.2. A necessary and sufficient condition for a function E_{n}, with positive shifts up to level p_{0}, to be a total difference is that

(3.1)  L_{u_{n}}^{(0)}(E_{n}) = 0.

L_{u_{n}}^{(0)} is the discrete variational derivative (Euler operator) [ASY00] defined by

(3.2)  L_{u_{n}}^{(0)} = \frac{∂}{∂u_{n}} (∑_{k=0}^{p_{0}} D^{−k}) = \frac{∂}{∂u_{n}} (I + D^{−1} + D^{−2} + ⋯ + D^{−p_{0}}).

A proof is given in e.g. [HH03].

Remark 3.3. To verify that an expression E(u_{n−q}, ⋯, u_{n}, ⋯, u_{n+p}) involving negative shifts is a total difference, one must first remove the negative shifts by replacing E_{n} by Ẽ_{n} = D^{q}E_{n}. Applied to Ẽ_{n}, (3.2) terminates at p_{0} = p + q.

We now introduce a tool to invert the total difference operator Δ = D − I.

4. The Discrete Homotopy Operator

We assume that L_{u_{n}}^{(0)}(E_{n}) = 0. So, ∃J_{n} such that E_{n} = ΔJ_{n}. To compute J_{n} = Δ^{−1}(E_{n}) one must invert the operator Δ = D − I. Working with the formal inverse, Δ^{−1} = D^{−1} + D^{−2} + D^{−3} + …, is impractical, perhaps impossible. We therefore present the (discrete) homotopy operator which circumvents the formal infinite series. In analogy to the continuous case [O93], we first introduce the discrete higher Euler operators.

Definition 4.1. The discrete higher Euler operators are defined by

(4.1)  L_{u_{n}}^{(i)} = \frac{∂}{∂u_{n}} (∑_{k=i}^{p_{i}} \binom{k}{i}D^{−k}), i ≥ 1.

The higher Euler operators all terminate at some maximal shifts p_{i}.

Example 4.2. For scalar component u_{n}, the first higher Euler operators are:

(4.2) L_{u_{n}}^{(1)} = \frac{∂}{∂u_{n}} (D^{−1} + 2D^{−2} + 3D^{−3} + 4D^{−4} + ⋯ + p_{1} D^{−p_{1}}),
(4.3) L_{u_{n}}^{(2)} = \frac{∂}{∂u_{n}} (D^{−2} + 3D^{−3} + 6D^{−4} + 10D^{−5} + ⋯ + \frac{1}{2}p_{2}(p_{2} − 1) D^{−p_{2}}),

where the p_{i} are the highest shifts in the arguments of the operators.

Similar formulae hold for L_{v_{n}}^{(i)}. Next, we introduce the homotopy operators. For simplicity we show the formulae for the two component case u_{n} = (u_{n}, v_{n}).

Definition 4.3. The homotopy operator is defined as,

(4.4)  H_{u_{n}} = ∫_{0}^{1}(f_{1,n}(u_{n})[λu_{n}] + f_{2,n}(u_{n})[λu_{n}])\frac{dλ}{λ},

where

(4.5) f_{1,n}(u_{n}) = ∑_{i=0}^{m_{1}−1}(D − I)^{i}(u_{n}L_{u_{n}}^{(i+1)}), f_{2,n}(u_{n}) = ∑_{i=0}^{m_{2}−1}(D − I)^{i}(v_{n}L_{v_{n}}^{(i+1)}),

where m_{1} and m_{2} are the highest shifts in u_{n} and v_{n} in the arguments of the operators. Furthermore, f_{1,n}(u_{n})[λu_{n}] means that in f_{1,n}(u_{n}) one replaces u_{n} → λu_{n}, and, of course, u_{n+1} → λu_{n+1}, u_{n+2} → λu_{n+2}, etc. Same for f_{2,n}(u_{n})[λu_{n}].

Theorem 4.4. J_{n} = Δ^{−1}(E_{n}) can be computed as J_{n} = H_{u_{n}}(E_{n}).

A similar theorem and proof for continuous homotopy operators is given in [O93, p. 372]. By constructing an appropriate discrete variational bicomplex one can prove the above theorem. See [HM04, MH02] and elsewhere in these proceedings [MQ04].

5. Computation of Polynomial Conserved Densities and Fluxes

As an example, we compute density-flux pair (2.19) for (2.4). Assuming that the weights (2.7) are computed and the rank of the density is selected (R = 3 here), the algorithm has three steps. See [GH98, GHE97, HGCM98] for details about the algorithm and alternate ways to compute fluxes. Step 1: Construct the form of the density List all monomials in u_{n} and v_{n} of rank 3 or less: M={u_{n}^{3}, u_{n}^{2}, u_{n}v_{n}, u_{n}, v_{n}}.

Next, for each monomial in M, introduce the correct number of t-derivatives so that each term has weight 3. Thus, using (2.4),

\frac{d^{0}}{dt^{0}} (u_{n}^{3}) = u_{n}^{3}, \frac{d^{0}}{dt^{0}} (u_{n}v_{n}) = u_{n}v_{n},

\frac{d}{dt}(u_{n}^{2}) = 2u_{n}u̇_{n} = 2u_{n}v_{n−1} − 2u_{n}v_{n}, \frac{d}{dt}(v_{n}) = v̇_{n} = u_{n}v_{n} − u_{n+1}v_{n},

\frac{d^{2}}{dt^{2}} (u_{n}) = \frac{d}{dt}(u̇_{n}) = \frac{d}{dt}(v_{n−1} − v_{n}) = u_{n−1}v_{n−1} − u_{n}v_{n−1} − u_{n}v_{n} + u_{n+1}v_{n}.

Gather all terms in N = {u_{n}^{3}, u_{n}v_{n−1}, u_{n}v_{n}, u_{n−1}v_{n−1}, u_{n+1}v_{n}}. Identify members belonging to the same equivalence classes and replace them by their main representatives. For example, u_{n}v_{n−1} ≡ u_{n+1}v_{n}, so the latter is replaced by u_{n}v_{n−1}. So, N is replaced by P = {u_{n}^{3}, u_{n}v_{n−1}, u_{n}v_{n}}, which has the building blocks of the density. Linear combination of the monomials in P with constant coefficients c_{i} gives the candidate density,

(5.1)  ρ_{n} = c_{1} u_{n}^{3} + c_{2} u_{n}v_{n−1} + c_{3} u_{n}v_{n}.

Step 2: Determine the coefficients Require that (2.11) holds. Compute D_{t}ρ_{n}. Use (2.4) to remove u̇_{n} and v̇_{n} and their shifts. Thus,

E_{n} = D_{t}ρ_{n} = (3c_{1} − c_{2})u_{n}^{2}v_{n−1} + (c_{3} − 3c_{1})u_{n}^{2}v_{n} + (c_{3} − c_{2})v_{n−1}v_{n}

+ c_{2} u_{n−1}u_{n}v_{n−1} + c_{2}v_{n−1}^{2} − c_{3}u_{n}u_{n+1}v_{n} − c_{3}v_{n}^{2}. (5.2)

Compute Ẽ_{n} = DE_{n}. First, apply (3.2) for component u_{n} to Ẽ_{n} :

L_{u_{n}}^{(0)}(Ẽ_{n}) = \frac{∂}{∂u_{n}} (I + D^{−1} + D^{−2})(Ẽ_{n})
= 2(3c_{1} − c_{2})u_{n}v_{n−1} + 2(c_{3} − 3c_{1})u_{n}v_{n}

+ (c_{2} − c_{3})u_{n−1}v_{n−1} + (c_{2} − c_{3})u_{n+1}v_{n}. (5.3)

Second, apply (3.2) for component v_{n} to Ẽ_{n} :

L_{v_{n}}^{(0)}(Ẽ_{n}) = \frac{∂}{∂v_{n}} (I + D^{−1})(Ẽ_{n})
= (3c_{1} − c_{2})u_{n+1}^{2} + (c_{3} − c_{2})v_{n+1} + (c_{2} − c_{3})u_{n}u_{n+1}

+ 2(c_{2} − c_{3})v_{n} + (c_{3} − 3c_{1})u_{n}^{2} + (c_{3} − c_{2})v_{n−1}. (5.4)

Both (5.3) and (5.4) must vanish identically. Solve the linear system

(5.5)  S = {3c_{1} − c_{2} = 0, c_{3} − 3c_{1} = 0, c_{2} − c_{3} = 0}.

The solution is 3c_{1} = c_{2} = c_{3}. Substituting c_{1} = \frac{1}{3}, c_{2} = c_{3} = 1, into (5.1)

(5.6)  ρ_{n} = \frac{1}{3}u_{n}^{3} + u_{n}(v_{n−1} + v_{n}).

Step 3: Compute the flux In view of (2.11), we will compute J_{n} = Δ^{−1}(−E_{n}) with the homotopy operator introduced in Section 4. Alternative methods are described in [GH98, HH03].

Insert c_{1} = \frac{1}{3}, c_{2} = c_{3} = 1 into (5.2) and compute

(5.7)  Ẽ_{n} = DE_{n} = u_{n}u_{n+1}v_{n} + v_{n}^{2} − u_{n+1}u_{n+2}v_{n+1} − v_{n+1}^{2}.

We apply formulae (4.5) to −Ẽ_{n}. The pieces are listed in Tables 1 and 2.

\begin{matrix}i & L_{u_{n}}^{(i+1)}(−Ẽ_{n}) & (D − I)^{i}(u_{n}L_{u_{n}}^{(i+1)}(−Ẽ_{n})) \\ 0 & u_{n−1}v_{n−1}+u_{n+1}v_{n} & u_{n}u_{n−1}v_{n−1}+u_{n}u_{n+1}v_{n} \\ 1 & u_{n−1}v_{n−1} & u_{n+1}u_{n}v_{n}−u_{n}u_{n−1}v_{n−1} \end{matrix}

Table 1. Computation of f_{1,n}(u_{n})(−Ẽ_{n})

\begin{matrix}i & L_{v_{n}}^{(i+1)}(−Ẽ_{n}) & (D − I)^{i}(v_{n}L_{v_{n}}^{(i+1)}(−Ẽ_{n})) \\ 0 & u_{n}u_{n+1}+2v_{n} & v_{n}u_{n}u_{n+1}+2v_{n}^{2} \end{matrix}

Table 2. Computation of f_{2,n}(u_{n})(−Ẽ_{n})

Adding the terms in the right columns in Table 1 and Table 2,

(5.8)  f_{1,n}(u_{n})(−Ẽ_{n}) = 2u_{n}u_{n+1}v_{n}, f_{2,n}(u_{n})(−Ẽ_{n}) = u_{n}u_{n+1}v_{n} + 2v_{n}^{2}.

Thus, the homotopy operator (4.4) gives

J̃_{n} = ∫_{0}^{1}(f_{1,n}(u_{n})(−Ẽ_{n})[λu_{n}] + f_{2,n}(u_{n})(−Ẽ_{n})[λu_{n}]) \frac{dλ}{λ}
= ∫_{0}^{1}(3λ^{2}u_{n}u_{n+1}v_{n} + 2λv_{n}^{2}) dλ
(5.9)  = u_{n}u_{n+1}v_{n} + v_{n}^{2}.

After a backward shift, J_{n} = D^{−1}J̃_{n}, we obtain the final result:

(5.10)  ρ_{n} = \frac{1}{3} u_{n}^{3} + u_{n}(v_{n−1} + v_{n}), J_{n} = u_{n−1}u_{n}v_{n−1} + v_{n−1}^{2}.

6. Computation of Polynomial Generalized Symmetries As an example, we compute symmetry (2.28) of rank (3, 4) for (2.4). The two

steps of the algorithm [GH99] are similar to those in Section 5. Step 1: Construct the form of the symmetry Start by listing all monomials in u_{n} and v_{n} of ranks 3 and 4, or less:

M_{1} = {u_{n}^{3}, u_{n}^{2}, u_{n}v_{n}, u_{n}, v_{n}}, M_{2} = {u_{n}^{4}, u_{n}^{3}, u_{n}^{2}v_{n}, u_{n}^{2}, u_{n}v_{n}, u_{n}, v_{n}^{2}, v_{n}}.

As in Step 1 in Section 5, for each monomial in M_{1} and M_{2}, introduce the necessary t-derivatives so that each term has rank 3 and 4, respectively. At the same time, use (2.4) to remove all t−derivatives. Doing so, based on M_{1} we obtain

(6.1)  N_{1} = {u_{n}^{3}, u_{n−1}v_{n−1}, u_{n}v_{n−1}, u_{n}v_{n}, u_{n+1}v_{n}}.

Similarly, based on M_{2}, we get

N_{2} = {u_{n}^{4}, u_{n−1}^{2}v_{n−1}, u_{n−1}u_{n}v_{n−1}, u_{n}^{2}v_{n−1}, v_{n−2}v_{n−1}, v_{n−1}^{2}, u_{n}^{2}v_{n},

u_{n}u_{n+1}v_{n}, u_{n+1}^{2}v_{n}, v_{n−1}v_{n}, v_{n}^{2}, v_{n}v_{n+1}}. (6.2)

In contrast to the strategy for densities, we do not introduce the main representatives, but linearly combine the monomials in N_{1} and N_{2} with constant coefficients c_{i} to get the form of the candidate symmetry:

G_{1} = c_{1} u_{n}^{3} + c_{2} u_{n−1}v_{n−1} + c_{3} u_{n}v_{n−1} + c_{4} u_{n}v_{n} + c_{5} u_{n+1}v_{n},
G_{2} = c_{6} u_{n}^{4} + c_{7} u_{n−1}^{2}v_{n−1} + c_{8} u_{n−1}u_{n}v_{n−1} + c_{9} u_{n}^{2}v_{n−1} + c_{10} v_{n−2}v_{n−1}
+ c_{11} v_{n−1}^{2} + c_{12} u_{n}^{2}v_{n} + c_{13} u_{n}u_{n+1}v_{n} + c_{14} u_{n+1}^{2}v_{n} + c_{15} v_{n−1}v_{n}

+ c_{16} v_{n}^{2} + c_{17} v_{n}v_{n+1}, (6.3)

Step 2: Determine the coefficients

To compute the coefficients c_{i} we require that (2.20) holds on solutions of (2.1). Compute D_{t}G_{1} and D_{t}G_{2} and use (2.4) to remove u̇_{n}, v̇_{n}, and their shifts. This gives the left hand sides of (2.20). Use (2.24) to get the right hand sides of (2.20):

F_{1}′(u_{n})[G] = D^{−1}G_{2} − IG_{2},
(6.4)  F_{2}′(u_{n})[G] = v_{n}IG_{1} − v_{n}DG_{1} + (u_{n} − u_{n+1})IG_{2}.

Substitute (6.3) into (6.4) and equate the corresponding left and right hand sides. Since all monomials in u_{n}, v_{n} and their shifts are independent, one obtains the linear system that determines the coefficients c_{i}. The solution is

c_{1} = c_{6} = c_{7} = c_{8} = c_{9} = c_{10} = c_{11} = c_{13} = c_{16} = 0,

c_{2} = c_{3} = −c_{4} = −c_{5} = c_{12} = −c_{14} = c_{15} = −c_{17}. (6.5)

With the choice c_{17} = 1, symmetry (6.3) finally becomes

G_{1} = v_{n}(u_{n} + u_{n+1}) − v_{n−1}(u_{n−1} − u_{n}),
(6.6)  G_{2} = v_{n}(u_{n+1}^{2} − u_{n}^{2} + v_{n+1} − v_{n−1}).

7. Computation of Scalar Recursion Operators

Recursion operators can be constructed as R = H_{2}H_{1}^{−1}, where H_{1} and H_{2} are appropriate Hamiltonian operators [OZF93]. We use a less general but direct method based on dilation invariance. Furthermore, we assume that the form of discrete recursion operators mimics its continuous counterpart [B93, M00, SW01a, SW01b, W98]. With the appropriate ansatz for R, defining equation (2.30) is then used to determine the coefficients.

In this section we show how to compute recursion operator (2.34) of (2.4). We observe that (2.34) in expanded form naturally splits into two pieces:

(7.1)  R = R_{0} + R_{1},

where R_{0} contains only terms with shift operators D^{−1}, I, and D, and R_{1} has terms involving (D − I)^{−1}.

Step 1: Determine the rank of the recursion operator In view of (2.29) and assuming that the symmetries are linked consecutively (s = 1), the recursion operator R has rank

(7.2)  R = rank R = rank G^{(2)} − rank G^{(1)} = 3 − 2 = 1,

where we used (2.25) and (2.26). If the assumption turns out to be correct, then the recursion operator has rank 1. If the assumption were false because symmetries are not linked consecutively, then R must be adjusted and the subsequent steps must be repeated. See [BHS04] for examples of PDEs for which that happens. Step 2: Determine the form of the recursion operator We split this into two sub-steps.

(i) Determine the form of R_{0} The candidate R_{0} is a sum of terms involving D^{−1}, I, and D. The coefficients of these terms are linear combinations of u_{n−1}, u_{n}, and u_{n+1}, so that all terms have the correct rank.

R_{0} = (c_{1}u_{n−1} + c_{2}u_{n} + c_{3}u_{n+1}) D^{−1} + (c_{4}u_{n−1} + c_{5}u_{n} + c_{6}u_{n+1}) I

+ (c_{7}u_{n−1} + c_{8}u_{n} + c_{9}u_{n+1}) D, (7.3)

where the c_{i} are constant coefficients.

A few remarks are in place. First, in R_{0} we moved the operators D^{−1}, I, and D all the way to the right. Second, 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 u_{n+p}, and the maximum up-shift in the next symmetry is u_{n+p+r}, then R_{0} must have D, D^{2}, ⋯, D^{r}. The same line of reasoning determines the minimum down-shift operator to be included. In our example, there is no need to include terms in D^{−2}, D^{2}, etc. Third, the coefficients of the operators can be restricted to linear combinations of the terms appearing in F. Hence, no terms in u_{n±2}, u_{n±3} and so on occur in (7.3). (ii) Determine the form of R_{1} As in the continuous case [HG99], R_{1} is a linear combination (with constant coefficients c̃_{jk}) of sums of all suitable products of symmetries and covariants (Fréchet derivatives of densities) sandwiching (D − I)^{−1}, i.e.

(7.4)  R_{1} = ∑_{j} ∑_{k} c̃_{jk}G^{(j)}(D − I)^{−1}ρ_{n}^{(k)′}.

For (2.3), G^{(1)} in (2.25) and ρ_{n}^{(0)} = ln(u_{n}) in (2.13) are the only suitable pair. Indeed, using (2.22) we have ρ_{n}^{(0)′} = ln(u_{n})′ = \frac{1}{u_{n}} I, which has rank -1. Combined with G^{(1)} of rank 2, we have a term of rank 1. Other combinations of symmetries and covariants would exceed rank 1. Therefore,

(7.5)  R_{1} = c̃_{10}u_{n}(u_{n+1} − u_{n−1})(D − I)^{−1}(\frac{1}{u_{n}})I

where c̃_{10} a constant coefficient. Using (7.1) and renaming c̃_{10} to c_{10},

R = (c_{1}u_{n−1} + c_{2}u_{n} + c_{3}u_{n+1})D^{−1} + (c_{4}u_{n−1} + c_{5}u_{n} + c_{6}u_{n+1})I

+ (c_{7}u_{n−1} + c_{8}u_{n} + c_{9}u_{n+1})D + c_{10}u_{n}(u_{n+1} − u_{n−1})(D − I)^{−1}(\frac{1}{u_{n}})I (7.6)

is the candidate recursion operator for (2.3).

Step 3: Determine the coefficients Starting from (7.6), we use (2.32) with F = u_{n}(u_{n+1} − u_{n−1}) to compute R′[F]. The partial results will not be shown due to length. Using (2.21), we compute

(7.7)  F′ = −u_{n}D^{−1} + (u_{n+1} − u_{n−1})I + u_{n}D.

Then we compute R ∘ F′ and F′ ∘ R. After substituting the pieces into (2.30) we simplify the resulting expression using rules such as

(D − I)^{−1}D = D(D − I)^{−1} = I + (D − I)^{−1},
(7.8)  (D − I)^{−1}D^{−1} = D^{−1}(D − I)^{−1} = −D^{−1} + (D − I)^{−1}.

We further simplify by recursively using formulas like

DU(u_{n})(D−I)^{−1}V (u_{n})I = U(u_{n+1})V (u_{n})I+U(u_{n+1})(D−I)^{−1}V (u_{n})I,
(7.9) (D−I)^{−1}U(u_{n})V (u_{n})D = U(u_{n−1})V (u_{n−1})I+(D−I)^{−1}U(u_{n−1})V (u_{n−1})I.

Finally, we equate like terms to obtain a linear system for the c_{i}. Substituting the solution

(7.10)  c_{1} = c_{3} = c_{4} = c_{7} = c_{9} = 0, c_{2} = c_{5} = c_{6} = c_{8} = c_{10} = 1

into (7.6) we obtain the final result

(7.11) R = u_{n}D^{−1} + (u_{n} + u_{n+1})I + u_{n}D + u_{n}(u_{n+1} − u_{n−1})(D − I)^{−1} \frac{1}{u_{n}} I.

A straightforward computation confirms that RG^{(1)} = G^{(2)} with G^{(1)} in (2.25) and G^{(2)} in (2.26).

8. Computation of Matrix Recursion Operators

We construct the recursion operator (2.35) for (2.4). Now all the terms in (2.30) are 2 × 2 matrix operators. Step 1: Determine the rank of the recursion operator The difference in rank of symmetries is again used to compute the rank of the elements of the recursion operator. Using (2.7), (2.27) and (2.28),

(8.1)  rank G^{(1)} = ( \begin{matrix}2 \\ 3 \end{matrix} ) , rank G^{(2)} = ( \begin{matrix}3 \\ 4 \end{matrix} ) .

Assuming that R G^{(1)} = G^{(2)}, we use the formula

(8.2)  rank R_{ij} = rank G_{i}^{(k+1)} − rank G_{j}^{(k)},

to compute a rank matrix associated to the operator

(8.3)  rank R = ( \begin{matrix}1 & 0 \\ 2 & 1 \end{matrix} ) .

Step 2: Determine the form of the recursion operator As in the scalar case, we build a candidate R_{0} :

(8.4)  R_{0} = ( \begin{matrix}(R_{0})_{11} & (R_{0})_{12} \\ (R_{0})_{21} & (R_{0})_{22} \end{matrix} ) ,

with (R_{0})_{11} = (c_{1}u_{n} + c_{2}u_{n+1}) I,
(R_{0})_{12} = c_{3}D^{−1} + c_{4}I,
(R_{0})_{21} = (c_{5}u_{n}^{2} + c_{6}u_{n}u_{n+1} + c_{7}u_{n+1}^{2} + c_{8}v_{n−1} + c_{9}v_{n}) I
+ (c_{10}u_{n}^{2} + c_{11}u_{n}u_{n+1} + c_{12}u_{n+1}^{2} + c_{13}v_{n−1} + c_{14}v_{n}) D,
(R_{0})_{22} = (c_{15}u_{n} + c_{16}u_{n+1}) I.

Analogous to the scalar case, the elements of matrix R_{1} are linear combinations with constant coefficients c̃_{jk}) of all suitable products of symmetries and covariants sandwiching (D − I)^{−1}. Hence,

(8.5)  ∑_{j} ∑_{k} c̃_{jk}G^{(j)}(D − I)^{−1} ⊗ ρ_{n}^{(k)′},

where ⊗ denotes the matrix outer product, defined as

( \begin{matrix}G_{1}^{(j)} \\ G_{2}^{(j)} \end{matrix} ) (D−I)^{−1} ⊗ (ρ_{n,1}^{(k)′} ρ_{n,2}^{(k)′}) = ( \begin{matrix}G_{1}^{(j)}(D−I)^{−1}ρ_{n,1}^{(k)′} & G_{1}^{(j)}(D−I)^{−1}ρ_{n,2}^{(k)′} \\ G_{2}^{(j)}(D−I)^{−1}ρ_{n,1}^{(k)′} & G_{2}^{(j)}(D−I)^{−1}ρ_{n,2}^{(k)′} \end{matrix} ) .

Only the pair (G^{(1)}, ρ_{n}^{(0)′}) can be used, otherwise the ranks in (8.3) would be exceeded. Using (2.24) and (2.16) we compute

(8.6)  ρ_{n}^{(0)′} = ( 0 \frac{1}{v_{n}} I ) ,

Therefore, using (8.5) and renaming c̃_{10} to c_{17},

(8.7)  R_{1} = ( \begin{matrix}0 & c_{17}(v_{n−1} − v_{n})(D−I)^{−1} \frac{1}{v_{n}} I \\ 0 & c_{17}v_{n}(u_{n} − u_{n+1})(D−I)^{−1} \frac{1}{v_{n}} I \end{matrix} ) .

Adding (8.4) and (8.7) we obtain

(8.8)  R = ( \begin{matrix}R_{11} & R_{12} \\ R_{21} & R_{22} \end{matrix} ) ,

with

R_{11} = (c_{1}u_{n} + c_{2}u_{n+1}) I
R_{12} = c_{3}D^{−1} + c_{4}I + c_{17}(v_{n−1} − v_{n})(D − I)^{−1} \frac{1}{v_{n}} I,
R_{21} = (c_{5}u_{n}^{2} + c_{6}u_{n}u_{n+1} + c_{7}u_{n+1}^{2} + c_{8}v_{n−1} + c_{9}v_{n}) I
+ (c_{10}u_{n}^{2} + c_{11}u_{n}u_{n+1} + c_{12}u_{n+1}^{2} + c_{13}v_{n−1} + c_{14}v_{n}) D,
R_{22} = (c_{15}u_{n} + c_{16}u_{n+1})I + c_{17}v_{n}(u_{n} − u_{n+1})(D − I)^{−1} \frac{1}{v_{n}} I.

Step 3: Determine the coefficients All the terms in (2.30) need to be computed. The strategy is similar to the scalar case, yet the computations are much more cumbersome. Omitting the details, the result is: c_{2} = c_{5} = c_{6} = c_{7} = c_{8} = c_{10} = c_{11} = c_{12} = c_{13} = c_{15} = 0, c_{1} = c_{3} = c_{4} = c_{9} = c_{14} = c_{16} = −1, and c_{17} = 1. Substitution the constants into (8.8) gives

(8.9)  R = ( \begin{matrix}−u_{n}I & −D^{−1} − I + (v_{n−1} − v_{n})(D − I)^{−1} \frac{1}{v_{n}} I \\ −v_{n}I − v_{n}D & −u_{n+1}I + v_{n}(u_{n} − u_{n+1})(D − I)^{−1} \frac{1}{v_{n}} I \end{matrix} ) .

It is straightforward to verify that RG^{(1)} = G^{(2)} with G^{(1)} in (2.27) and G^{(2)} in (2.28).

9. More Examples

Example 9.1. The modified Volterra lattice [ASY00, HH03],

(9.1)  u̇_{n} = u_{n}^{2}(u_{n+1} − u_{n−1}),

has two non-polynomial densities ρ_{n}^{(0)} = \frac{1}{u_{n}} and ρ_{n}^{(1)} = ln(u_{n}), and infinitely many polynomial densities. The first two symmetries,

(9.2)  G^{(1)} = u_{n}^{2}(u_{n+1} − u_{n−1}),
(9.3)  G^{(2)} = u_{n}^{2}u_{n+1}^{2}(u_{n} + u_{n+2}) − u_{n−1}^{2}u_{n}^{2}(u_{n−2} + u_{n}),

are linked by the recursion operator

(9.4)  R = u_{n}^{2}D^{−1} + 2u_{n}u_{n+1}I + u_{n}^{2}D + 2u_{n}^{2}(u_{n+1} − u_{n−1})(D − I)^{−1} \frac{1}{u_{n}} I.

Example 9.2. The AL lattice (2.2) has infinitely many densities [GH98] and symmetries [GH99]. The recursion operator is of the form (8.8) with

R_{11} = P_{n}D^{−1} − u_{n}Δ^{−1}v_{n+1}I − u_{n−1}P_{n}Δ^{−1} \frac{v_{n}}{P_{n}} I,
R_{12} = −u_{n}u_{n−1}I − u_{n}Δ^{−1}u_{n−1}I − u_{n−1}P_{n}Δ^{−1} \frac{u_{n}}{P_{n}} I,
R_{21} = v_{n}v_{n+1}I + v_{n}Δ^{−1}v_{n+1}I + v_{n+1}P_{n}Δ^{−1} \frac{v_{n}}{P_{n}} I,
(9.5) R_{22} = (u_{n}v_{n+1} + u_{n−1}v_{n})I + P_{n}D + v_{n}Δ^{−1}v_{n−1}I + v_{n+1}P_{n}Δ^{−1} \frac{u_{n}}{P_{n}} I,

where P_{n} = 1 + u_{n}v_{n} and Δ = D − I.

This recursion operator has an inverse which is quite exceptional. For further information see e.g. [ZTOF93].

Acknowledgements M. Hickman and B. Deconinck are gratefully acknowledged for valuable discussions. D. Baldwin is thanked for proof reading the manuscript.

References

[AC91] M. J. Ablowitz and P. A. Clarkson, Solitons, nonlinear evolution equations and inverse scattering, London Math. Soc. Lect. Note Ser. 149, Cambridge University Press, Cambridge, U. K., 1991.

[AL75] M. J. Ablowitz and J. F. Ladik, Nonlinear differential-difference equations, J. Math. Phys. 16 (1975), 598–603.

[ASY00] V. E. Adler, A. B. Shabat and R. I. Yamilov, Symmetry approach to the integrability problem, Theor. Math. Phys. 125 (2000), 1603–1661.

[BHS04] D. Baldwin, W. Hereman and J. Sayers, Symbolic algorithms for the Painlevé test, special solutions and recursion operators for nonlinear PDEs, present proceedings.

[B93] A. H. Bilge, On the equivalence of linearization and formal symmetries as integrability tests for evolution equations, J. Phys. A: Math. Gen. 26 (1993), 7511–7519.

[GH97a] Ü. Göktaş and W. Hereman, Symbolic computation of conserved densities for systems of nonlinear evolution equations, J. Symb. Comp. 24 (1997), 591–621.

[GH97b] , Mathematica package InvariantsSymmetries.m is available since 1997 from http://www.mathsource.com/cgi-bin/msitem?0208-932.

[GH98] , Computation of conserved densities for nonlinear lattices, Physica D 123 (1998), 425–436.

[GH99] , Algorithmic computation of higher-order symmetries for nonlinear evolution and lattice equations, Adv. Comp. Math. 11 (1999), 55–80.

[GHE97] Ü. Göktaş, W. Hereman, and G. Erdmann, Computation of conserved densities for systems of nonlinear differential-difference equations, Phys. Lett. A 236 (1997), 30– 38.

[HG99] W. Hereman and Ü. Göktaş, Integrability Tests for Nonlinear Evolution Equations. Computer Algebra Systems: A Practical Guide, (M. Wester, ed.), Wiley, New York, 1999, pp. 211–232.

[HGCM98] W. Hereman, Ü. Göktaş, M. D. Colagrosso and A. J. Miller, Algorithmic integrability tests for nonlinear differential and lattice equations, Comp. Phys. Comm. 115 (1998), 428–446.

[HH03] M. S. Hickman and W. A. Hereman, Computation of densities and fluxes of nonlinear differential-difference equations, Proc. Roy. Soc. Lond. A 459 (2003), 2705–2729.

[HM04] P. E. Hydon and E. L. Mansfield, A variational complex for difference equations, Found. Comp. Math. (2004), in press.

[KvM75] M. Kac and P. van Moerbeke, On an explicitly soluble system of nonlinear differential equations related to certain Toda lattices, Adv. Math. 16 (1975), 160–169.

[MH02] E. L. Mansfield and P. E. Hydon, On a variational complex for difference equations. Proc. NSF-CBMS Conf. on The Geometrical Study of Differential Equations, Comtemporary Mathematics 285, (J. A. Leslie and T. P. Robart, eds.), American Mathematical Society, Providence, Rhode Island, 2002, pp. 121–129.

[MQ04] E. L. Mansfield and G. R. W. Quispel, Towards a variational complex for the finite element method, present proceedings.

[M00] A. G. Meshkov, Computer package for investigation of the complete integrability, Proc. Inst. Maths. NAS Ukraine 30 (2000), 35–46.

[OZF93] W. Oevel, H. Zhang and B. Fuchssteiner, Mastersymmetries and multi-Hamiltonian formulations for some integrable lattice systems, Prog. Theor. Phys. 81 (1989), 294– 308.

[O93] P. J. Olver, Applications of lie groups to differential equations, 2nd. ed., Graduate Texts in Mathematics 107, Springer, New York, 1993.

[SW01a] J. A. Sanders and J. P. Wang, Integrable systems and their recursion operators, Nonlinear Anal. 47 (2001), 5213–5240.

[SW01b] , On recursion operators, Physica D 149 (2001), 1–10.

[T81] M. Toda, Theory of nonlinear lattices, Springer, Berlin, 1981.

[W98] J. P. Wang, Symmetries and conservation laws of evolution equations, Ph.D. Thesis, Thomas Stieltjes Institute for Mathematics, Amsterdam, 1998.

[ZTOF93] H. Zhang, G. Tu, W. Oevel, and B. Fuchssteiner, Symmetries, conserved quantities, and hierarchies for some lattice systems with soliton structure, J. Math. Phys. 32 (1991), 1908–1918.

Department of Mathematical and Computer Sciences, Colorado School of Mines, Golden, Colorado 80401-1887

E-mail address: whereman@mines.edu

Department of Mathematics, Faculty of Sciences, Free University, De Boelelaan 1081a, 1081 HV Amsterdam, The Netherlands

E-mail address: jansa@cs.vu.nl

Physics Department 103-33, California Institute of Technology, Pasadena, California 91125

E-mail address: jack@its.caltech.edu

Department of Mathematics, Brock University, St. Catharines, Ontario, Canada L2S 3A1

E-mail address: jwang@spartan.ac.BrockU.ca