Direct Methods and Symbolic Software for Conservation Laws of Nonlinear Equations

. Direct Metho ds and Sym b olic Soft w are for Conser v ation La w s of Nonlinear Equations Willy Hereman, P a ul J. Adams, Holly L. Eklund Departmen t of Mathematical and Computer Sciences Colorado Sc ho o l of Mines Golden, CO 804 0 1-1887, U.S.A. Mark S. Hic kman Departmen t of Mathematics and Statistics Univ ersit y of Canterbury , Priv at e Bag 4800 8140 Christc hurc h, New Zealand Barend M. Herbst Departmen t of Mathematical Sciences Applied Mathematics D ivision General Engineering Building Univ ersit y of Stellenbosc h Priv ate Bag X1, 7602 Stellen b osch, South Africa. In Memor y of Mar tin D. Kruskal (1925-2006) Researc h Supp orted in P art by the South African National Researc h F oundation (NRF) under Grant No. F A20070325 00003. Man uscript prepared for Adv ances in Nonlinear W av es and Sym b olic Computation Ed.: Zhen ya Y an, Nov a Science Publishers, New Y ork, U.S.A. (2008). 2000 Mat hematics Sub ject Classiﬁcation Primary: 37K05, 37J35, 37K10; Secondary: 37Q35, 35Q58, 37K40 , 37K60 Abstract W e presen t direct metho ds, algorithms, and sym b olic softw are for the computation o f con- serv ation law s of nonlinear pa rtial diﬀeren tia l equations (PDEs) and diﬀeren tial-diﬀerence equations (DD Es). Our me t ho d for PDEs is based on calculus, linear algebra, and v ariationa l calculus. Fir st, w e c ompute the dilation symmetries of the give n nonlinear system. Next, we build a candidate densit y as a linear com bination with undetermined co eﬃcien ts of terms t ha t are scaling in v ariant. The v ariational deriv ative (Euler op erator) is used t o deriv e a linear system for the undetermined coeﬃcien ts. This system is then analyzed and solv ed. Finally , w e compute the ﬂux with the ho mo t o p y op erat or. The metho d is applied to nonlinear PD Es in (1 + 1) dimens ions with p olynomial nonlinear- ities which include t he K o rtew eg-de V ries (K dV), Boussinesq , and Drinfel’d-Sok o lo v-Wilson equations. An adaptat io n of the metho d is applied to PDEs with tra nscenden ta l nonlineari- ties. Examples include the sine-Go rdon, sinh-Gordo n, and Liouville equations. F or equations in lab orato ry coo rdinates, the coeﬃcien ts of the candidate densit y are undetermined functions whic h m ust satisfy a mixed linear system o f a lg ebraic and ordinary diﬀerential equations. F or the computation of conserv a t io n laws of nonlinear DDEs we use a splitting of the iden tity op erat o r. This metho d is more eﬃcien t that an appro a c h based on the discrete Euler and homotop y op erators. W e apply the metho d of undetermined co eﬃcien ts to the Kac-v an Mo erb ek e, T o da, and Ablowitz-Ladik lattices. T o o ve rcome the shortcomings o f the undetermined co eﬃcien t tec hnique, w e designed a new metho d that ﬁrst calculates the leading order term and then the required terms of lo w er o rder. That metho d, whic h is no longer restricted to p olynomial conserv ation laws , is applied to discretizations of the KdV and mo diﬁed KdV equations, a nd a com bination thereof. Additional examples includ e lattices due to Bogoy avlens kii, Belo v-Chaltikian, and Blaszak-Marciniak. The undetermined co eﬃcien t metho ds for PD Es a nd DDEs ha v e b een implemen ted in Mathematica . The co de TransPDEDens ityFlux.m computes densities and ﬂuxes of systems of PDEs with or without transcenden tal nonlinearities. The co de DDEDens ityFlux.m do es the same for p olynomial nonlinear DD Es. Sta rting from the leading o r der terms, the new Maple library discrete computes densities and ﬂuxes of nonlinear DDEs. The soft w a re can be used to answ er in tegrabilit y questions and to gain insigh t in the ph ys- ical and mathematical prop erties of no nlinear mo dels. When applied to nonlinear systems with parameters, the soft w are computes the conditions on the parameters for conserv atio n la ws to exist. The existenc e o f a hierarch y of conserv ation la ws is a predictor for complete in tegrability of the system and its solv abilit y with the In ve rse Scattering T ra nsform. Con te n ts 1 In tr o duction 3 2 The Most F amous Example in Historic al Perspective 6 3 The Metho d of Undetermined Co eﬃcien ts 8 3.1 Dilation In v ariance of Nonlinear PDEs . . . . . . . . . . . . . . . . . . . . . . 8 1 3.2 The Metho d of Undetermined Co eﬃcien ts Applied to a Scalar Nonlinear PDE 9 4 T o ols from t he Calculus of V ariations and Diﬀeren tial Geometry 10 4.1 The Contin uous V a riational Deriv ativ e (Euler Op erator) . . . . . . . . . . . . 11 4.2 The Contin uous Homotop y Op erator . . . . . . . . . . . . . . . . . . . . . . . 12 5 Conserv ation La ws of N onlinear Systems of Polynomial PDE s 13 5.1 T o ols for Systems of Ev olution Equations . . . . . . . . . . . . . . . . . . . . . 13 5.2 The D r inf el’d-Sok olov-Wilson System: Dilation In v ariance and Conserv a tion La ws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 5.3 Computation of a Conserv atio n Law of the D r infel’d-Sok olov Wilson System . 15 5.4 The Bo ussines q Equation: Dilatio n Inv ariance and Conserv a t io n Laws . . . . . 16 5.5 Computation of a Conserv atio n Law for the Boussinesq System . . . . . . . . 17 6 Conserv ation La ws of Systems of PDEs with T r anscenden tal N onlinearities 19 6.1 The sine-G o rdon Equation: Dilatio n Inv ariance and Conserv a t io n Laws . . . . 19 6.2 Computation of a Conserv atio n Law for the sine-Go r do n System . . . . . . . . 20 7 Conserv ation La ws of Scalar E quations with T ranscenden t al and Mixed Deriv ative T erms 22 7.1 The sine-G o rdon Equation in Light-Cone Co ordinates . . . . . . . . . . . . . . 22 7.2 Examples of Equations with T ranscenden tal Nonlinearities . . . . . . . . . . . 23 8 Nonlinear DDEs and Conserv ation La ws 24 9 The Metho d of Undetermined Co eﬃcien ts for DDEs 27 9.1 A Classic Example: The Kac-v an Mo erb ek e Lattice . . . . . . . . . . . . . . . 27 9.2 The Metho d of Undetermined Co eﬃcien ts Applied to a Scalar Nonlinear D DE 28 10 Discret e E uler and Homotop y Op erators 30 11 Conserv ation La ws of Nonlinear Systems of DDE s 32 11.1 The T o da Latt ice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 11.2 Computation of a Conserv ation Law of the T o da Lattice . . . . . . . . . . . . 33 11.3 The Ablowitz -Ladik L a ttice . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 11.4 Computation of a Conserv ation Law of the Ablow itz-Ladik La t tice . . . . . . . 35 11.5 A “Divide a nd Conquer” Strategy . . . . . . . . . . . . . . . . . . . . . . . . . 36 12 A New Metho d t o Compute Conserv ation La ws of Nonlinear DDE s 37 12.1 Leading Order Ana lysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 12.2 A“Mo diﬁed” V olterra Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 13 The Gardner Lattice 40 13.1 D eriv ation of the G ardner Lat tice . . . . . . . . . . . . . . . . . . . . . . . . . 41 13.2 D ilation In v ariance of the Gardner Lattice . . . . . . . . . . . . . . . . . . . . 44 13.3 Conserv ation Laws of the G ardner La ttice . . . . . . . . . . . . . . . . . . . . 44 2 14 Additional Examples of Nonlinear DDE s 45 14.1 The Bogoy a vlenskii Lat t ice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 14.2 The Belov -Chaltikian La ttice . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 14.3 The Blaszak-Marciniak Lattices . . . . . . . . . . . . . . . . . . . . . . . . . . 46 15 Softw are to Compute Conserv ation La w s for PDEs and D D E s 47 15.1 Our Mathematica and Maple Softw are . . . . . . . . . . . . . . . . . . . . . . 47 15.2 Softw are P a c k ages of Other Researc hers . . . . . . . . . . . . . . . . . . . . . . 48 16 Summary 49 1 In tro ductio n This chapter fo cuses o n sym b o lic metho ds to compute p olynomial conserv ation laws o f par- tial diﬀerential equations (PDEs) in (1 + 1) dimensions and diﬀerential-diﬀerence equations (DDEs), whic h a r e semi-discrete lattices. F or the latter w e treat systems where time is con tinuous a nd the spatial v ariable has b een discretized. Nonlinear PDEs that admit conserv ation laws arise in man y disciplines of the applied sciences including ph ysical c hemistry , ﬂuid mec hanics, particle and quan t um ph ysics, plasma ph ysics, elasticity , gas dynamics, electromagnetism, magneto-hydro-dynamics, nonlinear o p- tics, a nd the bio-sciences. Conserv ation la ws are fundamen t al la ws of phy sics that maintain that a certain quan tity will not change in time during ph ysical pro cess es. F amiliar conser- v atio n laws include conserv ation o f momentum, mass (matter), electric charge, or energy . The con tinuit y equation of electromagnetic theory is an example of a conserv ation la w whic h relates c har ge to curren t. In ﬂuid dynamics, the con tin uit y equation express es conserv ation of mass, and in quantum mec hanics the conserv ation of proba bilit y of the densit y and ﬂux functions also yields a con tinuit y equation. There are many reasons to compute conserv ed densities and ﬂuxes of PDEs explicitly . In v arian ts often lead to new disco v eries as w as the case in soliton theory . One may w ant to v erify if conserv ed quantitie s of phy sical imp ortance (e.g. momentum , energy , Hamiltonians, en tropy , densit y , c harge) are in tact after constitutiv e relations ha ve b een added to close a system. F or PD Es with arbitrary parameters one may wish to compute conditions on the parameters so t ha t the mo del admits conserv ed quan tities. Conserv ed densities also facilitate the study of qualita t iv e prop erties of PDEs [86, 97], suc h as recursion o p erators, bi- or tri-Hamiltonian structures, and t he lik e. They often guide the c hoice of solution metho ds or rev eal the nature of sp ecial solutions. F or example, an inﬁnite sequence of conserv ed densities is a predictor of the existence of solitons [7] and complete integrabilit y [2] whic h means that the PDE can b e solv ed with t he In v erse Scattering T ransform (IST) metho d [2]. Conserv ed densities aid in the design o f num erical solv ers fo r PDEs [87, 88] and their sta- bilit y analysis (see references in [23]) . Indeed, semi-discretizations that conserv e discrete con- serv ed quantities lead to stable nume rical sc hemes that are free of no nlinear instabilities and blo wup. While solving D DEs, whic h ar ise in nonlinear net works and as semi-discretizations of PDEs, one should c hec k that their conserv ed quan tities indeed remain unc hanged as time steps are take n. 3 Computer algebra systems (CAS) like Math ematica , Maple , and REDU CE , can greatly assis t the computation of conserv ation law s of nonlinear PD Es and D DEs. Using CAS interactiv ely , one can mak e a judicious guess (ansatz) a nd ﬁnd a few simple densities and ﬂuxes. Y et, that approach is fruitless for complicated systems with non trivial conserv a t io n law s with increasing complexit y . F urthermore, completely in tegrable equations PDEs [2, 7, 74, 89] and DDEs [10, 75] admit inﬁnitely many indep enden t conserv ation la ws. Compu ting them is a c ha llenging task. It in volv es t edious computations whic h are prone to erro r if done with p en and pap er. Krusk al and collab ora tors demonstrated the complexities o f calculating conserv ation law s in their seminal pap ers [67 , 78, 79] on the Kortewe g-de V ries (KdV) equation from soliton theory [2, 7, 30]. W e use this historical example to introduce the metho d of undetermined co eﬃcien ts. In the ﬁrst part of this c hapter w e cov er the sym b olic computation of conserv ation laws of completely in tegrable PD Es in ( 1 + 1) dimensions (with indep enden t v ariables x and t ) . Our approac h [8, 38, 51, 53] uses the concept of dilation (scaling) in v aria nce and the metho d of undetermined co eﬃcien ts. Our metho d pro ceeds a s follows . First, build a candidate den- sit y as a linear com bination (with unde termined co eﬃcien ts) of “ building blo c ks” that are homogeneous under the scaling symmetry of the PDE. If no suc h symmetry exists, construct one by in tro ducing parameters with scaling. Next, use the Euler op era t o r (v ariatio na l deriv a- tiv e) to derive a linear alg ebraic system fo r the undetermined co eﬃcien t s. After the system is a nalyzed and solv ed, use the homotopy op erator to compute the ﬂux. When applied to systems with para meters, our co des can determine the conditions on the parameters so that a sequenc e of conserv ed densities exists. The metho d is a pplied to nonlinear PDEs in (1 + 1) dimensions with p olynomial terms whic h include the KdV, Boussinesq, and Drinfel’d-So kolo v-Wilson equations. An ada pt a tion of the metho d is applied to PDEs with transcenden tal nonlinearities. Examples include the sine-Gordon, sinh-Gordon, and L io uville equations. F or equations written in lab oratory co ordinates, the co eﬃcien t s of the candidate densit y are undetermined functions whic h mus t satisfy a mixed linear system o f a lg ebraic and ordinary diﬀerential equations (ODEs). Capitalizing on the analogy b et w een PDEs and D DEs, the second part o f this c hapter deals with the sym b olic computation of conserv ation la ws of nonlinear DDEs [34, 39 , 42, 51, 54, 57]. Again, w e use scaling symmetries a nd the metho d of undetermined co eﬃcien ts. One could use discrete v ersions of the Euler op erator (to verify exactness) and the homotop y op erator (to in v ert the forw a r d diﬀerenc e). Although these op erat ors migh t b e v aluable in theory , they are highly ineﬃcien t as to ols for the sym b olic computatio n of conserv atio n law s of DDEs. W e adv o cate the use of a “splitting and shifting” tec hnique, whic h allo ws us to compute densities and ﬂuxes sim ultaneously at minimal cost. The undetermined co eﬃcien t metho d for DDEs is illustrated with the Kac-v a n Mo erb eke, T o da, and Ablowitz-Ladik lattices. There is a fundamen tal diﬀerence b etw een the con tinuous and discrete cases in the w ay densities a re constructed. The total deriv ativ e has a weigh t whereas the shift op erator do es not. Consequen tly , a densit y of a PD E is b ounded in order with respect to the space v ariable. Unfortunately , there is no a priori b ound on the n um b er of shifts in the densit y , unless a leading order analysis is carried out. T o ov ercome this diﬃcult y and other shortcomings of the undetermined co eﬃcien t metho d, w e presen t a new metho d to compute conserv ed densities of DD Es. That metho d no longer uses dilation inv ariance and is no longer restricted to p olynomial conserv ation la ws. Instead of building a candidate densit y with undetermined 4 co eﬃcien ts, one ﬁrst computes the leading order term in the densit y and, secondly , generates the required terms o f lo wer o rder. The metho d is fast a nd eﬃcien t since unnecessary terms are nev er computed. The new metho d is illustrated using a mo diﬁed V olterra lattice as an example. The new metho d p erforms exceedin gly w ell when applied to lattices due to Bogo y av lenskii, Belo v-Chaltikian, and Blaszak-Marciniak. The new metho d is also applied to completely integrable discretizations of the KdV and mo diﬁed KdV ( mKdV) equations, and a combination thereof, known as the Ga r dner equation. Starting from a discretized eigen v alue problem, we ﬁrst deriv e the Gardner latt ice and then compute conserv ation laws . There are sev eral metho ds (see [51]) to compute conserv ation law s of nonlinear PDEs and D DEs. Some metho ds use a generating function [2, 7], whic h requires the knowledge of k ey pieces of the IST. Another common approac h uses the link b et we en conserv ation laws and symme tries as stated in No ether’s theorem [14, 15, 66, 81]. How ev er, the computation of generalized (v ariational) symmetries, though alg orithmic, is as daunting a task a s the direct computation of conserv atio n laws . Most of t he more algorithmic metho ds [12, 13, 20, 25, 6 3 , 101], require the solutio n o f a determining system of OD Es or PDEs. Despite their p o wer, only a few of t hese metho ds hav e b een implemen ted in CAS. W e dev ote a section t o sym b olic soft w are for the computation o f conserv ation la ws. Additional reviews can b e found in [38, 51, 101]. Ov er the past decade, in collab or a tion with studen ts and researc hers, w e ha ve designed and implemen ted direct algorithms for the computation of conserv atio n laws of nonlinear PDEs and D DEs. W e purp osely av oid No ether’s theorem, pre-know ledge of symmetries, and a La- grangian form ulat ion. Neither do w e use diﬀerential f orms or adv anced diﬀeren tial-geometric to ols. Instead, w e concen trate o n the undetermined co eﬃcien t metho d for PDEs a nd DD Es, whic h uses to ols fro m calculus , linear a lg ebra, and the v ariational calculus. Therefore, the metho d is easy to impleme n t in Mat hematica and easy to use b y scien t ists and enginee rs. The co de TransPDEDensityFlux .m computes densities and ﬂuxes of systems o f PDEs with or without tra nscenden tal nonlinearities. The co de DDEDensityFl ux.m do es the same fo r p olynomial nonlinear DDEs. Starting from the leading order terms, the new Maple libra ry discrete computes densities and ﬂuxes of nonlinear DDEs ve ry eﬃcien tly . The soft ware can th us b e used to answ er inte grabilit y ques tions and to gain insigh t in the ph ysical and mathematical prop erties o f nonlinear mo dels. Our soft w are is in the public domain. The Mathematica pac k ages and noteb o oks are a v ailable at [48] and Hickman’s co de in Maple is av ailable at [56]. W e are curren tly w orking on a comprehensiv e pac k age to compute conserv at io n law s of PDEs in multiple space dimensions [45, 51, 83]. P art I: P artial Diﬀe ren tial Equation s in (1+ 1) Dimen sions In this ﬁrst part w e cov er PD Es in (1 + 1) dimensions , that is, PDEs in one space v ariable and time. Starting from a historical example, w e in tro duce t he concept of dilation inv ariance and use the metho d of undetermined co eﬃcien ts to compute conserv ation la ws of ev olution equations. Later on, we adapt the metho d of undetermined co eﬃcien ts to cov er PDEs with transcenden tal terms. 5 2 The Most F amous Example in Hist orical P ersp ectiv e The story of conserv ation laws for nonlinear PDEs b egins with the disco v ery of an inﬁnite n umber of conserv ation laws o f the ubiquitous Kortewe g-de V ries equation whic h mo dels a v ariety of nonlinear w av e phenomena, including shallo w water w a ves [46] and ion-acoustic w av es in pla smas [2 , 7, 30]. The KdV equation can b e recast in dimensionless v ariables as u t + α uu x + u 3 x = 0 , (1) where the subscripts denote partial deriv ativ es, i.e. u t = ∂ u ∂ t , u x = ∂ u ∂ x , and u 3 x = ∂ 3 u ∂ x 3 . The parameter α can b e scaled to a n y real num b er. Commonly used v alues are α = ± 1 o r α = ± 6 . Equation (1) is an example of a scalar (1 + 1) − dimensional ev olution equation, u t = F ( x, t, u, u x , u 2 x , · · · , u nx ) , (2) of o rder n in the indep enden t space v ariable x and of ﬁrst o rder in time t. Obvious ly , the dep enden t v ariable is u ( x, t ) . If parameters are presen t in (2), they will b e denoted by low er- case G reek letters. A c onservation law of (2) is of the form D t ρ + D x J = 0 , (3) whic h is satisﬁed for all solutions u ( x, t ) o f the PDE. In ph ysics, ρ is called the c onserve d density (or c harg e); J is the associated ﬂux (or curren t). In general, b oth are diﬀeren tial functions (functionals), i.e. functions of x, t, u, a nd partia l deriv a tiv es of u with resp ect to x. In (3), D x denotes the total deriv ativ e with resp ect to x, that is, D x J = ∂ J ∂ x + N X k =0 ∂ J ∂ u k x u ( k +1) x , (4) where N is the order of J, and D t is the total deriv ativ e with resp ect to t, deﬁned b y D t ρ = ∂ ρ ∂ t + ρ ′ [ u t ] = ∂ ρ ∂ t + M X k =0 ∂ ρ ∂ u k x D k x u t , (5) where ρ ′ [ u t ] is t he F r ´ ec het deriv ative of ρ in the direction of u t and M is the order of ρ. The densities ρ (1) = u and ρ (2) = u 2 of (1) w ere long known. In 19 6 5, Whitham [98] ha d found a third densit y , ρ (3) = u 3 − 3 α u 2 x , whic h, in the con text of w ater w av es, corresp onds to Boussinesq’s momen t of instabilit y [7 6]. One can readily v erify that D t ( u ) + D x ( 1 2 αu 2 + u 2 x ) = 0 , (6) D t ( u 2 ) + D x ( 2 3 αu 3 − u 2 x + 2 uu 2 x ) = 0 , (7) D t ( u 3 − 3 α u 2 x ) + D x ( 3 4 αu 4 − 6 u u 2 x + 3 u 2 u 2 x + 3 α u 2 2 x − 6 α u x u 3 x ) = 0 . (8) 6 Indeed, (6) is the K dV equation written as a conserv a tion law; (7) is obtained after multiply ing (1) b y 2 u ; (8) requires more work. Hence, the ﬁrst three densit y-ﬂux pairs of (1) are ρ (1) = u , J (1) = 1 2 αu 2 + u 2 x , (9) ρ (2) = u 2 , J (2) = 2 3 αu 3 − u 2 x + 2 uu 2 x , (10) ρ (3) = u 3 − 3 α u 2 x , J (3) = 3 4 αu 4 − 6 u u 2 x + 3 u 2 u 2 x + 3 α u 2 2 x − 6 α u x u 3 x . (11) In tegrals of mot io n readily follow from the densities. Indeed, assuming that J v anishes at inﬁnit y (for example due to suﬃc ien tly fast deca y of u a nd its x deriv ativ es), up on in tegration of (3) with resp ect to x one obtains that P = Z ∞ −∞ ρ dx (12) is c ons tant in time. Such constan ts of motion also arise when u is p erio dic, in which case one in tegrates o v er the ﬁnite p erio d. Dep ending on the ph ysical setting, the ﬁrst f ew constan ts of motion (i.e. integrals (12 ) ) express conserv ation of mass, momen tum, and energy . Martin Krusk al and p ostdo ctoral fellow Nor man Za busky disco v ered the fo urth and ﬁfth densities for the KdV equation [111]. Ho w ev er, they failed in ﬁnding a sixth conserv ation la w due to an algebraic mistak e in their computat io ns. Krusk al ask ed Ro b ert Miura, also p ostdo ctoral fellow at the Princeton Plasma Phys ics La b oratory at New Jersey , to searc h for further conserv atio n la ws of the KdV equation. Miura [78] computed the sev enth conserv ation la w. After correcting the mistak e men t io ned b efor e, he also found the sixth and ev entually three a dditio nal conserv a tion laws . Rumor [80] has it that in the summer of 1966 Miura w en t up in to the Cana dia n Ro c kies and returned from the moun ta ins with the ﬁrst 10 conserv ation la ws of the KdV equation engrav ed in his noteb o ok. T his biblical metaphor probably do es not do justice to Miura’s intens e and t edious w ork with p en and pap er. With ten conserv ation la ws in hand, it w as conj ectured that the KdV equation had an inﬁnite seque nce of conserv ation law s, later pro v en to b e true [67, 79]. Aﬁcionados of explic it form ulas can ﬁnd the ﬁrst ten densities (and sev en of the a sso ciated ﬂuxes) in [79] and the elev enth densit y (with 45 t erms) in [6 7], where a recursion form ula is giv en to generate all furt her conserv ed densities. As an aside, in 1966 the ﬁrst ﬁv e conserv ed densities w ere computed on an IBM 7094 computer with F ORMAC , an early CAS. The sixth densit y could no longer b e computed b ecause the av ailable stora ge space w a s exceeded. In con trast, using a metho d of undetermined co eﬃcien ts, the ﬁrst elev en densities w ere computed in 1969 on a AEC CDC-66 0 0 computer in a record time of 2.2 seconds. D ue to limitations in handling large in tegers, the computer could no t correctly pro duce a n y further densities. Undoubtedly , the discov ery of conserv ation laws pla ye d a piv otal role in the comprehensiv e study of the prop erties and solutions of nonlinear completely inte grable PD Es (lik e the K dV equation) and the deve lopmen t of the IST (see e.g. [80] for the history). Cliﬀord G ardner, John Gr eene, Martin Krusk al, and Rob ert Miura receiv ed the 2006 L eroy P . Steele Prize [11 5 ], a warded b y the American Mathematical So ciet y , for their seminal con tribution to researc h on the KdV equation. In turn, Martin K r usk al has receiv ed n umerous honors and a w ards [114] for his fundamen tal con tributions to t he understanding of in tegrable systems and solito n theory . This c hapter is dedicated to Martin Krusk al (192 5-2006). 7 3 The Metho d of Undetermine d Co eﬃcien ts W e no w sk etch t he method of undetermined coeﬃcien t s to compute conserv ation la ws [38 , 96], whic h draws on ideas and observ ations in b efore mentioned w ork b y Krusk al and collab orators. 3.1 Dilation I nv ariance of Non linear PDEs Crucial to the computation of conserv ation la ws is that (3 ) m ust hold on the PDE. This is ac hiev ed by substituting u t (and u tx , u txx , etc.) from ( 1 ) in the ev aluation of (6)-( 8 ) and in all subse quen t conserv ation la ws of degree larger than 3. The elimination of a ll t − deriv at ives of u in f a vor of x deriv a tiv es ha s tw o imp ortan t conse quences : (i) an y symmetry of the PDE, in particular, the dilation symmetry , will b e a dopted b y the conserv ation law, (ii) once D t is computed and ev aluated on the PDE, t b ecomes a parameter in the computation of the ﬂux. W e will ﬁrst inv estigate the dilation (scaling) symmetry of ev olution equations. The KdV equation is dilation inv a riant under the scaling symmetry ( t, x, u ) → ( λ − 3 t, λ − 1 x, λ 2 u ) , (13) where λ is an arbitrary parameter. Indeed, after a c hange of v ariables with ˜ t = λ − 3 t, ˜ x = λ − 1 x, ˜ u = λ 2 u, and cancellation of a common factor λ 5 , the KdV f or ˜ u ( ˜ x, ˜ t ) arises. The dilation symmetry of (1) can b e expresse d as u ∼ ∂ 2 ∂ x 2 , ∂ ∂ t ∼ ∂ 3 ∂ x 3 , (14) whic h means that u correspo nds to t w o x − deriv ativ es and the time deriv ativ e corresponds to three x − deriv ativ es. If w e deﬁne the weight , W , of a v ariable (or op erator) as the exponent of λ in (13), then W ( x ) = − 1 or W ( ∂ ∂ x ) = 1 ; W ( t ) = − 3 or W ( ∂ ∂ t ) = 3 , and W ( u ) = 2 . All w eights of dep enden t v a riables and the w eights of ∂ /∂ x , ∂ /∂ t, are assumed to b e non-negativ e and rationa l. The r an k of a monomial is deﬁned as the total w eight o f the monomial. Such monomials ma y in v olv e the indep enden t and dependent v ariables and the op erators ∂ ∂ x , D x , ∂ ∂ t , and D t . Ranks m ust b e p o sitiv e in t egers or p o sitiv e ra t io nal num b ers. An express io n (or equation) is uniform in r ank if its monomial terms ha v e equal rank. F or example, (1) is uniform in ra nk since eac h o f the three terms has rank 5 . Con v ersely , if one do es not kno w the dilation symmetry of (1), then it can b e readily computed b y requiring that (1) is uniform in ra nk. Indeed , setting W ( ∂ /∂ x ) = 1 and equating the ranks of the three t erms in (1) giv es W ( u ) + W ( ∂ ∂ t ) = 2 W ( u ) + 1 = W ( u ) + 3 , (15) whic h yields W ( u ) = 2 , W ( ∂ /∂ t ) = 3 , and, in turn, conﬁrms ( 1 3). So , r equiring uniformity in r ank of a PDE allo ws one to compute the w eights of the v ariables (and thus the scaling symmetry) with linear a lgebra. Dilation symmetries, whic h are sp ecial Lie-p oint symmetries, are common to man y non- linear PDEs. Needless to sa y , not ev ery PDE is dilation in v arian t, but non-uniform PDEs can b e made uniform by extending the set of dep enden t v ar ia bles with auxiliary parameters with 8 appropriate w eigh ts. Up on completion of the computations one can set these parameters to one. In what f ollo ws, w e set W ( ∂ /∂ x ) = W ( D x ) = 1 and W ( ∂ /∂ t ) = W ( D t ) . Applied to (6), rank ρ (1) = 2, ra nk J (1) = 4 . Hence, rank ( D t ρ (1) ) = rank ( D x J (1) ) = 5 . Therefore, (6) is uniform of rank 5 . In (7), rank ρ (2) = 4 and rank J (2) = 6 , consequen tly , (7) is uniform o f rank 7 . In (11), eac h term in ρ (3) has rank 6 a nd each term in J (3) has rank 8 . Consequen tly , rank ( D t ρ (3) ) = rank ( D x J (3) ) = 9 , whic h mak es (8) is uniform of rank 9 . All densities of (1) are uniform in rank and so are the asso ciated ﬂuxes and the conserv ation law s. Equation (1) also has densit y-ﬂux pa irs that dep end explicitly on t and x ; for example, ˜ ρ = tu 2 + 2 α xu, ˜ J = t ( 2 3 αu 3 − u 2 x + 2 uu 2 x ) − x  u 2 − 2 α u 2 x  + 2 α u x . (16) Since W ( x ) = − 1 and W ( t ) = − 3 , one has r ank ˜ ρ = 1 , and rank ˜ J = 3 . The metho ds and algorithms discussed in subsequen t sections hav e b een adapted to compute densities and ﬂuxes explicitly dep enden t on x and t. Instead of addressing this issue in this chapter, w e refer the reader to [38, 53]. 3.2 The Metho d of Undetermined Co eﬃcien ts A p plied to a Scalar Nonlinear PDE W e outline ho w densities and ﬂuxes can b e constructed f o r a scalar ev olution equation (2). T o kee p matters transparen t, w e illustrate the steps for the KdV equation r esulting in ρ (3) of rank R = 6 with asso ciated ﬂux J (3) of rank 8 , b oth listed in (11). The to ols needed for the computations will b e presen ted in the next section. • Select the rank R of ρ. Mak e a list, R , of all monomials in u and its x -deriv at ives so that eac h mono mial has rank R. This can b e done as follo ws. Starting from the set V of dep enden t v aria bles (including parameters with w eight, when applicable), mak e a set M of all non- constan t monomials of rank R or less (but without x − deriv ativ es). Next, for eac h term in M , in tro duce the right n um b er o f x − deriv ativ es to adjust the rank of that term. Distribute the x − deriv ativ es, strip oﬀ the numerical co eﬃcien ts, and gather the resulting terms in a set R . F or the KdV equation and R = 6, V = { u } and M = { u 3 , u 2 , u } . Since u 3 , u 2 , a nd u ha ve r a nks 6 , 4 a nd 2 , resp ectiv ely , one computes ∂ 0 u 3 ∂ x 0 = u 3 , ∂ 2 u 2 ∂ x 2 = 2 u 2 x + 2 uu 2 x , ∂ 4 u ∂ x 4 = u 4 x . (17) Ignoring n umerical co eﬃcien ts in the right hand sides of the equations in (17), one gets R = { u 3 , u 2 x , uu 2 x , u 4 x } . • R emov e from R all monomials that are tot a l x − deriv a t iv es. Also remov e all “equiv alen t” monomials, i.e. the monomials that diﬀer fro m another b y a tota l x − deriv ativ e, keep ing the monomial of low est order. Call the resulting set S . In our example, u 4 x m ust b e remo v ed (b ecause u 4 x = D x u 3 x ) and uu 2 x m ust b e remo v ed since uu 2 x and u 2 x are equiv alent. Indeed, uu 2 x = D x ( uu x ) − u 2 x . Th us, S = { u 3 , u 2 x } . • Linearly combine the monomials in S with constan t undetermined co eﬃcien ts c i to obtain the candidate ρ. Con tin uing with t he example, ρ = c 1 u 3 + c 2 u 2 x , (18) 9 whic h is of ﬁrst order in x. • Using (5), compute D t ρ. Applied to (18) where M = 1 , one gets D t ρ = (3 c 1 u 2 I + 2 c 2 u x D x )[ u t ] . (19) As usual, D 0 x = I is the iden tity op erator. • Ev alua te − D t ρ on the PDE (2 ) b y replacing u t b y F . The result is a diﬀeren tial function E in whic h t is a para meter. F or the K dV equation (1), F = − ( αuu x + u 3 x ) . Af t er rev ersing the sign, the ev aluated form of (1 9 ) is E = (3 c 1 u 2 I + 2 c 2 u x D x )( αuu x + u 3 x ) = 3 c 1 αu 3 u x + 2 c 2 αu 3 x + 2 c 2 αuu x u 2 x + 3 c 1 u 2 u 3 x + 2 c 2 u x u 4 x . (20) • T o obtain a conserv a tion law, E m ust b e a total deriv ativ e. Starting with highest o rders, rep eatedly integrate E b y pa rts. Doing so, allow s one to write E as the sum o f a t otal x − deriv ativ e, D x J, and a non-integrable par t (i.e. the o bstructing terms). J is the ( candidat e) ﬂux with rank J = R + W ( D t ) − 1 . In tegrat ion b y pa rts of (20) giv es E = D x ( 3 4 c 1 αu 4 + 3 c 1 u 2 u 2 x + c 2 αuu 2 x + 2 c 2 u x u 3 x ) − 6 c 1 uu x u 2 x − 2 c 2 u 2 x u 3 x + c 2 αu 3 x = D x ( 3 4 c 1 αu 4 − 3 c 1 uu 2 x + c 2 αuu 2 x + 3 c 1 u 2 u 2 x + 2 c 2 u x u 3 x − c 2 u 2 2 x ) + (3 c 1 + c 2 α ) u 3 x . (2 1) The candidate ﬂux therefore is J = 3 4 c 1 αu 4 − 3 c 1 uu 2 x + c 2 αuu 2 x + 3 c 1 u 2 u 2 x + 2 c 2 u x u 3 x − c 2 u 2 2 x . (22) • Equate the coeﬃcien ts of the obstructing terms to zero. Solve the linear sy stem for the undetermined co eﬃcien ts c i . In the example (3 c 1 + c 2 α ) u 3 x is the only obstructing term whic h v anishes f or c 2 = − 3 α c 1 , where c 1 is arbitrary . • Substitute the co eﬃcien ts c i in to ρ and J to obt a in the ﬁnal fo rms of the densit y a nd asso ciated ﬂux (with a common arbitr a ry factor whic h can b e set to 1). Setting c 1 = 1 and substituting c 2 = − 3 α in to (18) and (22) yields ρ (3) and J (3) as giv en in (11). Constructing “minimal” densities, i.e. densities whic h are free of equiv alen t terms and to tal deriv atives terms, b ecomes c hallenging if the rank R is high. F urthermore, integration b y parts is cum b ersome and prone to mistak es if done b y hand. Moreo v er, it w ould b e adv an- tageous if the in tegratio n by parts could b e p ostp oned until the undetermined co eﬃcien ts c i ha ve b een computed and substituted in E . Ideally , the computations of the densit y and the ﬂux should b e decoupled. There is a need for computationa l to ols to address these issues, in particular, if one wan ts to compute conserv ation laws of systems of ev olution equations. 4 T o ols fro m the Calc u lus of V ariatio ns and Di ﬀe ren tial Geometry A sc alar diﬀeren tial function E of order M is called exact (in tegra ble) if and only if there exists a scalar diﬀeren t ial function J of order M − 1 suc h that E = D x J. Ob viously , J = 10 D − 1 x E = R E d x is then the primitiv e (or integral) of E . Tw o questions arise: (i) How can one test whether or not E is exact? (ii) If E is exact, how can one compute J without using standard in tegratio n by parts? T o answe r the ﬁrst question w e will use the v ariational deriv ative (Euler op era t or) fr om the calculus of v ar iations. T o p erform in tegration b y part s w e will use the homotop y op erato r fr o m diﬀeren tial g eometry . 4.1 The Cont in uous V ariational Deriv ativ e (Euler Op erator) The con t in uous variational derivative , also called the Euler op er ator of or der z e r o , L (0) u ( x ) , f o r v aria ble u ( x ) is deﬁned [81] by L (0) u ( x ) E = M X k =0 ( − D x ) k ∂ E ∂ u k x = ∂ E ∂ u − D x ∂ E ∂ u x + D 2 x ∂ E ∂ u 2 x − D 3 x ∂ E ∂ u 3 x + · · · + ( − 1) M D M x ∂ E ∂ u M x , (23) where E is a diﬀerential function in u ( x ) of or der M . A necessary and suﬃcien t condition for a diﬀerential function E to b e exact is that L (0) u ( x ) E ≡ 0 . A pro of o f this statemen t is giv en in e.g. [67]. If L (0) u ( x ) E 6 = 0 , t hen E is not a total x − deriv ative due to obstructing t erms. Application 1 . Returning to (20), w e no w use the v ariational deriv ative to determine c 1 and c 2 so that E of order M = 4 will b e exact. Using nothing but diﬀeren tiations, w e readily compute L (0) u ( x ) E = ∂ E ∂ u − D x ∂ E ∂ u x + D 2 x ∂ E ∂ u 2 x − D 3 x ∂ E ∂ u 3 x + D 4 x ∂ E ∂ u 4 x = 9 c 1 αu 2 u x + 2 c 2 αu x u 2 x + 6 c 1 uu 3 x − D x (3 c 1 αu 3 + 6 c 2 αu 2 x + 2 c 2 αuu 2 x + 2 c 2 u 4 x ) + D 2 x (2 c 2 αuu x ) − D 3 x (3 c 1 u 2 ) + D 4 x (2 c 2 u x ) = (9 c 1 αu 2 u x + 2 c 2 αu x u 2 x + 6 c 1 uu 3 x ) − (9 c 1 αu 2 u x + 14 c 2 αu x u 2 x + 2 c 2 αuu 3 x +2 c 2 u 5 x ) + (6 c 2 αu x u 2 x + 2 c 2 αuu 3 x ) − (18 c 1 u x u 2 x + 6 c 1 uu 3 x ) + (2 c 2 u 5 x ) = − 6(3 c 1 + c 2 α ) u x u 2 x . (24) Note that the terms in u 2 u x , uu 3 x , and u 5 x dropp ed out. He nce, requiring that L (0) u ( x ) E ≡ 0 leads to 3 c 1 + c 2 α = 0 . Substituting c 1 = 1 , c 2 = − 3 α , into (18) yields ρ (3) in (11). Application 2 . It is paramoun t that the candidate densit y is free of total x − deriv ativ es and equiv alen t terms. If such terms we re presen t, they could b e mov ed into the ﬂux J, and their co eﬃcien ts c i w ould b e ar bitr ary . ρ (1) and ρ (2) , are e quivalent if and only if ρ (1) + k ρ (2) = D x J, for some J and non-zero scalar k . W e write ρ (1) ≡ ρ (2) . Clearly ρ is equiv alent to an y no n- zero m ultiple of itself and ρ ≡ 0 if and only if ρ is exact. Instead of w orking with diﬀeren t densities, w e in v estigate the equiv alence of terms t i in the same densit y . F or example, returning to t he set R = { u 3 , u 2 x , uu 2 x , u 4 x } , terms t 2 = u 2 x and t 3 = uu 2 x are equiv alent b ecause t 3 + t 2 = uu 2 x + u 2 x = D x ( uu x ) . The v ariational deriv ativ e can b e used to detect equiv alen t and exact terms. Indeed, note that v 1 = L (0) u ( x ) ( uu 2 x ) = 2 u 2 x and v 2 = L (0) u ( x ) u 2 x = − 2 u 2 x are linearly dep enden t. Also, 11 for t 4 = u 4 x = D x u 3 x one gets v 3 = L (0) u ( x ) u 4 x = 0 . T o w eed out the terms t i in R that are equiv alen t or total deriv ativ es, it suﬃces to c hec k the linear indep endence of their images v i under the Euler op erator. One can optimize this pr o cedure b y starting from a set R where some of the equiv alent and total deriv ativ es terms hav e b een remo ve d a priori . Indeed, in view of (17), one can ignore the highest-order terms (typically the last terms) in eac h of the righ t hand sides. Therefore, R = { u 3 , u 2 x } and, fo r this example, no furt her reduction would b e necessary . V a r io us algorithms are p ossible to construct minimal densities. D eta ils are g iven in [38, 51]. 4.2 The Cont in uous Homotop y O p erator W e no w discus s the homotopy op erator [12 , 13, 52, 81] whic h will allo w one to reduce the computation of J = D − 1 x E = R E d x to a single integral with resp ect to an auxiliary v a riable denoted b y λ (not to b e confused with λ in Section 3.1). Hence, the homotopy op erator circum v ents in tegration b y parts and reduces t he inv ersion of D x to a problem of single- v aria ble calculus. The ho m otopy op er a tor [81, p. 372] f or v ariable u ( x ) , acting on an exact expression E of order M , is giv en b y H u ( x ) E = Z 1 0 ( I u E ) [ λu ] dλ λ , (25) where the in tegrand I u E is given b y I u E = M X k =1 k − 1 X i =0 u ix ( − D x ) k − ( i +1) ! ∂ E ∂ u k x . (26) In (25), ( I u E )[ λu ] means that in I u E one replaces u → λu, u x → λu x , etc . This is a sp ecial case of the homotopy , λ ( u (1) − u (0) ) + u (0) , b etw een t wo p oints, u (0) = ( u 0 , u 0 x , u 0 2 x , · · · , u 0 M x ) and u (1) = ( u 1 , u 1 x , u 1 2 x , · · · , u 1 M x ) , in the jet space. F or our purp oses w e set u (0) = (0 , 0 , · · · , 0) and u (1) = ( u, u x , u 2 x , · · · , u M x ) . F ormula (26 ) is equiv alen t to the one in [52], whic h in turn is equiv alen t to the formula in terms of higher Euler op erators [45, 51 ]. Giv en an exact diﬀeren tial function E o f order M one has J = D − 1 x E = R E d x = H u ( x ) E . A pro of of this statemen t can b e fo und in [52]. Application . After substituting c 1 = 1 and c 2 = − 3 α in to (20) we obt a in the exact expression E = 3 αu 3 u x − 6 u 3 x − 6 uu x u 2 x + 3 u 2 u 3 x − 6 α u x u 4 x , (27) of order M = 4 . First, using (26), w e compute I u E = 4 X k =1 k − 1 X i =0 u ix ( − D x ) k − ( i +1) ! ∂ E ∂ u k x = ( u I ) ( ∂ E ∂ u x ) + ( − u D x + u x I )( ∂ E ∂ u 2 x ) + ( u D 2 x − u x D x + u 2 x I )( ∂ E ∂ u 3 x ) +( − u D 3 x + u x D 2 x − u 2 x D x + u 3 x I )( ∂ E ∂ u 4 x ) . (28) 12 After substitution of (2 7), o ne gets I u E = ( u I )(3 α u 3 + 18 u 2 x − 6 uu 2 x − 6 α u 4 x ) + ( − u D x + u x I )( − 6 uu x ) +( u D 2 x − u x D x + u 2 x I )(3 u 2 ) + ( − u D 3 x + u x D 2 x − u 2 x D x + u 3 x I )( − 6 α u x ) = 3 α u 4 − 18 uu 2 x + 9 u 2 u 2 x + 6 α u 2 2 x − 12 α u x u 3 x , (29) whic h has the correct terms of J (3) but incorrect co eﬃcien ts. Finally , using (25), J = H u ( x ) E = Z 1 0 ( I u E )[ λu ] dλ λ = Z 1 0  3 αλ 3 u 4 − 18 λ 2 uu 2 x + 9 λ 2 u 2 u 2 x + 6 α λu 2 2 x − 12 α λu x u 3 x  dλ = 3 4 αu 4 − 6 u u 2 x + 3 u 2 u 2 x + 3 α u 2 2 x − 6 α u x u 3 x , (30) whic h matc hes J (3) in (11). The crux of the homotopy op erator metho d [12 , 13, 26, 8 1] is that the in tegrat io n by parts of a diﬀ erential expression lik e (27), whic h inv olv es an arbitrary function u ( x ) and its x − deriv ativ es, can b e r educed to a standard integration of a p olynomial in λ. 5 Conser v ation La w s of Nonlinear Systems of P olyno- mial PDEs Th us far w e hav e dealt with the computation of densit y-ﬂux pairs of scalar ev olution equa- tions, with the KdV equation a s the leading example. In this section w e sho w how the metho d and to ols can b e generalized to cov er systems of ev olutio n eq uations. W e will use the Drinfel’d-Sokolo v-Wilson system and t he Boussinesq equation to illustrate the steps. 5.1 T o ols for Systems of Ev olution Equations F or diﬀeren tial functions (lik e densities and ﬂuxes) of tw o dep enden t v ariables ( u, v ) and their x − deriv ativ es, the total deriv at ives ar e D t ρ = ∂ ρ ∂ t + M 1 X k =0 ∂ ρ ∂ u k x D k x u t + M 2 X k =0 ∂ ρ ∂ v k x D k x v t , (31) D x J = ∂ J ∂ x + N 1 X k =0 ∂ J ∂ u k x u ( k +1) x + N 2 X k =0 ∂ J ∂ v k x v ( k +1) x , (32) where M 1 and M 2 are the (highest) orders of u and v in ρ, and N 1 and N 2 are the (highest) orders of u and v in J. 13 T o accommo date t w o dep enden t v ariables, we need Euler operato rs L (0) u ( x ) and L (0) v ( x ) for eac h dep enden t v ariable separately . F or brevit y , w e will use v ector no t ation, that is, L (0) u ( x ) E = ( L (0) u ( x ) E , L (0) v ( x ) E ) . Like wise, the homoto p y op erat o r in (25) m ust b e replaced by H u ( x ) E = Z 1 0 ( I u E + I v E )[ λ u ] dλ λ , (33) where I u E = M 1 X k =1 k − 1 X i =0 u ix ( − D x ) k − ( i +1) ! ∂ E ∂ u k x , (34) and I v E = M 2 X k =1 k − 1 X i =0 v ix ( − D x ) k − ( i +1) ! ∂ E ∂ v k x , (35) where M 1 , M 2 are the orders of E in u, v , resp ectiv ely . In (33), u → λu, u x → λu x , · · · , v → λv , v x → λv x , etc . 5.2 The Drinfel’d-Sok olo v-Wilson System: Dilation In v ariance and Conserv ation La ws W e consider a parameterized family of the Drinfel’d-Sok olo v-Wilson (DSW) equations u t + 3 v v x = 0 , v t + 2 uv x + α u x v + 2 v 3 x = 0 , (36) where α is a nonzero parameter. The system with α = 1 w as ﬁrst prop osed b y Drinf el’d and Sok olov [31, 32] and Wilson [9 9]. It can b e obtained [59] as a reduction of the Kadom tsev- P etviash vili equation (i.e. a t wo-dimensional v ersion of the KdV equation) and is a completely in tegrable system. In [109], Y ao and Li computed conserv ation la ws of (36 ), where they had in tro duced four a rbitrary co eﬃcien t s. Using sc ales on x, t, u and v , all but one co eﬃcien ts in (36) can b e scaled to an y r eal n umber. Therefore, to cov er the entire family of DSW equations it suﬃces to leav e one co eﬃcien t arbitr a ry , e.g. α in fro n t of u x v . T o compute the dilation symmetry of (36), w e assign w eights, W ( u ) and W ( v ) , to b oth dep enden t v aria bles and express that eac h equation separately mu st b e uniform in rank (i.e. the ranks of the equations in (36) ma y diﬀer from each other). F or the D SW equations (36), one has W ( u ) + W ( ∂ /∂ t ) = 2 W ( v ) + 1 , W ( v ) + W ( ∂ /∂ t ) = W ( u ) + W ( v ) + 1 = W ( v ) + 3 , (37) whic h yields W ( u ) = W ( v ) = 2 , W ( ∂ /∂ t ) = 3 . The DSW syste m (36) is th us inv arian t under the scaling symmetry ( x, t, u, v ) → ( λ − 1 x, λ − 3 t, λ 2 u, λ 2 v ) , (38) where λ is an arbitra ry scaling parameter. 14 The ﬁrst three densit y-ﬂux pairs for the DSW equations (36) a r e ρ (1) = u, J (1) = 3 2 v 2 , (39) ρ (2) = v , J (2) = 2( u v + v 2 x ) , if α = 2 , (40) ρ (3) = ( α − 1) u 2 + 3 2 v 2 , J (3) = 3( α uv 2 − v 2 x + 2 v v 2 x ) , (41) Both ρ (1) and ρ (2) ha ve rank 2; their ﬂuxes ha v e rank 4 . The pair ( ρ (1) , J (1) ) exists for any α, whereas ( ρ (2) , J (2) ) only exists if α = 2 . Densit y ρ (3) of rank 4 a nd ﬂux J (3) of rank 6 a re v alid for an y α. A t rank R = 6 , ρ (4) = ( α + 1)( α − 2) u 3 − 9 2 ( α + 1) uv 2 − 3 2 ( α − 2) u 2 x − 27 2 v 2 x . (42) The corresp onding ﬂux (not sho wn) has 7 terms. A t rank R = 8 , ρ (5) = u 4 − 9 2 u 2 v 2 − 27 8 v 4 − 9 2 uu 2 x + 3 4 u 2 2 x + 45 2 v u x v x + 27 u v 2 x − 81 4 v 2 2 x , (43) pro vided α = 1 . The correspo nding ﬂux (not sho wn) has 15 terms. There exists a densit y-ﬂux pair for all ev en ranks R ≤ 10 provided α = 1 , for whic h (36) is completely in tegrable. 5.3 Computation of a Conserv ation La w of the Drinfel’d-Sok olo v Wilson S ystem T o illustrate ho w the presence of a parameter, lik e α , aﬀects the computation of dens ities, w e compute ρ (1) and ρ (2) of rank R = 2 g iv en in (39) and (40). Step 1 : C onst r uct the form of t he densit y The set of dependen t v ar iables is V = { u, v } . Both elemen ts are of rank 2 so, no x − de riv atives are needed. Th us, M = R = S = { u, v } . Linearly combining the elemen ts in S gives ρ = c 1 u + c 2 v . Step 2 : C ompute t he undetermined co eﬃcien ts c i Ev aluating E = − D t ρ = − ( c 1 u t + c 2 v t ) on (36), yields E = 3 c 1 v v x + c 2 (2 uv x + α u x v + 2 v 3 x ) , (44) whic h will b e ex act if L (0) u ( x ) E = ( L (0) u ( x ) E , L (0) v ( x ) E ) ≡ (0 , 0) . Since E is of order M 1 = 1 in u and order M 2 = 3 in v , L (0) u ( x ) E = ∂ E ∂ u − D x ∂ E ∂ u x = 2 c 2 v x − D x ( c 2 αv ) = (2 − α ) c 2 v x , (45) and L (0) v ( x ) E = ∂ E ∂ v − D x ∂ E ∂ v x + D 2 x ∂ E ∂ v 2 x − D 3 x ∂ E ∂ v 3 x = 3 c 1 v x + c 2 αu x − D x (3 c 1 v + 2 c 2 u ) − D 3 x (2 c 2 ) = ( α − 2) c 2 u x . (46) 15 Both (45 ) and (46) will v anish iden tically if and only if ( α − 2) c 2 = 0 . This equation (with unkno wns c 1 and c 2 ) is parameterized b y α 6 = 0 . The solution algorithm [38] consid ers all branc hes of the solution and p ossible compatibilit y conditions. Setting c 1 = 1 , leads to either (i) c 2 = 0 if α 6 = 2 , or (ii) c 2 arbitrary if α = 2 . Setting c 2 = 1 leads to the compatibilit y condition, α = 2 , a nd c 1 arbitrary . Substituting the solutions in to ρ = c 1 u + c 2 v giv es ρ = u whic h is v a lid for an y α ; and ρ = u + c 2 v or ρ = c 1 u + v pro vided α = 2 . In other w ords, ρ (1) = u is the only dens it y of rank 2 for ar bitr ary v alues o f α . F or α = 2 there exist tw o indep enden t densities, ρ (1) = u and ρ (2) = v . Step 3 : C ompute t he asso ciated ﬂux J As an example, w e compute the ﬂux in ( 4 0) asso ciated with ρ (2) = v and α = 2 . In this case, c 1 = 0 , c 2 = 1 , for whic h (44 ) simpliﬁes in to E = 2 ( uv x + u x v + v 3 x ) , (47) whic h is of order M 1 = 1 in u a nd order M 2 = 3 in v . Using (34) and (35 ), we obtain I u E = u I ∂ E ∂ u x = u I (2 v ) = 2 uv , (48) and I v E = ( v I )( ∂ E ∂ v x ) + ( − v D x + v x I )( ∂ E ∂ v 2 x ) + ( v D 2 x − v x D x + v 2 x I )( ∂ E ∂ v 3 x ) = ( v I )(2 u ) + ( v D 2 x − v x D x + v 2 x I )(2) = 2 uv + 2 v 2 x . (49) Hence, using (33), J = H u ( x ) E = Z 1 0 ( I u E + I v E )[ λ u ] dλ λ = Z 1 0 (4 λuv + 2 v 2 x ) dλ = 2( u v + v 2 x ) , (50) whic h is J (2) in (40). The in tegratio n of (47) could easily b e done by hand. The homoto p y op erator metho d pays o ﬀ if the expression to b e inte grated has a large n umber of terms. 5.4 The Boussinesq Equation: Dilation In v ariance and Conserv a- tion L a ws The w av e equation, u 2 t − u 2 x + 3 u 2 x + 3 uu 2 x + α u 4 x = 0 , (51) for u ( x, t ) with real parameter α , was prop osed by Boussinesq t o describe surface w a ves in shallo w w ater [2]. F or what f o llo ws, w e rewrite (51) as a system of ev olution equations, u t + v x = 0 , v t + u x − 3 u u x − αu 3 x = 0 , (5 2) where v ( x, t ) is an auxiliary dep enden t v ariable. The Boussines q system (52) is not uniform in rank b ecause the terms u x and α u 3 x lead to an inconsisten t system of we igh t equations. T o circum v en t the pro blem we in tro duce an auxiliary parameter β with (unkno wn) w eight, and replace (52) by u t + v x = 0 , v t + β u x − 3 uu x − α u 3 x = 0 . (53) 16 Requiring uniformit y in rank, w e obtain (after some algebra) W ( u ) = 2 , W ( v ) = 3 , W ( β ) = 2 , W ( ∂ ∂ t ) = 2 . (54) Therefore, (53) is inv arian t under the scaling symmetry ( x, t, u, v , β ) → ( λ − 1 x, λ − 2 t, λ 2 u, λ 3 v , λ 2 β ) . (55) As the ab o ve example sho ws, a PDE that is not dilatio n inv arian t can b e made so by extending the se t of dep enden t v ariables with one or more auxiliary parameters with w eigh ts. Up on completion of the computatio ns one can set eac h of these para meters equal to 1 . The Bo ussinesq equation (51 ) has inﬁnitely man y conserv atio n la ws and is completely in tegrable [2, 7]. The ﬁrst four densit y-ﬂux pairs [8] for (53 ) a r e ρ (1) = u , J (1) = v , (56) ρ (2) = v , J (2) = β u − 3 2 u 2 − α u 2 x , (57) ρ (3) = u v , J (3) = 1 2 β u 2 − u 3 + 1 2 v 2 + 1 2 αu 2 x − α uu 2 x , (58) ρ (4) = β u 2 − u 3 + v 2 + α u 2 x , J (4) = 2 β uv − 3 u 2 v + 2 αu x v x − 2 α u 2 x v . (59) These densities are of ranks 2 , 3 , 5 and 6 , resp ectiv ely . The corresp onding ﬂuxes are of one rank higher. After setting β = 1 w e obtain the conserv ed quan tit ies of ( 5 2) ev en though initially this system was no t unifo r m in rank. 5.5 Computation of a Conserv ation La w for the Boussinesq Sys tem W e sho w the computation of ρ (4) and J (4) of ra nks 6 and 7 , resp ectiv ely . The presence of the auxiliary parameter β with w eigh t complicates matters. At a ﬁxed rank R, conserv ed densities corresp onding to lo w er ranks migh t app ear in the r esult. These lo w er-r a nk densities are easy to recognize for they are m ultiplied with arbitrary co eﬃcien ts c i . Conseq uen tly , when parameters with w eight are intro duced, the densities corresp onding to distinct r anks are no longer linearly independen t. As the example below will sho w, densities m ust b e split in to indep enden t pieces. Step 1 : C onst r uct the form of t he densit y Augmen t the set o f dep enden t v aria bles with the parameter β (with non-zero w eight). Hence, V = { u, v , β } . Construct M = { β 2 u, β u 2 , β u, β v , u 3 , u 2 , u, v 2 , v , uv } , whic h con tains all non- constan t monomials of (c hosen) ra nk 6 or less (without deriv ativ es). Next, for eac h term in M , in tro duce the right num b er of x -deriv ativ es so that eac h term has ra nk 6 . F or example, ∂ 2 β u ∂ x 2 = β u 2 x , ∂ 2 u 2 ∂ x 2 = 2 u 2 x + 2 uu 2 x , ∂ 4 u ∂ x 4 = u 4 x , ∂ ( u v ) ∂ x = v u x + uv x , etc.. (60) Gather the terms in t he righ t hand sides of the equations in (6 0) to get R = { β 2 u, β u 2 , u 3 , v 2 , v u x , u 2 x , β v x , uv x , β u 2 x , uu 2 x , v 3 x , u 4 x } . (61) 17 Using ( 2 3) and a similar form ula for v , for eve ry term t i in R w e compute v i = L (0) u ( x ) t i = ( L (0) u ( x ) t i , L (0) v ( x ) t i ) . If v i = (0 , 0 ) then t i is discarded and so is v i . If v i 6 = (0 , 0 ) w e v erify whether or no t v i is linearly indep enden t of the non-zero ve ctors v j , j = 1 , 2 , · · · , i − 1 . If independen t, the term t i is k ept, otherwise, t i is discarded and so is v i . Up on a pplicatio n of L (0) u ( x ) , the ﬁrst six terms in R lead to linearly indep enden t v ectors v 1 through v 6 . Therefore, t 1 through t 6 are kept (and so are the corresp onding ve ctors). F or t 7 = β v x w e compute v 7 = L (0) u ( x ) ( β v x ) = (0 , 0) . So, t 7 is discarded and so is v 7 . F or t 8 = uv x w e get v 8 = L (0) u ( x ) ( uv x ) = ( v x , − u x ) = − v 5 . So , t 8 is discarded and so is v 8 . Pro ceeding in a similar f ashion, t 9 , t 10 , t 11 and t 12 are discarded. Th us, R is replaced by S = { β 2 u, β u 2 , u 3 , v 2 , v u x , u 2 x } , (62) whic h is free of div ergences a nd div erg ence-equiv alen t terms. Ignor ing the highest-order terms (t ypically the last terms) in each of the righ t hand sides of the equations in (60) optimizes the pro cedure. Indeed, R w o uld ha v e had six instead of tw elv e terms. Coinciden tally , in this example no further eliminations w o uld b e needed to obtain S . Next, linearly com bine the terms in S to get ρ = c 1 β 2 u + c 2 β u 2 + c 3 u 3 + c 4 v 2 + c 5 v u x + c 6 u 2 x . (63) Step 2 : C ompute t he undetermined co eﬃcien ts c i Compute D t ρ. Here, ρ is of or der M 1 = 1 in u and order M 2 = 0 in v . Hence, a pplication of (31) giv es D t ρ = ∂ ρ ∂ u I u t + ∂ ρ ∂ u x D x u t + ∂ ρ ∂ v I v t = ( c 1 β 2 + 2 c 2 β u + 3 c 3 u 2 ) u t + ( c 5 v + 2 c 6 u x ) u tx + (2 c 4 v + c 5 u x ) v t . (64) Use ( 5 3) t o eliminate u t , u tx , and v t . Then, E = − D t ρ ev aluates to E = ( c 1 β 2 + 2 c 2 β u + 3 c 3 u 2 ) v x + ( c 5 v + 2 c 6 u x ) v 2 x +(2 c 4 v + c 5 u x )( β u x − 3 u u x − α u 3 x ) , (65) whic h must b e exact. Th us, require that L (0) u ( x ) E = ( L (0) u ( x ) E , L (0) v ( x ) E ) ≡ (0 , 0 ) . Gr o up lik e terms. Set their co eﬃcien ts equal to zero to obtain the par ameterized system β ( c 2 − c 4 ) = 0 , c 3 + c 4 = 0 , c 5 = 0 , α c 5 = 0 , β c 5 = 0 , α c 4 − c 6 = 0 , (66) where α 6 = 0 and β 6 = 0 . In ves tigate the eliminan t of the system. Set c 1 = 1 and obtain the solution c 1 = 1 , c 2 = c 4 , c 3 = − c 4 , c 5 = 0 , c 6 = αc 4 , (67) whic h holds without condition on α and β . Substitute (67) in to (63) to get ρ = β 2 u + c 4 ( β u 2 − u 3 + v 2 + α u 2 x ) . (68) The densit y m ust b e split into independen t pieces. Indeed, since c 4 is arbitrary , set c 4 = 0 or c 4 = 1 , thus splitting (68) in to t w o indep enden t densities, ρ = β 2 u ≡ u and ρ = β u 2 − u 3 + v 2 + α u 2 x , (69) 18 whic h are ρ (1) and ρ (4) in (5 6)-(59). Step 3 : C ompute t he ﬂux J Compute the ﬂux cor r esponding to ρ in (69). Substitute (67 ) in to (65 ) . T ak e t he terms in c 4 and set c 4 = 1 . Th us, E = 2 β uv x + 2 β v u x − 3 u 2 v x − 6 u v u x + 2 α u x v 2 x − 2 αv u 3 x , ( 7 0) whic h is of order M 1 = 3 in u and o rder M 2 = 2 in v . Using (34) and (35), one readily obtains I u E = 2 β uv − 6 u 2 v + 2 αu x v x − 2 αu 2 x v , (71) and I v E = 2 β uv − 3 u 2 v + 2 αu x v x − 2 αu 2 x v . (72) Hence, using (33), J = H u ( x ) E = Z 1 0 ( I u E + I v E )[ λ u ] dλ λ = Z 1 0  4 β λuv − 9 λ 2 u 2 v + 4 αλu x v x − 4 α λu 2 x v  dλ = 2 β u v − 3 u 2 v + 2 αu x v x − 2 α u 2 x v , (73) whic h is J (4) in (59). One can set β = 1 at the end of the computations. 6 Conser v ation La ws of Syste ms of PDEs wit h T rans- cendental Nonli nearities W e no w turn to the sym b olic computation of conserv ation la ws of certain classes of PDEs with transcenden tal no nlinearities. W e o nly consider PD Es where the transcenden tal func- tions act on one dependen t v a r iable u (and not on x − deriv ativ es of u ) . In con tra st to the examples in t he previous sections, the candidate densit y will no longer ha v e c onstant unde- termined co eﬃcien ts but functional co eﬃcien ts whic h dep end o n the v ariable u. F urthermore, w e consider only PDEs whic h hav e one type of nonlinearit y . F or example, sine, or cosine, or exp o nen tial terms are ﬁne but not a mixture of these f unctions. 6.1 The sine-Gordon Equation: Dilation In v ariance and Conser- v ation L a ws The sine-Gordon (sG) equation app ears in the literat ure [17, 69] in tw o diﬀeren t w ays: • In ligh t- cone co ordinates the sG equation, u xt = sin u, has a mixed deriv a t ive term, whic h complicates matters. W e r eturn to this type of equation in Section 7.1. • The sG equation in lab oratory co ordinates, u 2 t − u 2 x = sin u, can b e recast as u t + v = 0 , v t + u 2 x + sin u = 0 , (74) 19 where v ( x, t ) is an auxiliary v ariable. System (74) is amenable to o ur approach, sub ject to mo diﬁcations to a ccommo da t e the transcenden tal nonlinearit y . The sG equation describ es t he propagation of crystal dislocations, sup erconductivit y in a Josephson junction, and ultra-short optical pulse propagation in a resonan t medium [69]. In mathematics, the sG equation is long kno wn in the diﬀeren tia l geometry of surfaces of constan t negat iv e G aussian curv ature [30, 80]. The sine-Gordon equation (74) is not uniform in rank unless w e replace it by u t + v = 0 , v t + u 2 x + α s in u = 0 , (75) where α is a real para meter with w eigh t . Indee d, substituting t he Maclaurin series, sin u = u − u 3 3! + u 5 5! − · · · , and requiring uniformity in rank yields W ( u ) + W ( ∂ /∂ t ) = W ( v ) , W ( v ) + W ( ∂ /∂ t ) = W ( u ) + 2 = W ( α ) + W ( u ) = W ( α ) + 3 W ( u ) = W ( α ) + 5 W ( u ) = · · · . (7 6 ) This forces us to set W ( u ) = 0 and W ( α ) = 2 . Consequen tly , (75) is scaling inv arian t under the symmetry ( x, t, u, v , α ) → ( λ − 1 x, λ − 1 t, λ 0 u, λ 1 v , λ 2 α ) , (77) corresp onding to W ( ∂ /∂ x ) = W ( ∂ /∂ t ) = 1 , W ( u ) = 0 , W ( v ) = 1 , W ( α ) = 2 . The ﬁrst a nd second equations in (75) are uniform of ranks 1 a nd 2, resp ectiv ely . The ﬁrst few (of inﬁnitely man y) densit y-ﬂux pairs [8, 29] for the sG equation (75) are ρ (1) = 2 α cos u + v 2 + u 2 x , J (1) = 2 v u x , (78) ρ (2) = 2 v u x , J (2) = − 2 α cos u + v 2 + u 2 x , (79) ρ (3) = 6 α v u x cos u + v 3 u x + v u 3 x − 8 v x u 2 x , (80) ρ (4) = 2 α 2 cos 2 u − 2 α 2 sin 2 u + 4 α v 2 cos u + v 4 + 20 α u 2 x cos u + 6 v 2 u 2 x + u 4 x − 16 v 2 x − 16 u 2 2 x , (81) J (3) and J (4) are not show n due to length. Again, all densities and ﬂuxe s are uniform in rank (b efore α is set equal to 1). 6.2 Computation of a Conserv ation La w for the sine-Gordon S ys- tem W e show ho w to compute densities ρ (1) and ρ (2) , b oth of rank 2, and their asso ciated ﬂuxes J (1) and J (2) . The candidate densit y will no longer ha ve c onstant undetermined co eﬃcien ts c i but functional co eﬃcien ts h i ( u ) whic h dep end on the v ariable with weigh t zero [8]. T o a void ha ving to solv e PDEs, we tacitly assume that there is only one dep enden t v ariable with w eigh t zero. As b efore, the algorithm pro ceeds in three steps: Step 1 : C onst r uct the form of t he densit y Augmen t the set of dep enden t v ariables with α (with non-zero w eigh t) and replace u by u x (since W ( u ) = 0 ) . Hence, V = { α , u x , v } . Compute R = { α , v 2 , v 2 , u 2 x , v u x , u 2 x } and remov e div ergences and equiv alent terms to get S = { α, v 2 , u 2 x , v u x } . The candidate densit y ρ = αh 1 ( u ) + h 2 ( u ) v 2 + h 3 ( u ) u 2 x + h 4 ( u ) v u x , (82) 20 with undetermined functional co eﬃcien ts h i ( u ) . Step 2 : C ompute t he functions h i ( u ) Compute D t ρ = ∂ ρ ∂ u I u t + ∂ ρ ∂ u x D x u t + ∂ ρ ∂ v I v t = ( αh ′ 1 + v 2 h ′ 2 + u 2 x h ′ 3 + v u x h ′ 4 ) u t + (2 u x h 3 + v h 4 ) u tx + (2 v h 2 + u x h 4 ) v t , (83) where h ′ i denotes dh i du . Aft er replacing u t and v t from (75), E = − D t ρ b ecomes E = ( αh ′ 1 + v 2 h ′ 2 + u 2 x h ′ 3 + v u x h ′ 4 ) v + (2 u x h 3 + v h 4 ) v x + (2 v h 2 + u x h 4 )( α sin u + u 2 x ) . (84) E mus t b e exact. Therefore, require tha t L (0) u ( x ) E ≡ 0 and L (0) v ( x ) E ≡ 0 . Set the co eﬃcien ts of lik e terms equal to zero to get a mixed linear system of algebraic equations and ODEs: h 2 ( u ) − h 3 ( u ) = 0 , h ′ 2 ( u ) = 0 , h ′ 3 ( u ) = 0 , h ′ 4 ( u ) = 0 , h ′′ 2 ( u ) = 0 , h ′′ 4 ( u ) = 0 , 2 h ′ 2 ( u ) − h ′ 3 ( u ) = 0 , 2 h ′′ 2 ( u ) − h ′′ 3 ( u ) = 0 , (85) h ′ 1 ( u ) + 2 h 2 ( u ) sin u = 0 , h ′′ 1 ( u ) + 2 h ′ 2 ( u ) sin u + 2 h 2 ( u ) cos u = 0 . Solv e the system [8] and substitute the solution h 1 ( u ) = 2 c 1 cos u + c 3 , h 2 ( u ) = h 3 ( u ) = c 1 , h 4 ( u ) = c 2 , (86) (with arbitrary constants c i ) into (82) to obtain ρ = c 1 (2 α cos u + v 2 + u 2 x ) + c 2 v u x + c 3 α. (87) Step 3 : C ompute t he ﬂux J Compute the ﬂux corresp onding to ρ in (87). Substitute (8 6 ) into ( 84), to get E = c 1 (2 u x v x + 2 v u 2 x ) + c 2 ( αu x sin u + v v x + u x u 2 x ) . (88) Since E = D x J, one m ust in tegra t e to obtain J. Using ( 2 6) and (35) one gets I u E = 2 c 1 v u x + c 2 ( αu sin u + u 2 x ) a nd I v E = 2 c 1 v u x + c 2 v 2 . Using (33), J = H u ( x ) E = Z 1 0 ( I u E + I v E )[ λ u ] dλ λ = Z 1 0  4 c 1 λv u x + c 2 ( αu sin( λu ) + λv 2 + λu 2 x )  dλ = c 1 (2 v u x ) + c 2  − α cos u + 1 2 v 2 + 1 2 u 2 x  . (89) Finally , split densit y (87) and ﬂux (89) in to indep enden t pieces (fo r c 1 and c 2 ) t o g et ρ (1) = 2 α cos u + v 2 + u 2 x , J (1) = 2 v u x , (90) ρ (2) = v u x , J (2) = − α cos u + 1 2 v 2 + 1 2 u 2 x . (91) F or E in ( 8 8), J in (8 9) can easily b e computed by hand [8]. Ho wev er, the computation o f ﬂuxes corr esp onding to densities of ranks ≥ 2 is cum b ersome and requires in tegration with the homotopy op erator. 21 7 Conser v ation La w s of Scalar Equatio ns with T rans- cendental and Mixe d Deriv ativ e T erms Our metho d to compute densities and ﬂuxes of scalar equations with transcenden tal terms and a mixed deriv ativ e term (i.e. u xt ) is an adaptat io n of the tec hnique sho wn in Section 5. W e o nly consider single PDEs with one t yp e of transcenden tal nonlinearity . Since w e are no longer dealing with ev olution equations, densities and ﬂuxes could dep enden t on u t , u 2 t , etc. W e do not co ver suc h cases; instead, w e refer the reader to [101 ]. 7.1 The sine-Gord on Equation in Ligh t-Cone Co ordinates In ligh t - cone co ordinates (or c haracteristic co ordinates) the sG equation, u xt = sin u, (92) has a mixed deriv ativ e as w ell as a transcenden tal term. A change of v ariables, Φ = u x , Ψ = − 1 + cos u, allo ws o ne to replace (92) by Φ xt − Φ − ΦΨ = 0 , 2Ψ + Ψ 2 + Φ 2 t = 0 , (93) without transcenden ta l terms. Unfortunately neither (92) nor (93) can b e written as a system of evolution equations. As sho wn in Section 6.1, to deal with the transcenden tal nonlinearity , whic h imp oses W ( u ) = 0 , one has to replace (92) by u xt = α sin u , (94) where α is an auxiliary parameter with w eight. Indeed, (94) is dilation in v ariant under the scaling symmetry ( x, t, u, α ) → ( λ − 1 x, λ − 1 t, λ 0 u, λ 2 α ) , (95) corresp onding to W ( ∂ /∂ x ) = W ( ∂ /∂ t ) = 1 , W ( u ) = 0 , and W ( α ) = 2 . The densit y-ﬂux pairs [8, 29] of ranks 2 , 4 , 6 , and 8 (whic h are indep enden t of u t , u 2 t , etc.), are ρ (1) = u 2 x , J (1) = 2 α cos u, (96) ρ (2) = u 4 x − 4 u 2 2 x J (2) = 4 αu 2 x cos u , (97) ρ (3) = u 6 x − 20 u 2 x u 2 2 x + 8 u 2 3 x , J (3) = 2 α (3 u 4 x cos u + 8 u 2 x u 2 x sin u − 4 u 2 2 x cos u ) , (98) ρ (4) = 5 u 8 x − 28 0 u 4 x u 2 2 x − 11 2 u 4 2 x + 224 u 2 x u 2 3 x − 64 u 2 4 x , (99) J (4) = 8 α  5 u 6 x cos u + 40 u 4 x u 2 x sin u + 20 u 2 x u 2 2 x cos u + 16 u 3 2 x sin u − 16 u 3 x u 3 x cos u − 48 u x u 2 x u 3 x sin u + 8 u 2 3 x cos u  . (100) There are inﬁnitely many densit y- ﬂux pairs (all of ev en rank). Since u xt = ( u x ) t , one can view (94) as an ev olution equation in a new v ariable, U = u x , and construct densities as linear com binations with constant co eﬃcien ts of monomia ls in U and its x − deriv at ives . As b efore, eac h monomial has a (pre-selected) rank. T o accommo date the transcenden tal term(s) o ne migh t b e incorrectly tempted to linearly com bine suc h monomials with functional co eﬃcien ts h i ( u ) instead of constant co eﬃcien ts c i . F or example, ho w ev er, for r a nk R = 2 , one should ta k e ρ = c 1 u 2 x instead of ρ = h 1 ( u ) u 2 x , b ecause the lat t er w ould lead to D t ρ = h ′ 1 u t u 2 x + 2 h 1 u x u 2 x and u t cannot b e replaced from (94). 22 7.2 Examples of Equations with T r ans c end en tal Nonlinearities In this section w e consider additional PDEs of the form u xt = f ( u ) , where f ( u ) has trans- cenden tal terms. Using the P ainlev ´ e in tegrabilit y test, researc hers [12] ha ve concluded that the only PDEs of that t yp e that a re completely in tegrable are equiv alen t to one of the standard forms of the nonlinear Klein-Gordon equation [2, 12]. These standard forms (in ligh t-cone co ordinates) include the sine-Gordon equation, u xt = sin u, discussed in Section 6.1, the sinh-Gordon equation u xt = sinh u, the Liouville equation u xt = e u , and the double L io uville equations, u xt = e u ± e − 2 u . The latter is also referred to in the literature as the Tzetzeica and Mikhailo v equations. F or eac h of these equations o ne can compute conserv a t io n laws with the metho d discussed in Section 7.1 . Alternat iv ely , if these equations we r e tr ansformed into lab oratory co ordinates, one w ould apply the metho d of Section 6 .2 . The multiple sine-Gordon equations, e.g. u xt = sin u + sin 2 u, hav e only a ﬁnite num b er of conserv ation law s and are not completely integrable, as supp orted b y other evidence [2]. The sinh-Gordon equation, u xt = sinh u , arises as a sp ecial case o f t he T o da lattice discusse d in Section 11.1. It also describ es the dynamics of strings in constan t curv ature space-times [70]. In thermo dynamics, the sinh-Gor don equation can b e used to calculate partition and correlation f unctions, and th us supp ort Langevin sim ula tions [64]. In T able 1, w e sho w a few densit y-ﬂux pa ir s for the sinh-Gordon equation in ligh t- cone co ordinates, u xt = α sinh u. As with the sG equation (92), W ( ∂ / ∂ x ) = W ( ∂ /∂ t ) = 1 , W ( u ) = 0 , and W ( α ) = 2 . The ranks in the ﬁrst column of T able 1 corresp ond to the ranks of the densities, whic h are p olynomial in U = u x and its x − deriv ativ es. The sinh-Gor don equation has inﬁnitely man y conserv a t io n la ws and is know n to b e completely in tegrable [2]. T able 1 : Conserv ation La ws of the sinh-Gordon equation, u xt = α sinh u Rank Densit y ( ρ ) Flux ( J ) 2 u 2 x − 2 α cosh u 4 u 4 x + 4 u 2 2 x − 4 αu 2 x cosh u 6 u 6 x + 20 u 2 x u 2 2 x + 8 u 2 3 x − 2 α [(3 u 4 x + 4 u 2 2 x ) cosh u + 8 u 2 x u 2 x sinh u ] 8 5 u 8 x + 280 u 4 x u 2 2 x − 112 u 4 2 x − 8 α [(5 u 6 x − 20 u 2 x u 2 2 x + 16 u 3 x u 3 x + 8 u 2 3 x ) cosh u +224 u 2 x u 2 3 x + 64 u 2 4 x +(40 u 4 x u 2 x − 16 u 3 2 x + 48 u x u 2 x u 3 x ) sinh u ] . The Liouville equation, u xt = e u , pla ys an imp ort a n t role in mo dern ﬁeld t heory [68], e.g. in the theory of strings, where the quantum Liouville ﬁeld app ears as a conformal anomaly [65]. The ﬁrst few (of inﬁnitely man y) densit y-ﬂux pairs for the Liouville equation in light- cone co ordinates, u xt = α e u , are giv en in T able 2. As b efore, W ( ∂ / ∂ x ) = W ( ∂ /∂ t ) = 1 , W ( u ) = 0 , and W ( α ) = 2 . The ranks in t he ta ble refer to the ranks of the densities. Do dd and Bullough [29 ] ha ve sho wn that the Liouville equation has no densities of ranks 3 , 5 , and 7 . As sho wn in T able 2, there are tw o densities of rank 6 , and three densities o f ra nk 8 . Our results agree with those in [29], where one can also ﬁnd the unique densit y o f rank 9 a nd the four indep enden t densities o f ra nk 1 0 . 23 T able 2 : Conserv ation La ws of the Liouville equation, u xt = α e u R Densit y ( ρ ) Flux ( J ) 2 u 2 x − 2 α e u 4 u 4 x + 4 u 2 2 x − 4 αu 2 x e u 6 c 1 ( u 6 x − 20 u 3 2 x − 12 u 2 3 x ) − α [6 c 1 ( u 4 x − 4 u 2 x u 2 x − 2 u 2 2 x ) + c 2 ( u 2 x u 2 2 x + u 3 2 x + u 2 3 x ) + c 2 u 2 x (2 u 2 x + u 2 x )] e u 8 c 1 ( u 8 x − 56 u 2 x u 3 2 x − 168 u 2 x u 2 3 x − α [8 c 1 ( u 6 x − 6 u 4 x u 2 x + 3 u 2 x u 2 2 x − 20 u 3 2 x − 672 u 2 x u 2 3 x − 14 4 u 2 4 x ) − 36 u 3 x u 3 x − 108 u x u 2 x u 3 x − 18 u 2 3 x ) + c 2 ( u 4 x u 2 2 x + u 2 x u 3 2 x + 5 u 2 x u 2 3 x + c 2 (2 u 4 x u 2 x − u 2 x u 2 2 x + 4 u 3 2 x + 8 u 3 x u 3 x +18 u 2 x u 2 3 x + 4 u 2 4 x ) +24 u x u 2 x u 3 x + 4 u 2 3 x ) + c 3 ( u 4 2 x + 3 u 2 x u 2 3 x + 15 u 2 x u 2 3 x + 3 u 2 4 x ) + c 3 (4 u 3 2 x + 6 u 3 x u 3 x + 18 u x u 2 x u 3 x + 3 u 2 3 x )] e u The double Liouville equations, u xt = e u ± e − 2 u , (101) arise in the ﬁeld of “laser-induced vibrational predesorption of molecules phy sisorb ed on insulating subs trates.” More precise ly , (101) is used to in v estigate the dynamics of energy ﬂo w of excited admolecules on insulating substrates [82]. Double Liouville equations ar e also relev ant in studies of global prop erties of scalar-v acuum conﬁgurations in general relativit y and similarly systems in alt ernat ive theories of gravit y [22]. In T able 3, w e show some density - ﬂux pa irs of u xt = α (e u − e − 2 u ) . There are no densit y- ﬂux pairs for ranks 4 and 10 . W e computed a densit y-ﬂux pair for rank 1 2 (not sho wn due to length). The results fo r (1 01) with the plus sign are similar. P art I I: Nonlinear Diﬀere n tial-Diﬀer e nce Equ ati ons In the second part o f this c hapter we discuss tw o distinct metho ds to construct conserv ation la ws of nonlinear DDEs. The ﬁrst metho d follows closely the tec hnique for PDEs discussed in P art I. It is quite eﬀectiv e for certain classes of DDEs, including the Kac-v an Mo erb ek e and T o da latt ices, but f ar less eﬀectiv e for more complicated lattices, suc h a s the Bogoy a vlenskii and the Gardner la ttices. The latter examples are treated with a new metho d based on a leading order analysis prop osed b y Hic kman [5 5]. 8 Nonlinear DDEs and Con serv atio n La ws W e consider autonomous nonlinear systems of DDEs of the form ˙ u n = F ( u n − l , ..., u n − 1 , u n , u n +1 , ..., u n + m ) , (102) 24 T able 3 : Conserv ation Laws of the double Liouville equation, u xt = α (e u − e − 2 u ) Rank Densit y ( ρ ) Flux ( J ) 2 u 2 x − α (2e u + e − 2 u ) 4 —— —— 6 u 6 x + 15 u 2 x u 2 2 x − 5 u 3 2 x + 3 u 2 3 x − 3 α [(2 u 4 x + 2 u 2 x u 2 x + u 2 2 x )e u +( u 4 x − 8 u 2 x u 2 x + 2 u 2 2 x )e − 2 u ] 8 u 8 x + 42 u 4 x u 2 2 x − 14 u 2 x u 3 2 x − 7 u 4 2 x − α [(8 u 6 x + 36 u 4 x u 2 x − 18 u 2 x u 2 2 x − 20 u 3 2 x +21 u 2 x u 2 3 x − 21 u 2 x u 2 3 x + 3 u 2 4 x +6 u 3 x u 3 x + 18 u x u 2 x u 3 x + 3 u 2 3 x )e u +(4 u 6 x − 72 u 4 x u 2 x − 18 u 2 x u 2 2 x − 4 u 3 2 x +48 u 3 x u 3 x − 72 u x u 2 x u 3 x + 6 u 2 3 x )e − 2 u ] 10 —— —— where u n and F are v ector-v alued functions with N comp onents. W e only consider DDEs with one discrete v ariable, denoted by in t eger n. In man y applications, n comes from a discretization of a space v a r ia ble. The dot stands for diﬀerentiation with resp ect to the con tinuous v ariable whic h frequen tly is time, t. W e assume that F is po lynomial with constan t co eﬃcien ts, although t his restriction can b e waiv ed for the metho d presen ted in Section 12. No restrictions a re imp osed on the degree of nonlinearity of F . If parameters are presen t in (102), they will b e denoted b y low er-case G reek letters. F dep ends on u n and a ﬁnite num b er of fo r w ard and bac kward shifts o f u n . W e iden tify l with the furthest negative shift of a ny v ariable in the system, a nd m with the furthest p ositiv e shift of any v ariable in the system. No restrictions are imp osed on the in tegers l and m, whic h measure the degree of non-lo cality in (102). By analogy with D x , w e deﬁne the shift op erator D by D u n = u n +1 . The op erato r D is often called the up-shift op er ator o r forw a r d- o r righ t-shift op erator. Its in v erse, D − 1 , is the down-shift op er a tor or backw ar d- or left-shift op erator, D − 1 u n = u n − 1 . The action of the shift op erators is extended to functions b y acting on their arguments . F or example, D F ( u n − l , ..., u n − 1 , u n , u n +1 , ..., u n + m ) = F ( D u n − l , ..., D u n − 1 , D u n , D u n +1 , ..., D u n + m ) = F ( u n − l +1 , ..., u n , u n +1 , u n +2 , ..., u n + m +1 ) . (103) F ollo wing [57], we g enerate (102) f r o m ˙ u 0 = F ( u − l , u − l +1 , ..., u − 1 , u 0 , u 1 , ..., u m − 1 , u m ) , (104) where ˙ u n = D n ˙ u 0 = D n F . T o further simplify the nota tion, we denote the zero-shifted dep enden t v ariable, u 0 , b y u . Shifts of u a re generated by repeated application o f D . F or instance, u k = D k u . The iden tity op erator is denoted by I , where D 0 u = I u = u . A c onservation law of (104), D t ρ + ∆ J = 0 , (105) 25 links a c onserve d den sity ρ to an a s s o ci a te d ﬂux J, where b oth are scalar functions dep ending on u and its shifts. In (1 05), whic h holds on solutions of ( 1 04), D t is the tota l deriv ativ e with resp ect to time and ∆ = D − I is the fo rw ard diﬀerence op erato r. F or readability (in particular, in the examples), the comp o nents of u will b e denoted b y u, v , w , etc. In what fo llo ws w e consider only auto nomous functions, i.e. F , ρ, and J do not explicitly dep end on t. The time deriv ativ es are deﬁned in a similar wa y as in the con tinuous case, see (5 ) and (31). W e sho w t he discrete analog of (31 ) . F o r a densit y ρ ( u p , u p +1 , · · · , u q , v r , v r +1 , · · · , v s ), in volv ing tw o dep enden t v ariables ( u, v ) and t heir shifts, the time deriv ativ e is computed a s D t ρ = q X k = p ∂ ρ ∂ u k ˙ u k + s X k = r ∂ ρ ∂ v k ˙ v k = q X k = p ∂ ρ ∂ u k D k ! ˙ u + s X k = r ∂ ρ ∂ v k D k ! ˙ v , (106) since D and d / d t comm ute. Obvious ly , the diﬀerence op erato r extends to functions. F or example, ∆ J = D J − J for a ﬂux, J. A density is trivial if there exists a function ψ so that ρ = ∆ ψ . Similar to the con tin uous case, we sa y tha t tw o densities, ρ (1) and ρ (2) , are e quivalent if and only if ρ (1) + k ρ (2) = ∆ ψ , for some ψ and some non-zero scalar k . It is paramoun t that the densit y is free of equiv alen t terms for if suc h terms w ere presen t, they could b e mov ed into the ﬂux J. Instead of w orking with diﬀeren t densities, w e will use the equiv alence of monomial terms t i in the same densit y (of a ﬁxed rank). Comp ositions of D or D − 1 deﬁne an e quivalenc e r elation ( ≡ ) on monomial terms. Simply stated, all shifted terms are e quivalent , e.g. u − 1 v 1 ≡ uv 2 ≡ u 2 v 4 ≡ u − 3 v − 1 since u − 1 v 1 = uv 2 − ∆ ( u − 1 v 1 ) = u 2 v 4 − ∆ ( u 1 v 3 + uv 2 + u − 1 v 1 ) = u − 3 v − 1 + ∆ ( u − 2 v + u − 3 v − 1 ) . ( 1 07) This equiv alence relat io n holds for a ny function of the dep enden t v ariables, but for the con- struction of conserv ed densities we will a pply it only to monomials. In the alg orithm used in Sections 9.2, 11.2, and 11.4, w e will use the follow ing e quivalenc e criterion : t w o monomial terms , t 1 and t 2 , are equiv alen t, t 1 ≡ t 2 , if and only if t 1 = D r t 2 for some in teger r. Ob viously , if t 1 ≡ t 2 , then t 1 = t 2 + ∆ J for some p olynomial J, whic h dep ends on u a nd its shifts. F or example, u − 2 u ≡ u − 1 u 1 b ecause u − 2 u = D − 1 u − 1 u 1 . Hence, u − 2 u = u − 1 u 1 + [ − u − 1 u 1 + u − 2 u ] = u − 1 u 1 + ∆ J, with J = − u − 2 u. F or eﬃciency we need a criterion to c ho ose a unique represen ta tiv e from eac h equiv a- lence class. There a r e a n um b er of wa ys to do this. W e deﬁne the c anon ic a l represen tativ e as that mem b er that has (i) no negative shifts and (ii) a non-tr ivial dep endence on the lo- c al (that is, zero-shifted) v ariable. F or example, uu 2 is t he canonical represe n tative of the class {· · · , u − 2 u, u − 1 u 1 , uu 2 , u 1 u 3 , · · · } . In t he case of e.g. tw o v a riables ( u and v ), u 2 v is the canonical represen tativ e of the class {· · · , u − 1 v − 3 , uv − 2 , u 1 v − 1 , u 2 v , u 3 v 1 , · · · } . Alternativ ely , one could c ho o se a v ariable ordering and then choo se the mem b er that dep ends on the zero-shifted v ariable of low est lexicographical order. The co de in [48 ] uses lexicographical ordering of the v ariables, i.e. u ≺ v ≺ w , etc. Th us, uv − 2 (instead of u 2 v ) is c hosen a s the canonical represen tative of {· · · , u − 1 v − 3 , uv − 2 , u 1 v − 1 , u 2 v , u 3 v 1 , · · · } . 26 It is easy to sho w [55] that if ρ is a densit y then D k ρ is also a densit y . Hence, using an appropria te “ up-shift” all negative shifts in a densit y can b een remo v ed. Th us, without loss of generalit y , w e ma y assume that a densit y that dep ends o n q shifts has c anonic al form ρ ( u , u 1 , · · · , u q ) . 9 The Metho d of Un d etermined C o e ﬃci e n ts for DDEs In this section w e sho w how p olynomial conserv ation laws can b e computed for a scalar DDE, ˙ u = F ( u − l , u − l +1 , · · · , u, · · · , u m − 1 , u m ) . (108) The Kac-v an Mo erb ek e example is used to illustrate the steps. 9.1 A Classic Example: The Kac-v an Mo erb ek e Lattice The Kac-v an Mo erb ek e (KvM) lattice [60, 62 ], also kno wn as the V olterra lattice, ˙ u n = u n ( u n +1 − u n − 1 ) , (109) arises in the study of Langmuir o scillations in plasmas, p opulation dynamics, etc. Eq. (1 0 9) app ears in the literature in v arious forms, including ˙ R n = 1 2 (exp( − R n − 1 ) − exp( − R n +1 )) , and ˙ w n = w n ( w 2 n +1 − w 2 n − 1 ) , whic h relate to (10 9) by simple tra nsformations [92]. W e contin ue with (109) and, adhering to the simpliﬁed notation, write it as ˙ u = u ( u 1 − u − 1 ) , (110) or, with the conv entions adopted a b o ve , ˙ u = u ( D u − D − 1 u ) . Lattice (110) is in v ariant under the scaling symmetry ( t, u ) → ( λ − 1 t, λu ) . Hence, u corre- sp onds to one deriv a tiv e with resp ect to t, i.e. u ∼ d dt . In analogy to the con t inuous case, w e deﬁne the weig h t W of a v ariable as t he exp onen t of λ that m ultiplies the v ariable [39, 4 0]. W e assume that shifts of a v aria ble hav e the same w eights , t ha t is, W ( u − 1 ) = W ( u ) = W ( u 1 ) . W eights of dep enden t v ariables are nonnegative and rational. The r ank of a monomial equals the total we ig h t of the mo no mial. An ex pression (or equation) is unifo rm in r ank if all its monomial terms hav e equal ra nk. Applied to (1 10), W (d / d t ) = W ( D t ) = 1 and W ( u ) = 1 . Con v ersely , the scaling symmetry can b e computed with linear algebra as follows . Setting W (d / d t ) = 1 and requiring that ( 110) is uniform in rank yields W ( u ) + 1 = 2 W ( u ) . Th us, W ( u ) = 1 , whic h agrees with the scaling symmetry . Man y in tegrable nonlinear DD Es are scaling in v arian t . If no t , they can b e made so by extending the set of dep endent v ariables with parameters with w eigh ts. Examples of suc h cases are give n in Sections 11 .3 a nd 13.2. The KvM lattice has inﬁnitely man y p olynomial density -ﬂux pairs. W e giv e the conserv ed 27 densities of rank R ≤ 4 with asso ciated ﬂuxes ( J (4) is omitted due to length): ρ (1) = u , J (1) = − uu − 1 , (111 ) ρ (2) = 1 2 u 2 + uu 1 , J (2) = − ( u − 1 u 2 − u − 1 uu 1 ) , (112) ρ (3) = 1 3 u 3 + uu 1 ( u + u 1 + u 2 ) , J (3) = − ( u − 1 u 3 + 2 u − 1 u 2 u 1 + u − 1 uu 2 1 + u − 1 uu 1 u 2 ) , (113) ρ (4) = 1 4 u 4 + u 3 u 1 + 3 2 u 2 u 2 1 + uu 2 1 ( u 1 + u 2 ) + uu 1 u 2 ( u + u 1 + u 2 + u 3 ) . (114 ) In addition to inﬁnitely man y p olynomial conserv ed densities, (110) has a non-p o lynomial densit y , ρ (0) = ln u with ﬂux J (0) = − ( u + u − 1 ) . W e discuss the computation of non-p olynomial densities in Section 12 . 9.2 The Metho d of Undetermined Co eﬃcien ts A p plied to a Scalar Nonlinear DDE W e outline how densities and ﬂuxes can b e constructed fo r a scalar DD E (108). Using (110) as an example, w e compute ρ (3) of rank R = 3 and asso ciated ﬂux J (3) of rank R = 4 , b ot h listed in (113). • Select the rank R of ρ. Star t from the set V of dep enden t v ariables (including parameters with weigh t, when applicable), and fo rm a set M of a ll non-constant monomials of rank R or less (without shifts). F or eac h monomial in M intro duce the righ t n umber of t − deriv ative s to adjust the rank of that term. Using the DDE, ev aluate the t − deriv atives , strip oﬀ the n umerical co eﬃcien ts, and gather the resulting t erms in a set R . F or the KvM lattice (1 10), V = { u } and M = { u 3 , u 2 , u } . Since u 3 , u 2 , and u ha ve r anks 3 , 2 , and 1 , respectiv ely , o ne computes d 0 u 3 d t 0 = u 3 , d u 2 d t = 2 u ˙ u = 2 u 2 ( u 1 − u − 1 ) = 2 u 2 u 1 − 2 u − 1 u 2 , (115) and d 2 u d t 2 = d ˙ u d t = d ( u ( u 1 − u − 1 )) d t = ˙ u ( u 1 − u − 1 ) + u ( ˙ u 1 − ˙ u − 1 ) = u ( u 1 − u − 1 ) 2 + u ( u 1 ( u 2 − u ) − u − 1 ( u − u − 2 )) = u u 2 1 − 2 u − 1 uu 1 + u 2 − 1 u + u u 1 u 2 − u 2 u 1 − u − 1 u 2 + u − 2 u − 1 u, (116) where (110) (and it s shifts) has b een used to remo ve the time deriv atives . Build R using the terms from the righ t hand sides of the equations in (1 15) and igno r ing n umerical co eﬃcien t s, R = { u 3 , u 2 u 1 , u − 1 u 2 , uu 2 1 , u − 1 uu 1 , u 2 − 1 u, uu 1 u 2 , u − 2 u − 1 u } . (117) • Iden tify the elemen ts in R that b elong to the same equiv alence classes, replace t hem b y their canonical represen ta t ives , and remo ve all duplicates. Call the resulting set S , whic h has the building blo cks of a candidate densit y . Con tin uing with (1 17), u − 2 u − 1 u ≡ u − 1 uu 1 ≡ uu 1 u 2 . Lik ewise, u − 1 u 2 ≡ uu 2 1 and u 2 − 1 u ≡ u 2 u 1 . Thus , S = { u 3 , u 2 u 1 , uu 2 1 , uu 1 u 2 } . 28 • F orm an a rbitrary linear combination of the elemen ts in S . This is the candidate ρ. Con- tin uing with the example, ρ = c 1 u 3 + c 2 u 2 u 1 + c 3 uu 2 1 + c 4 uu 1 u 2 . (118) • Compute D t ρ = q X k =0 ∂ ρ ∂ u k ˙ u k = q X k =0 ∂ ρ ∂ u k D k ! ˙ u, (119) where q is the highest shift in ρ. Using (118) where q = 2 , D t ρ =  ∂ ρ ∂ u I + ∂ ρ ∂ u 1 D + ∂ ρ ∂ u 2 D 2  ˙ u =  (3 c 1 u 2 + 2 c 2 uu 1 + c 3 u 2 1 + c 4 u 1 u 2 ) I + ( c 2 u 2 + 3 c 3 uu 1 + c 4 uu 2 ) D + c 4 uu 1 D 2  ˙ u. ( 120) • Ev aluate D t ρ on the D DE (108) b y replacing ˙ u b y F . Call the result E . In (110), F = u ( u 1 − u − 1 ) . The ev alua ted form of (120 ) is E = (3 c 1 u 2 + 2 c 2 uu 1 + c 3 u 2 1 + c 4 u 1 u 2 ) u ( u 1 − u − 1 ) + ( c 2 u 2 + 2 c 3 uu 1 + c 4 uu 2 ) u 1 ( u 2 − u ) +( c 4 uu 1 ) u 2 ( u 3 − u 1 ) = ( 3 c 1 − c 2 ) u 3 u 1 − 3 c 1 u − 1 u 3 + 2( c 2 − c 3 ) u 2 u 2 1 − 2 c 2 u − 1 u 2 u 1 + c 3 uu 3 1 − c 3 u − 1 uu 2 1 − c 4 u − 1 uu 1 u 2 + ( c 2 − c 4 ) u 2 u 1 u 2 + 2 c 3 uu 2 1 u 2 + c 4 uu 1 u 2 2 + c 4 uu 1 u 2 u 3 . (121) • Set J = 0 . T ransform E in to its canonical form. In doing so modify J so that E + ∆ J remains unc hanged. F or example in (121), replace − 3 c 1 u − 1 u 3 in E by − 3 c 1 uu 3 1 and add − 3 c 1 u 1 u 3 to J since uu 3 1 − [ u u 3 1 − u − 1 u 3 ] . D o the same fo r all the other terms whic h are not in canonical form. After gr o uping like terms, (121) b ecomes E = (3 c 1 − c 2 ) u 3 u 1 + ( c 3 − 3 c 1 ) uu 3 1 + 2( c 2 − c 3 ) u 2 u 2 1 + 2( c 3 − c 2 ) uu 2 1 u 2 +( c 4 − c 3 ) uu 1 u 2 2 + ( c 2 − c 4 ) u 2 u 1 u 2 , (122) with J = − (3 c 1 u − 1 u 3 + 2 c 2 u − 1 u 2 u 1 + c 3 u − 1 uu 2 1 + c 4 u − 1 uu 1 u 2 ) . (123) • E is no w the obstruction to ρ b eing a densit y . Set E = 0 and solv e for the undetermined co eﬃcien ts c i . Thus , 3 c 1 − c 2 = 0 , 3 c 1 − c 3 = 0 , c 2 − c 3 = 0 , c 3 − c 4 = 0 , c 2 − c 4 = 0 , (124) whic h yields c 2 = c 3 = c 4 = 3 c 1 , where c 1 is arbitrary . • Substitute the solution for the c i in to the candidates for ρ and J to obta in the ﬁnal densit y and asso ciated ﬂux (up to a common arbitrary factor whic h can b e set to 1 or an y o ther nonzero v a lue). F or the example, setting c 1 = 1 3 and substituting c 2 = c 3 = c 4 = 1 into (118) and (123) yields ρ (3) and J (3) as giv en in (113 ). 29 10 Discre t e Eule r and Ho motop y Op erators F or simplicit y , we will consider the case of only one discrete (dep enden t) v ariable u. F irst, w e remo ve negative shifts f r o m E . Thus , ˜ E = D l E where l corresp onds to the low est shift in E . The discrete v ariational deriv ative ( discrete Euler oper at or ) An expres sion E is exact if and only if it is a total diﬀerence. The following exactness test is w ell-know n [10, 57 ]: A necessary and suﬃcien t condition for a function E , with p ositive shifts, to b e exact is that L (0) u E ≡ 0 , where L (0) u is the discr e te variational derivative (discrete Euler op erator of order zero) [10] deﬁned b y L (0) u E = ∂ ∂ u m X k =0 D − k ! E = ∂ ∂ u  I + D − 1 + D − 2 + D − 3 + · · · + D − m  E , = I ∂ E ∂ u + D − 1 ∂ E ∂ u 1 + D − 2 ∂ E ∂ u 2 + D − 3 ∂ E ∂ u 3 + · · · + D − m ∂ E ∂ u m , (125 ) where m is highest shift o ccurring in E . Application . W e return to (121) where l = 1. Therefore, ˜ E = D E = (3 c 1 − c 2 ) u 3 1 u 2 − 3 c 1 uu 3 1 + 2( c 2 − c 3 ) u 2 1 u 2 2 − 2 c 2 uu 2 1 u 2 + c 3 u 1 u 3 2 − c 3 uu 1 u 2 2 − c 4 uu 1 u 2 u 3 + ( c 2 − c 4 ) u 2 1 u 2 u 3 + 2 c 3 u 1 u 2 2 u 3 + c 4 u 1 u 2 u 2 3 + c 4 u 1 u 2 u 3 u 4 , (126) whic h is free of negativ e shifts. Applying (12 5) to (126), where m = 4 , giv es L (0) u ˜ E = ∂ ∂ u  I + D − 1 + D − 2 + D − 3 + D − 4  ˜ E = 3 (3 c 1 − c 2 ) u 2 u 1 + 3( c 3 − 3 c 1 ) u − 1 u 2 + 2( c 2 − c 4 ) uu 1 u 2 + 4( c 2 − c 3 ) uu 2 1 +4( c 3 − c 2 ) u − 1 uu 1 + 2( c 3 − c 2 ) u 2 1 u 2 + ( c 3 − 3 c 1 ) u 3 1 + ( c 4 − c 3 ) u − 1 u 2 1 +( c 4 − c 3 ) u 1 u 2 2 + (3 c 1 − c 2 ) u 3 − 1 + ( c 2 − c 4 ) u 2 − 1 u 1 + 4( c 2 − c 3 ) u 2 − 1 u +2( c 3 − c 2 ) u − 2 u 2 − 1 + 2( c 4 − c 3 ) u − 2 u − 1 u + ( c 2 − c 4 ) u 2 − 2 u − 1 , (127) whic h, as exp ected, v anishes iden t ically when c 1 = 1 3 , c 2 = c 3 = c 4 = 1 . Due to the larg e amoun t of terms generated b y the Euler op erator, this metho d for ﬁnding the undetermined co eﬃcien ts is m uc h less eﬃcie nt than the “splitting and shifting” approach illustrated on the same example in Section 9.2. The discrete homotop y op erator As in the contin uous case, the discrete homotop y op erator reduces the in v ersion of the dif- ference op erator, ∆ = D − I , to a problem of single-v ariable calculus. Indeed, assuming that E is exact and free of negativ e shifts, the ﬂux J = ∆ − 1 E can b e computed without “summa- tion by pa rts.” Instead, o ne computes a single in tegral with resp ect to an auxiliary v aria ble denoted by λ (not to b e confused with λ in Section 9.1). Consider an exact expression E (of one v a riable u ) , free of negativ e shifts, and with highest shift m. The discr ete homotopy op er ator is deﬁned [61, 71, 72 ] by H u E = Z 1 0 ( I u E ) [ λu ] dλ λ , (128) 30 with I u E = m X k =1 m − k X i =0 u i ∂ ∂ u i ! D − k E = m X k =1 k X i =1 D − i ! u k ∂ E ∂ u k , (1 29) where ( I u E )[ λu ] means t hat in I u E o ne replaces u → λu, u 1 → λu 1 , u 2 → λu 2 , etc . The form ulas in (129) are equiv alen t to the one in [52], which in turn is equiv alen t to t he form ula in terms of discrete higher Euler op era t ors [51, 54]. Giv en an exact function E without negativ e shifts one has J = ∆ − 1 E = H u E . A pro of can b e found in [61, 72]. Application. Up on substitution o f c 1 = 1 3 , c 2 = c 3 = c 4 = 1 into (126), we obta in ˜ E = D E = − uu 3 1 − 2 uu 2 1 u 2 + u 1 u 3 2 − uu 1 u 2 2 − uu 1 u 2 u 3 + 2 u 1 u 2 2 u 3 + u 1 u 2 u 2 3 + u 1 u 2 u 3 u 4 , (130) where the highest shift is m = 4 . Using (1 29), I u ˜ E = 4 X k =1 k X i =1 D − i ! u k ∂ ˜ E ∂ u k = ( D − 1 ) u 1 ∂ ˜ E ∂ u 1 + ( D − 1 + D − 2 ) u 2 ∂ ˜ E ∂ u 2 +( D − 1 + D − 2 + D − 3 ) u 3 ∂ ˜ E ∂ u 3 + ( D − 1 + D − 2 + D − 3 + D − 4 ) u 4 ∂ ˜ E ∂ u 4 = D − 1 u 1 ( − 3 uu 2 1 − 4 u u 1 u 2 + u 3 2 − u u 2 2 − u u 2 u 3 + 2 u 2 2 u 3 + u 2 u 2 3 + u 2 u 3 u 4 ) +( D − 1 + D − 2 ) u 2 ( − 2 uu 2 1 + 3 u 1 u 2 2 − 2 u u 1 u 2 − u u 1 u 3 + 4 u 1 u 2 u 3 + u 1 u 2 3 + u 1 u 3 u 4 ) +( D − 1 + D − 2 + D − 3 ) u 3 ( − uu 1 u 2 + 2 u 1 u 2 2 + 2 u 1 u 2 u 3 + u 1 u 2 u 4 ) +( D − 1 + D − 2 + D − 3 + D − 4 ) u 4 ( u 1 u 2 u 3 ) = 4 ( uu 3 1 + 2 uu 2 1 u 2 + uu 1 u 2 2 + uu 1 u 2 u 3 ) , (131) whic h has the correct terms of J (3) but incorrect co eﬃcien ts. Finally , using (128 ) ˜ J = H u ( − ˜ E ) = − Z 1 0 ( I u ˜ E )[ λu ] dλ λ = − 4 Z 1 0  uu 3 1 + 2 uu 2 1 u 2 + uu 1 u 2 2 + uu 1 u 2 u 3  λ 3 dλ = − ( u u 3 1 + 2 uu 2 1 u 2 + uu 1 u 2 2 + uu 1 u 2 u 3 ) . (132) Recall that ˜ E = D E . Hence, J = D − 1 ˜ J = − ( u − 1 u 3 + 2 u − 1 u 2 u 1 + u − 1 uu 2 1 + u − 1 uu 1 u 2 ) , (133) whic h corresp o nds to J (3) in (11). The ho mot op y metho d is computationally ineﬃcien t. Ev en for a simple example, lik e (131), the integrand has a la r g e num b er of terms, mo st of whic h ev en t ua lly cancel. T o compute the ﬂux we recommend the “splitting and shifting” approa c h whic h w as illustrated (on the same example) in Section 9.2. The generalization of the exactness test to an expression E with multiple dep enden t v aria bles ( u, v , · · · ) is straigh t f orw ard. F or example, an expression E of discrete v ariables u, v and their forward shifts will b e exact if and o nly if L (0) u E = ( L (0) u E , L (0) v E ) ≡ (0 , 0) , where L (0) v is deﬁned analogously to (125). Similar to the con t in uous case, the homo t op y op erator form ulas (128) and (129) straightforw ardly g eneralize to m ultiple dependent v ariables. The reader is referred to [5 1, 5 2, 5 4 ] fo r details. 31 11 Conse rv ation La w s of Nonlinear Syst e ms of DDEs W e use the metho d discussed in Section 9.2 to compute conserv atio n la ws for the T o da and Ablo witz-Ladik la t tices. Using the latter lattice, w e illustrate a “divide and conquer” strategy , based on multiple scales, whic h allows on to circum ve n t diﬃculties in the computation o f densities of DD Es that a re not dilation in v ariant. 11.1 The T o da Lattice One of the earliest and most famous examples o f completely in tegrable DDEs is the T o da lattice [93, 94], ¨ y n = exp ( y n − 1 − y n ) − exp ( y n − y n +1 ) , (134) where y n is the displacemen t from equilibrium of the n th particle with unit mass under an exp onential decay ing in teraction force b etw een nearest neighbors. With the c hang e of v aria bles, u n = ˙ y n , v n = exp ( y n − y n +1 ) , lattice (134 ) can b e written in algebraic form ˙ u n = v n − 1 − v n , ˙ v n = v n ( u n − u n +1 ) . (135) Adhering to the simpliﬁed no tation, w e con tinue with ˙ u = v − 1 − v , ˙ v = v ( u − u 1 ) . (136) As b efore, w e set W (d / d t ) = 1 , assign unknow n w eigh ts, W ( u ) , W ( v ) , to the dep enden t v aria bles, and require that eac h equation in (136) is uniform in rank. This yields W ( u ) + 1 = W ( v ) , W ( v ) + 1 = W ( u ) + W ( v ) . (137) The solution W ( u ) = 1 , W ( v ) = 2 rev eals that (136 ) is in v ariant under the scaling symmetry ( t, u, v ) → ( λ − 1 t, λu, λ 2 v ) , (138) where λ is an arbitra ry parameter. The T o da lattice has inﬁnitely many conserv ation la ws [44]. The ﬁrst t w o density - ﬂux pairs are easy to compute by hand. Here w e giv e the densities of rank R ≤ 4 with asso ciated ﬂuxes, J (4) b eing omitted due to length: ρ (1) = u , J (1) = v − 1 , (139) ρ (2) = 1 2 u 2 + v , J (2) = uv − 1 , (140) ρ (3) = 1 3 u 3 + u ( v − 1 + v ) , J (3) = u − 1 uv − 1 + v 2 − 1 , (141) ρ (4) = 1 4 u 4 + u 2 ( v − 1 + v ) + uu 1 v + 1 2 v 2 + v v 1 . (142) 32 11.2 Computation of a Conserv ation L a w of the T o da Lattice As an example, we compute densit y ρ (3) (of rank R = 3) and asso ciated ﬂux J (3) (of rank 4) in (141). Step 1: C onst r uct the form of the densit y Start from V = { u, v } , i.e. the set of dep enden t v ariables with w eight. List all monomials in u and v of rank R = 3 or less: M = { u 3 , u 2 , uv , u , v } . Next, for eac h mono mia l in M , in tro duce the correct n umber of t -deriv ativ es so that each term has ra nk 3 . Using (136) , compute d 0 u 3 d t 0 = u 3 , d 0 uv d t 0 = uv , d u 2 d t = 2 u ˙ u = 2 uv − 1 − 2 uv , d v d t = ˙ v = uv − u 1 v , d 2 u d t 2 = d ˙ u d t = d( v − 1 − v ) d t = u − 1 v − 1 − u v − 1 − u v + u 1 v . (143) Gather the terms in t he righ t hand sides in (143) to get R = { u 3 , uv − 1 , uv , u − 1 v − 1 , u 1 v } . Iden tify mem b ers b elonging to the same equiv a lence classes and replace them by their canonical repres en tative s. F or example, uv − 1 ≡ u 1 v . Adhering to lexicographical ordering, w e will use uv − 1 instead of u 1 v . D oing so, replace R by S = { u 3 , uv − 1 , uv } , whic h has the building blo cks of the densit y . Linearly com bine the monomials in S with undetermined co eﬃcien ts c i to get the candidate densit y o f ra nk 3 : ρ = c 1 u 3 + c 2 uv − 1 + c 3 uv . (144) Step 2: C ompute t he undetermined co eﬃcien ts c i Compute D t ρ and use (136) to eliminate ˙ u a nd ˙ v and their shifts. Th us, E = (3 c 1 − c 2 ) u 2 v − 1 + ( c 3 − 3 c 1 ) u 2 v + ( c 3 − c 2 ) v − 1 v + c 2 u − 1 uv − 1 + c 2 v 2 − 1 − c 3 uu 1 v − c 3 v 2 . (145) Step 3: Find the asso ciated ﬂux J T ransform (145) in to canonical for m to obtain E = ( 3 c 1 − c 2 ) u 2 1 v + ( c 3 − 3 c 1 ) u 2 v + ( c 3 − c 2 ) v v 1 + c 2 uu 1 v + c 2 v 2 − c 3 uu 1 v − c 3 v 2 (146) with J = (3 c 1 − c 2 ) u 2 v − 1 + ( c 3 − c 2 ) v − 1 v + c 2 u − 1 uv − 1 + c 2 v 2 − 1 . (147) Set E = 0 to get the linear system 3 c 1 − c 2 = 0 , c 3 − 3 c 1 = 0 , c 2 − c 3 = 0 . (148) Set c 1 = 1 3 and substitute the solution c 1 = 1 3 , c 2 = c 3 = 1 , into (144) and (147) to obtain ρ (3) and J (3) in (1 41). 33 11.3 The A blo witz-Ladik Lattice In [4, 5, 6], Ablowitz and Ladik derived and studied the fo llo wing completely integrable discretization of the nonlinear Sc hr¨ odinger equation: i ˙ u n = u n +1 − 2 u n + u n − 1 ± u ∗ n u n ( u n +1 + u n − 1 ) , ( 1 49) where u ∗ n is the complex conjugate of u n . W e con tinue with (149) with the plus sign; the case with the negativ e sign is analogo us. Instead of splitting u n in to its real and imaginary parts, w e treat u n and v n = u ∗ n as independen t v ariables and augmen t (149) with its complex conjugate equation, ˙ u n = ( u n +1 − 2 u n + u n − 1 ) + u n v n ( u n +1 + u n − 1 ) , ˙ v n = − ( v n +1 − 2 v n + v n − 1 ) − u n v n ( v n +1 + v n − 1 ) , (150) where i has b een absorb ed into a scale on t. Since v n = u ∗ n w e ha ve W ( v n ) = W ( u n ) . Neither equation in (150) is dilation inv arian t. T o circum v en t this problem w e introduce a n auxiliary parameter α with we igh t, and replace (150) b y ˙ u = α ( u 1 − 2 u + u − 1 ) + uv ( u 1 + u − 1 ) , ˙ v = − α ( v 1 − 2 v + v − 1 ) − uv ( v 1 + v − 1 ) , (151) presen ted in the simpliﬁed nota t io n. Both equations in (151) are unif o rm in rank pro vided W ( u ) + 1 = W ( α ) + W ( u ) = 2 W ( u ) + W ( v ) , W ( v ) + 1 = W ( α ) + W ( v ) = 2 W ( v ) + W ( u ) , (152) whic h holds when W ( u ) + W ( v ) = W ( α ) = 1 . Since W ( u ) = W ( v ) w e ha ve W ( u ) = W ( v ) = 1 2 , and W ( α ) = 1 , whic h expresses that (151) is inv arian t under the scaling symmetry ( t, u, v , α ) → ( λ − 1 t, λ 1 2 u, λ 1 2 v , λα ) . (153) W e giv e t he conserv ed densities of (151) of ra nks 2 t hro ugh 4 . Only the ﬂux of rank 3 asso ciated to ρ (1) is sho wn. The o thers a re omitted due to length. ρ (1) = α ( c 1 uv − 1 + c 2 uv 1 ) , (154) J (1) = − α ( c 1 ( αuv − 2 − αu − 1 v − 1 + u − 1 uv − 2 v − 1 ) + c 2 ( αuv − α u − 1 v 1 − u − 1 uv v 1 )) , (155) ρ (2) = α  c 1 ( 1 2 u 2 v 2 − 1 + uu 1 v − 1 v + αuv − 2 ) + c 2 ( 1 2 u 2 v 2 1 + uu 1 v 1 v 2 + α uv 2 )  , (156) ρ (3) = α  c 1  1 3 u 3 v 3 − 1 + uu 1 v − 1 v ( uv − 1 + u 1 v + u 2 v 1 ) + αuv − 1 ( uv − 2 + u 1 v − 1 ) + αuv ( u 1 v − 2 + u 2 v − 1 ) + α 2 uv − 3  + c 2  1 3 u 3 v 3 1 + uu 1 v 1 v 2 ( uv 1 + u 1 v 2 + u 2 v 3 ) + αuv 2 ( uv 1 + u 1 v 2 ) + αuv 3 ( u 1 v 1 + u 2 v 2 ) + α 2 uv 3  , (157) where c 1 and c 2 are ar bit r a ry constan ts. Our results conﬁrm those in [6]. Since our metho d is restricted to p olynomial densities w e cannot compute the densit y with a logarithmic term, ρ n = ( u ∗ n ( u n − 1 + u n +1 ) − 2 ln(1 + u n u ∗ n )) , (158) whic h corresp o nds [3, 6] to the Ha milto nian of (149), viz. H = − i P n ρ n . 34 11.4 Computation of a Conserv ation L a w of the Ablo witz-Ladik Lattice T o make (150) scaling inv arian t we had to intro duce a n auxiliary pa r ameter α . This compli- cates matters in tw o w a ys as we will sho w in this section. First, to compute rather simple conserv ed densities, lik e ρ = uv − 1 and ρ = uv 1 , we will hav e to select r a nk R = 2 fo r whic h the candidate density has t wen ty terms. Ho w ev er, eigh teen of these terms ev en t ua lly drop out. Second, conserv ed densities correspo nding to low er ranks migh t app ear in the result. These lo wer-rank densities are easy to recognize for t hey are multiplied with arbitrary co eﬃcien ts c i . Consequen tly , when par a meters with w eight are in tro duced, the densities corresp onding to distinct ranks are no long er linearly indep enden t. W e compute ρ (1) of rank R = 2 in (154) with asso ciated ﬂux J (1) in (1 55). Note that ρ (1) and J (1) cannot b e computed with the steps b elo w when R = 1 . Step 1 : C onst r uct the form of t he densit y Augmen t the set of dep enden t v ar iables with the parameter α (with non-zero w eigh t) . Hence, V = { u , v , α } . Construct M = { u, v , u 2 , uv , v 2 , αu, u 3 , αv , u 2 v , uv 2 , v 3 , αu 2 , u 4 , αuv , u 3 v , α v 2 , u 2 v 2 , uv 3 , v 4 } , (159) whic h con tains all non-constan t mono mials of (c hosen) rank 2 or less (without shifts). As with the previous examples, for eac h elemen t in M add the righ t n um b er of t − deriv ativ es. Use (1 5 1) to ev aluate the t − deriv atives , gather the terms in the right hand sides, intro duce the cano nical represen tatives (based o n lexicographical ordering), a nd remo ve duplicates to get S = { α u 2 , u 4 , αuu 1 , αuv − 1 , αuv , u 3 v , u 2 u 1 v , u 2 v − 1 v , α v 2 , u 2 v 2 , uu 1 v 2 , uv − 1 v 2 , uv 3 , v 4 , αuv 1 , uu 2 1 v 1 , αv v 1 , u 2 v v 1 , uv 2 v 1 , uu 1 v 2 1 } . (160) Linearly com bine the monomials in S with undetermined co eﬃcien ts c i to get t he candidate densit y of rank 2: ρ = c 1 αu 2 + c 2 u 4 + c 3 αuu 1 + c 4 αuv − 1 + c 5 αuv + c 6 u 3 v + c 7 u 2 u 1 v + c 8 u 2 v − 1 v + c 9 αv 2 + c 10 u 2 v 2 + c 11 uu 1 v 2 + c 12 uv − 1 v 2 + c 13 uv 3 + c 14 v 4 + c 15 αuv 1 + c 16 uu 2 1 v 1 + c 17 αv v 1 + c 18 u 2 v v 1 + + c 19 uv 2 v 1 + c 20 uu 1 v 2 1 . (161) Step 2: C ompute t he undetermined co eﬃcien ts c i The computations pro ceed as in the examples in Sections 9.2 and 1 1.2. Th us, compute D t ρ a nd use (151) to eliminate ˙ u and ˙ v and their shifts. Ne xt, bring the expression E in t o canonical form to obta in the linear system for the undetermined co eﬃcien ts c i . Without showing the lengthy computations, one ﬁnds tha t all constants c i = 0 , except c 4 and c 15 whic h are arbitrary . Substitute the co eﬃcien ts into (16 1) to get ρ (1) in (154). Step 3 : Find the asso ciated ﬂux J The asso ciated ﬂux comes for free when E is transformed into canonical fo r m. Alternative ly , one could apply the homotop y approac h for m ultiple dependent v ariables [51, 54] to compute the ﬂux. In either case, one g ets J (1) in (155). 35 11.5 A “Divide and Conquer” Strategy It should b e clear from the example in Section 11.4 that our metho d is not eﬃcien t if the densities are not o f the form (111 )-(114) f o r the KvM lattice and (139)-(142 ) for the T o da lattice. Indeed, the densities for t he AL lattice in (154) - (157) are quite diﬀeren t in structure. Therefore, in [34], Eklund presen ted alternate strategies t o deal more eﬃcien tly with D DEs, in particular with tho se inv olving we ig h ted parameters suc h as (15 1). A ﬁrst alternativ e is to w ork with multiple sc ales by setting either W ( D t ) = 0 or W ( D t ) = 1 , the latter c hoice is what w e ha v e used th us far. A second p ossibilit y is to leav e W ( D t ) unsp eciﬁed and, if needed, intro duce extra parameters with w eight in to the D DE. Let us explore these ideas for the AL lattice (15 1), which we therefore replace b y ˙ u = α β ( u 1 − 2 u + u − 1 ) + β uv ( u 1 + u − 1 ) , ˙ v = − α β ( v 1 − 2 v + v − 1 ) − β uv ( v 1 + v − 1 ) , (162) where β is a second auxiliary parameter with we igh t. Requiring uniformity in rank leads to W ( u ) + W ( D t ) = W ( α ) + W ( β ) + W ( u ) = W ( β ) + 2 W ( u ) + W ( v ) , W ( v ) + W ( D t ) = W ( α ) + W ( β ) + W ( v ) = W ( β ) + 2 W ( v ) + W ( u ) , (163) whic h are satisﬁed when W ( α ) = W ( u ) + W ( v ) and W ( D t ) = W ( u )+ W ( v )+ W ( β ) . Therefore, w e can set W ( D t ) = a, W ( u ) = b, and W ( β ) = c with a, b, c rational n um b ers so that W ( v ) = a − b − c and W ( α ) = a − c are strictly p ositive . Thus (162) is dilation in v arian t under a three-parameter fa mily o f scaling symmetries, ( t, u, v , α , β ) → ( λ − a t, λ b u, λ a − b − c v , λ a − c α, λ c β ) . (164) Scaling symme try (153) corresp onds to the case where a = 1 , b = 1 2 , and c = 0 , more precisely , β = 1 . The fact that (162) is inv ariant under m ultiple scales is adv an tageous. Indeed, one can use the inv ariance under one scale to construct a candidate densit y and, subsequen tly , use additiona l scale(s) to split ρ into smaller densities. This “divide a nd conquer” strategy drastically reduces the complexit y of the computations. The use of m ultiple scales has pro v en to b e successful in the computation of conserv ation la ws for PDEs with more than one space v aria ble [51]. i Rank Candidate ρ i Final ρ i Final J i 1 a + 2 b − c c 1 αu 2 + c 3 αuu 1 + c 6 u 3 v + c 7 u 2 u 1 v + c 16 uu 2 1 v 1 0 0 2 4 b c 2 u 4 0 0 3 2( a − c ) c 4 αuv − 1 + c 5 αuv + c 8 u 2 v − 1 v + c 10 u 2 v 2 c 4 αuv − 1 J (1) in (155 ) + c 11 uu 1 v 2 + c 15 αuv 1 + c 18 u 2 v v 1 + c 20 uu 1 v 2 1 + c 15 αuv 1 4 3 a − 2 b − 3 c c 9 αv 2 + c 12 uv − 1 v 2 + c 13 uv 3 + c 17 αv v 1 + c 19 uv 2 v 1 0 0 5 4( a − b − c ) c 14 v 4 0 0 T able 4 : Candidate densities for the Ablowitz-Ladik la t t ice (162). 36 Candidate density (161) is uniform of r ank 2 under (153) but can be split in to smaller pieces, ρ i , using (164 ), ev en without sp ecifying v alues for a, b, a nd c. Indeed, as sho wn in T able 4, ρ in (1 6 1) can b e split in to ρ 1 through ρ 5 of distinct ranks under (16 4). Steps 2 and 3 of the algorithm are then carried out fo r eac h of t hese ρ i , ( i = 1 , · · · , 5) separately . As sho wn in the table, all but one densit y lead to a trivial result. The longest densit y , ρ 3 , with 8 terms, leads to ρ (1) and J (1) in (154) and ( 1 55), resp ectiv ely . 12 A New Metho d t o Compute Conse r v ation La ws of Nonlinear DDEs In the con tinuous case, the total deriv ative D x has we ig h t one. Consequen tly , a n y densit y is b ounded with respect to the order in x. F or example, as shown in Section 3.2 for the KdV case, the candidat e densit y ρ of r ank 6 in (18) is of ﬁrst order (after remo ving the sec ond and fourth order terms.) Obviously , a densit y of rank 6 could nev er hav e leading order terms of ﬁfth or higher order in x b ecause t he rank of suc h terms w ould exceed 6 . In the discrete case, ho w eve r, the shift op erator D has no w eight. Th us, any shift D k u = u k of a dep enden t v ariable u has the same w eight a s tha t dep enden t v ariable. Simply put, W ( u k ) = W ( u ) for any in t eger k . Consequen tly , using the tota l deriv a tiv e D t as a t o ol to construct a (candidate) density has a ma jor shortcoming: the densit y ma y la ck t erms in volv ing suﬃcien tly high shifts of the v ariables. As sho wn in Section 9.2 for the KvM lattice, the candidate density ρ of rank 3 in (118) has leading order term u u 1 u 2 , i.e. the term with the highest shift (2 in this example). It is a priori not excluded that ρ mig ht hav e terms in volving higher shifts. F or example, uu 1 u 3 has rank 3 and so do inﬁnitely man y other cubic terms. Note that for this example w e constructed ρ starting from p ow ers of u, viz. u 3 , u 2 , u ; and, by rep eated diﬀeren tiation, w ork ed our wa y “down” tow ar ds the leading order term uu 1 u 2 . In this section w e outline the k ey features of a new metho d which go es in the opp osite direction: (i) ﬁrst compute the leading order t erm and subsequen tly (ii) compute the terms in volving lo w er shifts. In step (ii) only the necessary terms are computed, nothing more, nothing less. This metho d is fast and p ow erful for it circum v ents the use of the dilation in v ariance and the metho d of undetermined co eﬃcien ts. More imp ortantly , the new metho d is not restricted to densities and ﬂuxes of p o lynomial form. 12.1 Leading O r d er Analysis Consider a densit y , ρ, that dep ends on q shifts. Since D q ρ is a lso a densit y , w e may , without loss of generalit y , a ssume that ρ ha s canonical fo rm ρ ( u , u 1 , . . . , u q ) with ∂ 2 ρ ∂ u ∂ u q 6 = 0 . In [55], Hic kman deriv ed necessary conditions o n this leading term (whic h, in the system case, is a matrix). First all terms in the candidate density ρ t hat con tribute directly to the ﬂux are remo ve d. Rat her than applying the Euler op erator on the r emaining terms in ρ , the nece ssary condition [57], ∂ 2 g ∂ u ∂ u q = 0 , (165) for g to b e a total diﬀerenc e is applied to obtain a sy stem of equations for the terms that dep end on the maximal shift, u q , in ρ. This system is rewritten as a matrix equation. Solutions 37 to this system will give us the form of the leading term in ρ. W e a pply a splitting of the identit y op erator, I = ( D − I + I ) D − 1 = ∆ D − 1 + D − 1 , (166) to the part, ρ ∗ , o f the candidat e ρ that is indep enden t of the v ariables with the low est order shift. The ﬁrst term ∆ D − 1 ρ ∗ con tributes to the ﬂux while the second term D − 1 ρ ∗ has a strictly low er shift than ρ ∗ . Applying this split rep eatedly w e get I = ( D k − I + I ) D − k = ∆  D k − 1 + D k − 2 + · · · + D + I  D − k + D − k , (167) where, again, the ﬁrst term con t ributes to the ﬂux and the remainder has strictly lo w er shift. This decomp osition is rep eatedly applied to terms that do not in v olv e the low est order shifted v ariables. An y terms that remain will in v olve the low est order shifted v ariable. These terms yield the constraints on the undetermined co eﬃcien ts or unkno wn functions in the densit y ρ. As sho wn in [55], the r esult o f this split is that ρ is a densit y if a nd only if σ = D l  F ∂ ∂ u  l X j = 0 D j ρ + q X j = l +1 D j  F ∂ ∂ u  ρ (168) is a total diﬀerence, where the op erator F ∂ ∂ u = X α F α ∂ ∂ u α (169) in volv es a summation ov er the comp onen t s in the system of DDEs. As b efore, in the examples w e will use u, v , w , etc. to denote t he dep enden t v aria bles u α . No w, σ dep ends on the shifted v ariables u , . . . , u q + L where L = ma x ( l , m ). F or q > L , (165) giv es ∂ 2 σ ∂ u ∂ u q + L = D L  ∂ 2 ρ ∂ u ∂ u q ∂ F ∂ u − L  + D q  ∂ F ∂ u L  T ∂ 2 ρ ∂ u ∂ u q = 0 , (170) where T stands f or transp ose. The system case is treated in detail in [55]. F or brevity , we con tinue with t he scalar case. Let λ = ∂ F ∂ u − L , µ = ∂ F ∂ u L , (171) then the leading term will satisfy S ∂ 2 ρ ∂ u ∂ u q = 0 , (172) with S = D L λ D L + D q µ. W e immediately see that if l 6 = m then either λ or µ is zero and (172) has no non-trivial solutions. Let q = pL + r with p and r in tegers, 0 ≤ r < L, and c = p − 1 Y k =1 D k L λ ! D pL ζ . (173) 38 Then S c = p − 1 Y k =1 D k L λ D pL  λ D L ζ + ζ D r µ  . (17 4) Th us c 6 = 0 lies in the kernel o f S if and o nly if λ D L ζ + ζ D r µ = 0 (175) has a non-zero solution ζ . Supp ose (175 ) has t wo non-zero solutions, sa y , ζ 1 and ζ 2 . Then, ζ 2 ζ 1 = D L  ζ 2 ζ 1  . (176) So, since L 6 = 0, ζ 2 = aζ 1 for some constant a . Therefore, the k ernel of S is, at m ost, one dimensional . This implies that a scalar DDE can hav e, at most, one conserv ed densit y that dep ends on q shifts for q > L. The leading term, ˜ ρ, will satisfy ∂ 2 ρ ∂ u ∂ u q = c, (177) whic h, up on in tegration, yields ˜ ρ = Z Z c du du q . (178) The densit y (if it exists ) ma y no w b e computed b y a “split and shift” strat egy on this leading term. Start ing with ρ = ˜ ρ. the ob jectiv e is to succes siv ely compute the terms (of low er shift!) that m ust b e added to ρ un til D t ρ ≡ 0 . F irst, D t ρ is computed and ev aluated on the DDE. Next, all terms are shifted so that the resulting expression dep ends on u (and not on lo wer shifts of u ) . Then the leading t erms, ξ , in D t ρ are isolated. Last, we solv e D t ρ (1) = ξ + terms of low er shift . (1 79) If (17 9) has no solution then a densit y with q shifts do es no t exist. On the o t her hand, if (179) has a solutio n, the “corr ection” term ρ (1) is then subtracted from ρ and we recompute D t ρ. By construction, the highest shift t ha t o ccurs in the result will no w b e lo w er than b efore and w e rep eat the en tire pro cedure to obta in a new correction term ρ (2) . After a ﬁnite n umber of steps, w e will either ﬁnd an that the correction term do es not exist (and so the density do es not exist) or w e will obtain D t ρ ≡ 0 and ρ will b e a densit y . This algo rithm has b een co ded [56] in Maple . F urther details and w o r ked examples will b e presen ted in [58]. F or no w, w e illustrate the algor ithm with a simple example. 12.2 A“Mo diﬁed” V olterra Lattice Consider the DD E, ˙ u = u 2 ( u 2 − u − 2 ) , (180) whic h is related to the w ell-known mo diﬁed V olterra L a ttice [11]. Here L = 2 and, using (171), λ = ∂ F ∂ u − 2 = − u 2 , µ = ∂ F ∂ u 2 = u 2 . (181) 39 Th us, the condition (175) for a non- t rivial densit y b ecomes ζ u 2 r = u 2 D 2 ζ for r = 0 , 1 . F or r = 0 , w e ha ve ζ = D 2 ζ and so we may c ho ose ζ = 1 . F or r = 1 , w e ha v e ζ u 2 1 = u 2 D 2 ζ , whic h has no non-zero solutions. Th us, densities only exist for r = 0 . Since q = pL = 2 p, with p in teger, w e conclude that q m ust b e ev en. In these cases the leading term will satisfy ∂ 2 ρ ∂ u ∂ u q = c = p − 1 Y k =1 u 2 k ! = u u 2 · · · u q − 2 . (182) Therefore, after scaling to remo ve constants , the leading term is ˜ ρ = u 2 u 2 · · · u q − 2 u q . F or example, let us compute the densit y tha t depends on q = 4 shifts. W e set ρ = ˜ ρ = u 2 u 2 u 4 . Using (180), D t ρ = u 2 u 2 2 u 4 ( u 2 − u − 2 ) + 2 u u 3 2 u 4 ( u 4 − u ) + u u 2 2 u 2 4 ( u 6 − u 2 ) ≡ u u 3 2 u 2 4 − u 2 u 3 2 u 4 . (183) The highest shift is u 4 and so the leading terms a re ξ = u u 3 2 u 2 4 − u 2 u 2 3 u 4 . (184) Note that the terms in u 6 must cancel by t he construction of ˜ ρ . Next, w e note that u 4 as a leading term m ust arise as either ˙ u 2 = u 2 2 ( u 4 − u ) or u 2 ˙ u = u 2 u 2 ( u 2 − u − 2 ) ≡ u 2 u 2 2 − u u 2 2 u 4 . ( 1 85) The quadratic term m ust arise from ( 185) (as w e hav e a lr eady determined all terms that in volv e u 4 , so, w e cannot hav e a quadratic term given by u 4 ˙ u 2 ) . Now , w e solv e (179) to get ρ (1) = − 1 2 u 2 u 2 2 + · · · . (186) Indeed, D t  − 1 2 u 2 u 2 2  = − u u 2 2 ˙ u − u 2 u 2 ˙ u 2 ≡ u u 3 2 u 2 4 − u 2 u 3 2 u 4 + low er shift terms = ξ + low er shift terms, (187) with ξ in (184). Th us, w e up date ρ = u 2 u 2 u 4 b y subtracting ρ (1) . This yields, ρ = u 2 u 2 u 4 + 1 2 u 2 u 2 2 . (188) W e readily v erify t ha t D t ρ ≡ 0 a nd, thu s, ρ in (18 8 ) is t he unique densit y o f (180) that dep ends on 4 shifts. 13 The Gardner Lattice In this section we consider the D DE describ ed in [91], ˙ u =  1 + αh 2 u + β h 2 u 2   1 h 3  1 2 u − 2 − u − 1 + u 1 − 1 2 u 2  + α 2 h  u − 1 u − 2 + u 2 − 1 + u ( u − 1 − u 1 ) − u 2 1 − u 1 u 2  + β 2 h  u 2 − 1 ( u − 2 + u ) − u 2 1 ( u + u 2 )   , (189) 40 whic h is a completely integrable discretization of the Gar dner equation, u t + 6 α uu x + 6 β u 2 u x + u 3 x = 0 . (190) Therefore, w e call (189 ) the Gardner lattice. Note t hat (190) is a com binatio n of the K dV equation ( β = 0) and the mKdV equation ( α = 0) . Cons equen tly , (189 ) includes the com- pletely in t egra ble discretizations of the K dV and mKdV equations as sp ecial cases. 13.1 Deriv ation of the Gardner Lattice Based on w or k b y T aha [91], w e outline the deriv ation o f the Gardner lattice f r o m a discretize d v ersion [4] of the eigen v alue problem of Zakharov and Shabat. Consider the discrete system [1, 7], V 1 ,n +1 = z V 1 ,n + Q n ( t ) V 2 ,n , V 2 ,n +1 = 1 z V 2 ,n + R n ( t ) V 1 ,n , (191) ˙ V 1 ,n = A n ( z ) V 1 ,n + B n ( z ) V 2 ,n , ˙ V 2 ,n = C n ( z ) V 1 ,n + D n ( z ) V 2 ,n , (192) where, in general, the co eﬃcien ts A n through D n dep end on the p oten tials, R n and Q n . The eigen v alue z is assumed to b e time-indep enden t. Step 1 : Expressin g the compatibility conditions, ˙ V i,n +1 = D ˙ V i,n , for i = 1 , 2 , a nd equating the co eﬃcien ts of V 1 ,n and V 2 ,n , w e get [7] C n Q n − B n +1 R n = z ( A n +1 − A n ) , B n R n − C n +1 Q n = 1 z ( D n +1 − D n ) , (193) ˙ Q n + D n Q n − A n +1 Q n = 1 z B n +1 − z B n , ˙ R n + A n R n − D n +1 R n = 1 z C n +1 − z C n . (194) Step 2 : Substituting the expansions [91 ], A n = 2 X j = − 2 z 2 j A (2 j ) n , D n = 2 X j = − 2 z 2 j D (2 j ) n , (195) B n = 2 X j = − 1 z 2 j − 1 B (2 j − 1) n , C n = 2 X j = − 1 z 2 j − 1 C (2 j − 1) n , (196) in to (193) and (194) and setting the co eﬃcien ts of lik e p ow ers of z equal to zero, w e ob- tain a system of tw en ty equations fo r t he eightee n unknow ns functions A (2 j ) n , D (2 j ) n with j = − 2 , − 1 , 0 , 1 , 2 ; and B (2 j − 1) n , C (2 j − 1) n with j = − 1 , 0 , 1 , 2 . The simplest equations arise from the co eﬃcien ts of the terms in z ± 5 and z ± 4 : A (4) n +1 − A (4) n = 0 , D ( − 4) n +1 − D ( − 4) n = 0 , (197) Q n ( D (4) n − A (4) n +1 ) + B (3) n = 0 , R n ( A ( − 4) n − D ( − 4) n +1 ) + C ( − 3) n = 0 , (198) Q n ( D ( − 4) n − A ( − 4) n +1 ) − B ( − 3) n +1 = 0 , R n ( A (4) n − D (4) n +1 ) − C (3) n +1 = 0 , (199) Step 3 : W e outline the solution pro cess whic h fo llo ws the strategy in [7 , Section 2.2a]. F rom (197), w e conclude that A (4) n and D ( − 4) n are independen t of n. Hence, A (4) n = ˜ A (4) and 41 D ( − 4) n = ˜ D ( − 4) are constan ts. The tilde will remind us that w e a re dealing with c onstants . Solving (198) and ( 199), we get B (3) n = ( ˜ A (4) − D (4) n ) Q n , C ( − 3) n = ( ˜ D ( − 4) − A ( − 4) n ) R n , (200) B ( − 3) n = ( ˜ D ( − 4) − A ( − 4) n ) Q n − 1 , C (3) n = ( ˜ A (4) − D (4) n ) R n − 1 . (201) Substituting these results into t w o of the equations coming from z ± 3 , w e ﬁnd that A ( − 4) n = ˜ A ( − 4) and D (4) n = ˜ D (4) are constan ts. F rom equations corresp onding to z ± 3 , z ± 1 , w e o btain ∆ A (2) n = A (2) n +1 − A (2) n = ˜ α ( Q n R n − 1 − Q n +1 R n ) , (202) ∆ D ( − 2) n = D ( − 2) n +1 − D ( − 2) n = ˜ β ( R n Q n − 1 − R n +1 Q n ) , (203) ∆ A ( − 2) n = A ( − 2) n +1 − A ( − 2) n = ˜ β ( Q n R n +1 − Q n − 1 R n ) , (204) ∆ D (2) n = D (2) n +1 − D (2) n = ˜ α ( R n Q n +1 − R n − 1 Q n ) , (205) where ˜ α = ˜ A (4) − ˜ D (4) and ˜ β = ˜ D ( − 4) − ˜ A ( − 4) . Equations (202)- (205) ar e satisﬁed f o r A (2) n = ˜ A (2) − ˜ αQ n R n − 1 , D ( − 2) n = ˜ D ( − 2) − ˜ β R n Q n − 1 , (206) A ( − 2) n = ˜ A ( − 2) + ˜ β Q n − 1 R n , D (2) n = ˜ D (2) + ˜ αR n − 1 Q n , (207) where ˜ A (2) , ˜ D ( − 2) and ˜ A ( − 2) , ˜ D (2) are constan ts. Next, w e solv e equations (from the terms in z ± 2 ) for B ( ± 1) n and C ( ± 1) n , yielding B (1) n = ˜ δ Q n + ˜ αQ n +1 − ˜ αQ n ( Q n +1 R n + Q n R n − 1 ) , (208) C ( − 1) n = ˜ γ R n + ˜ β R n +1 − ˜ β R n ( R n +1 Q n + R n Q n − 1 ) , (209) B ( − 1) n = ˜ γ Q n − 1 + ˜ β Q n − 2 − ˜ β Q n − 1 ( R n − 1 Q n − 2 + R n Q n − 1 ) , (210) C (1) n = ˜ δ R n − 1 + ˜ αR n − 2 − ˜ αR n − 1 ( Q n − 1 R n − 2 + Q n R n − 1 ) , (211) where ˜ δ = ˜ A (2) − ˜ D (2) and ˜ γ = ˜ D ( − 2) − ˜ A ( − 2) . Next, w e solve equations coming fr o m z ± 1 , whic h in volv e ∆ A (0) n and ∆ D (0) n . These equations, which are similar in structure to (20 2)-(205), yield A (0) n = ˜ α  − Q n +1 R n − 1 − Q n R n − 2 + Q n Q n +1 R n − 1 R n + Q n − 1 Q n R n − 2 R n − 1 + Q 2 n R 2 n − 1  + ˜ A (0) − ˜ δ R n − 1 Q n , (212) D (0) n = ˜ β  − R n +1 Q n − 1 − R n Q n − 2 + R n R n +1 Q n − 1 Q n + R n − 1 R n Q n − 2 Q n − 1 + R 2 n Q 2 n − 1  + ˜ D (0) − ˜ γ Q n − 1 R n , (213) where ˜ A (0) and ˜ D (0) are constants. All the diﬀerence equations for the unknow n co eﬃcien ts are no w satisﬁed. W e a re left with the t wo D DEs (coming from the terms in z 0 ) , ˙ Q n − ˜ κQ n + (1 − Q n R n ) h − ˜ αQ n +2 − ˜ δ Q n +1 + ˜ γ Q n − 1 + ˜ β Q n − 2 + ˜ αQ n +1 ( Q n +2 R n +1 + Q n +1 R n ) − ˜ β Q n − 1 ( Q n − 2 R n − 1 + Q n − 1 R n ) + ˜ αQ n Q n +1 R n − 1 − ˜ β Q n Q n − 1 R n +1 i = 0 , (214) ˙ R n + ˜ κR n + (1 − R n Q n ) h − ˜ β R n +2 − ˜ γ R n +1 + ˜ δ R n − 1 + ˜ αR n − 2 + ˜ β R n +1 ( R n +2 Q n +1 + R n +1 Q n ) − ˜ αR n − 1 ( R n − 2 Q n − 1 + R n − 1 Q n ) + ˜ β R n R n +1 Q n − 1 − ˜ αR n R n − 1 Q n +1 i = 0 , (215) 42 where ˜ κ = ˜ A (0) − ˜ D (0) . Step 4 : The terms − ˜ κQ n in (214) and ˜ κR n in (215) could b e remov ed with a suitable transformation. Accomplishing t he same, w e set ˜ κ = 0 . Next, we substitute Q n = u n , R n = − h 2 ( α + β u n ) , (216) in to (214) and (215). Next, w e set ˜ β = ˜ α and ˜ γ = ˜ δ to remov e all constant terms. Doing so, (214) and (215) collapse in to a single DDE, ˙ u n + (1 + α h 2 u n + β h 2 u 2 n ) n ˜ α ( u n − 2 − u n +2 ) + ˜ δ ( u n − 1 − u n +1 ) + α ˜ αh 2 ( u n − 2 u n − 1 + u 2 n − 1 + u n ( u n − 1 − u n +1 ) − u 2 n +1 − u n +1 u n +2 ) + β ˜ αh 2 ( u 2 n − 1 ( u n − 2 + u n ) − u 2 n +1 ( u n + u n +2 ))  = 0 , (217) Step 5 : T o ﬁx the scale on t, as w ell as the constan ts ˜ α, ˜ β , and ˜ δ , we consider the limit of (217) for h → 0 . Using, lim h → 0 u n ( t ) = u ( x, t ) , lim h → 0 ˙ u n ( t ) = u t ( x, t ) , and substituting lim h → 0 u n + m ( t ) = u ( x, t ) + mhu x ( x, t ) + 1 2 ( mh ) 2 u 2 x ( x, t ) + 1 6 ( mh ) 3 u 3 x ( x, t ) + . . . , (218) in to (217), w e obtain u t − 2 h ( ˜ δ + 2 ˜ α ) u x − 2 h 3 ( ˜ δ + 8 ˜ α )  αuu x + β u 2 u x + 1 6 u 3 x  + O ( h 5 ) = 0 . (219) The term in h disapp ears when ˜ δ = − 2 ˜ α. Substituting ˜ α = − 1 2 , ˜ δ = 1 , into (217) and (21 9), w e get ˙ u =  1 + αh 2 u + β h 2 u 2   1 2 u − 2 − u − 1 + u 1 − 1 2 u 2 + 1 2 αh 2  u − 1 u − 2 + u 2 − 1 + u ( u − 1 − u 1 ) − u 2 1 − u 1 u 2  + 1 2 β h 2  u 2 − 1 ( u − 2 + u ) − u 2 1 ( u + u 2 )  , (220) and u t + h 3 (6 αuu x + 6 β u 2 u x + u 3 x ) = 0 , (221) where the O ( h 5 ) terms w ere ignored. Using a scale, t → h 3 t, w e get (189) and (190). Remarks. Step 2 is w ell suited for a CAS. Solving the system, as outlined in Step 3, is a c hallenging task, in particular, if attempted with p en and pap er. The system consists of t wo DDEs and eightee n diﬀerence equations. No ne of the CAS has a build-in solv er for suc h mixed systems . A fully automated solution is therefore imp ossible. W e used a feedback mec hanism whic h mimics what one would do by hand: Solv e the simplest diﬀerence equations; en ter that partial solution; let CAS simplify the entire sys tem; rep eat the pro cess un til all diﬀerence equations are satisﬁed and the t wo D DEs are simpliﬁed as far as p ossible. Con tin uing with (214) and (215), steps 4 and 5 w ere straightforw ard to implemen t. The ﬁv e steps hav e b een implemen ted in Mathematica [50]. Starting from (19 1), the co de generates the Gardner lat t ice (220). The co de could b e mo diﬁed to assist in the deriv ation o f other completely in tegrable DDEs, suc h a s discrete v ersions of the nonlinear Schr¨ odinger and sine-Gordon equations. 43 13.2 Dilation Inv ariance of the Gardner Lattice Since (189) is not uniform in rank w e m ust introduce auxiliary parameters with w eigh t. This can b e done in sev eral wa ys. One of the p ossibilities is to replace (1 89) by ˙ u =  γ + α u + β u 2   γ  1 2 u − 2 − u − 1 + u 1 − 1 2 u 2  + α 2  u − 1 u − 2 + u 2 − 1 + u ( u − 1 − u 1 ) − u 2 1 − u 1 u 2  + β 2  u 2 − 1 ( u − 2 + u ) − u 2 1 ( u + u 2 )   , ( 2 22) where w e hav e set h = 1 (by scaling) and introduced a parameter γ . Expressing uniformity of rank, setting W ( D t ) = 1 , and solving the linear system for the w eigh ts, one ﬁnds that W ( u ) = W ( α ) = 1 4 , W ( β ) = 0 , and W ( γ ) = 1 2 . So, w e do no t need a scale on β and (22 2) is in v ariant under the scaling symmetry ( t, u, α, γ ) → ( λ − 1 t, λ 1 4 u, λ 1 4 α, λ 1 2 γ ) . (223 ) F or β = 0 , (222 ) reduces to ˙ u = ( γ + αu )  γ  1 2 u − 2 − u − 1 + u 1 − 1 2 u 2  + α 2  u − 1 u − 2 + u 2 − 1 + u ( u − 1 − u 1 ) − u 2 1 − u 1 u 2  , (224) whic h is a completely integrable discretization of the KdV equation, u t + 6 αu u x + u 3 x = 0 . Computing the weigh ts, o ne can set W ( α ) = 0 , whic h leads to W ( u ) = W ( γ ) = 1 2 . So, (22 4 ) is in v aria nt under the scaling sy mmetry ( t, u , γ ) → ( λ − 1 t, λ 1 2 u, λ 1 2 γ ) . F or α = 0 , (222) reduces to ˙ u =  γ + β u 2   γ  1 2 u − 2 − u − 1 + u 1 − 1 2 u 2  + β 2  u 2 − 1 ( u − 2 + u ) − u 2 1 ( u + u 2 )   , (2 25) whic h is a completely inte grable discretization of t he mKdV equation, u t + 6 β u 2 u x + u 3 x = 0 . In this case, one can set W ( β ) = 0 . Then W ( u ) = 1 4 and W ( γ ) = 1 2 . Th us, (225) is inv arian t under the scaling symmetry ( t, u , γ ) → ( λ − 1 t, λ 1 4 u, λ 1 2 γ ) . 13.3 Conserv ation La ws of the G ardner Lattice One can either apply the metho d of Section 12 directly to (189) or , alternative ly , apply to (222) the techniq ue based on dilation inv ariance outlined in Sections 9.2, 11.2, and 11.4. In particular, one can use t he “ divide and conquer” strategy of Section 1 1.5 t o split candidate densities in to smaller pieces. Computational details can b e f ound in [34] The results b elow w ere computed [58] with the metho d in Section 12. F or q = 0 shifts there are tw o (non- p olynomial) densit y-ﬂux pairs: ρ (0) 1 = ln (1 + αh 2 u + β h 2 u 2 ) , (226) J (0) 1 = − 1 2 n α h ( u − 2 − u − 1 − u + u 1 ) + β h (2 u − 2 u − 4 u − 1 u + 2 u − 1 u 1 ) + α 2 h  u − 1 ( u − 2 + u − 1 + u ) + u ( u − 1 + u + u 1 )  + α β h  u 2 − 1 u − 2 + 2 u − 2 u − 1 u + 3 u 2 − 1 u +3 u − 1 u 2 + 2 u − 1 uu 1 + u 2 u 1  + 2 β 2 h u − 1 u  u − 2 u − 1 + u − 1 u + uu 1  o , ( 2 27) 44 and ρ (0) 2 = arctanh h ( α + 2 β u ) p α 2 h 2 − 4 β ! , (228) J (0) 2 = 1 4 p α 2 h 2 − 4 β n 1 h 2 (2 u − 1 − u − 2 − 2 u 1 + u 2 ) + α  u − 2 u − 1 + u 2 − 1 + 2 u − 1 u + u 2 + uu 1  + β  u 2 − 1 ( u − 2 + u ) + u 2 ( u − 1 + u 1 )  o . (229) The next t w o (o f inﬁnitely man y p olynomial) densities are ρ (1) = u u 1 + α β u, (230) ρ (2) = u u 2  1 + αh 2 u 1 + β h 2 u 2 1  + α h 2 uu 1 ( u + u 1 ) + 1 2 β h 2 u 2 u 2 1 + α β (1 − α 2 β h 2 ) u + α 2 2 β h 2 u 2 , (231) where the asso ciated ﬂuxes hav e b een omitted due to length. Sp ecial Cases . W e consider t wo imp o rtan t sp ecial cases. The ﬁrst few densities for (189) with β = 0 are: ρ (0) 1 = ln ( 1 + αh 2 u ) , ρ (0) 2 = u, (232) ρ (1) = 1 2 u 2 + uu 1 , ρ (2) = uu 2 (1 + αh 2 u 1 ) + αh 2 u ( u 2 1 + uu 1 + 1 3 u 2 ) . (233) The ﬁrst few densities for (18 9 ) with α = 0 are ρ (0) 1 = ln ( 1 + β h 2 u 2 ) , ρ (0) 2 = arctan ( p β hu ) , (234) ρ (1) = u u 1 , ρ (2) = uu 2 (1 + β h 2 u 2 1 ) + 1 2 β h 2 u 2 u 2 1 . (235) 14 Additio nal Example s of Non linear DDEs 14.1 The B ogo ya vlenskii Lattice The Bogoy a vlenskii lattice [21] and [90, Eq. (17.1.2)], ˙ u = u p Y j = 1 u j − p Y j = 1 u − j ! , (236) is a generalization o f the KvM lattice (110). F or p = 2 , lattice (23 6) b ecomes ˙ u = u ( u 1 u 2 − u − 1 u − 2 ) , (237) whic h is in v ariant under the follow ing scaling symmetry ( t, u ) → ( λ − 1 t, λ 1 2 u ) . (238) 45 Lattice (237) has the follo wing densit y- ﬂux pair s (of inﬁnitely many): ρ (0) = ln u , J (0) = − ( u − 1 u − 2 + u − 1 u + uu 1 ) , (239) ρ (1) = u , J (1) = − uu − 1 ( u − 2 + u 1 ) , (240) ρ (2) = u u 1 , J (2) = − u − 1 uu 1 ( u − 2 + u + u 2 ) , (241) ρ (3) = u u 1 ( 1 2 uu 1 + u 1 u 2 + u 2 u 3 ) , (242) J (3) = − u − 1 uu 1 ( u − 2 uu 1 + u 2 u 1 + u − 2 u 1 u 2 + 2 uu 1 u 2 + u 1 u 2 2 + u − 2 u 2 u 3 + uu 2 u 3 + u 2 2 u 3 + u 2 u 3 u 4 ) . (243) F or ( 237), w e also computed the densities ρ (4) through ρ (9) . Ev ery time the rank increases b y one, the n um b er of terms in the density increases by a factor three. F or example, ρ (9) has 2187 terms and the highest shift is 15 . 14.2 The B elo v-Ch altikian L attice The Belo v-Chaltikian lattice [18, Eq. (12)], ˙ u = u ( u 1 − u − 1 ) + v − 1 − v , ˙ v = v ( u 2 − u − 1 ) , (244) in in v ariant under the scaling symmetry ( t, u, v ) → ( λ − 1 t, λu, λ 2 v ) . (245) The ﬁrst few densit y-ﬂux pairs (of inﬁnitely man y) are ρ (1) = u, J (1) = − u − 1 u + v − 1 , (246) ρ (2) = 1 2 u 2 + uu 1 − v , J (2) = − u − 1 u 2 − u − 1 uu 1 + uv − 1 + u 1 v − 1 + u − 1 v , (247) ρ (3) = u ( 1 3 u 2 + uu 1 + u 2 1 + u 1 u 2 − v − 2 − v − 1 − v − v 1 ) , (248) where J (3) has b een omitted due to length. Our results matc h these in [84]. 14.3 The B laszak-M arc iniak Lattices In [1 9], Blaszak and Marciniak used the R matrix appro a c h to deriv e families of in tegrable lattices inv olving three and f our ﬁelds. Below w e consider tw o cases inv olving three ﬁelds. Examples based on fo ur ﬁelds could b e dealt with in a similar fashion [113]. The Blaszak-Marciniak three ﬁeld lattice I [84, Eq. (2 ) ], ˙ u = w 1 − w − 1 , ˙ v = u − 1 w − 1 − u w , ˙ w = w ( v − v 1 ) , (249) is in v ariant under the scaling symmetry ( t, u, v , w ) → ( λ − 1 t, λ 1 2 u, λv , λ 3 2 w ) . (250) 46 W e computed the following densit y-ﬂux pairs o f ( 249), whic h is a completely in tegrable lattice: ρ (0) = ln w , J (0) = v , (251) ρ (1) = u, J (1) = − w − 1 − w , (252) ρ (2) = v , J (2) = u − 1 w − 1 , (2 53) ρ (3) = 1 2 v 2 + uw , J (3) = u − 1 v w − 1 − w − 1 w , (254) ρ (4) = 1 3 v 3 + uv w + uv 1 w − w w 1 , (255) J (4) = w − 1 ( u − 1 v 2 + u − 1 uw − v w − v 1 w ) . (256) Our results conﬁrm those in [84] and [112]. The Blaszak-Marciniak three ﬁeld lattice I I [112, Eq. (1.4)], ˙ u = v 1 − v + u ( w − 1 − w ) , ˙ v = v ( w − 2 − w ) , ˙ w = u 1 − u, (257) is in v ariant under the scaling symmetry ( t, u, v , w ) → ( λ − 1 t, λ 2 u, λ 3 v , λw ) . (258) The ﬁrst few densit y-ﬂux pairs for (257), whic h is completely inte grable, a re ρ (1) = w , J (1) = − u, (2 59) ρ (2) = 1 2 w 2 − u , J (2) = v − uw − 1 , (260) ρ (3) = 1 3 w 3 + v − uw − 1 − u w , J (3) = u − 1 u + v w − 2 + v w − 1 − u w 2 − 1 , (261) ρ (4) = 1 4 w 4 + 1 2 u 2 + uu 1 + v w − 2 + v w − 1 − u w 2 − 1 + v w − uw − 1 w − uw 2 , (262) J (4) = − u − 2 v − uv − u − 1 v 1 + v w 2 − 2 + 2 u − 1 uw − 1 + v w − 2 w − 1 + v w 2 − 1 − uw 3 − 1 + u − 1 uw . ( 2 63) Our results agr ee with those in [8 4 ] and [112]. 15 Soft wa re to Compute Cons erv ation La ws for PDEs and DDEs W e ﬁrst discuss our pac k ages for conserv ation la ws of PD Es and DDEs, follo we d b y a brief summary of sym b olic co des deve lop ed by other research ers. 15.1 Our M athematica and Maple Soft w are The pac k age TransPDE DensityFlux .m [9], automates the computatio n of conserv ation laws of nonlinear PDEs in (1 + 1) dimensions. In addition to p o lynomial PDEs, the softw are can han- dle PDEs with transcenden tal nonlinearities. The results in Sections 3 and 5 w ere computed with TransPDEDensity Flux.m a nd cross-c hec ke d with t he newe st vers ion of condens.m , in- tro duced in [38]. W e used TransPDEDen sityFlux.m to compute the densit y-ﬂux pairs for the examples in Sections 6 and 7. Details ab out the algorithm and a discussion of implemen tation issues can b e found in [8]. The co de DDEDensityFlux.m [3 5] was used to compute the conserv ation la ws in Sections 9 and 11. The results we re cross-c hec k ed with the latest v ersion of diff dens.m , featured in 47 [39]. Using m ultiple scales, the eﬃciency of DDEDensityFlux.m w a s dra stically improv ed. Nonetheless, the algorithms [34] within DDEDensityFlu x.m are impractical for ﬁnding den- sities and ﬂuxes of high rank. Therefore, we used t he new Maple library discrete [5 6] to compute the results in Sections 13 and 14. Some of the features of earlier v ersions of condens.m and diffdens.m w ere combined in to the InvariantsSym metries.m [39, 4 1], whic h a llo ws one to compute generalized symmetries as w ell as conserv ed densities (but no ﬂuxes). InvariantsSymmetrie s.m is a v ailable from MathSource, the Mathematica pro g ram ba nk of W olfram Researc h, Inc. Our Mathematica pac k ages and noteb o oks are av ailable at [48] and Hic kman’s Maple co de is a v ailable at [56]. Our Mathematica co des for the contin uous and discrete Euler a nd homotop y op erators in one dimension are av ailable at [49]. W e a r e curren t ly designing a comprehensiv e pack age to compute conserv ation laws o f PDEs in multiple space dimensions [45, 51, 83]. Our co des ha ve b een used in a v ariety of researc h pro jects. F or example, condens.m [38] was used by Sako vic h a nd Tsuc hida [85, 95] to compute conserv at io n laws of nonlinear Sc hr¨ odinger equations. In [84], Sahadev an and K housaly a use the algorithms of diffdens.m [42] and InvariantsSymm etries.m [39, 41] to compute conserv ed densities of t he Belov- Chaltikian and Blaszak-Marciniak lattices. Ergen¸ c and Karas¨ ozen [33] used our soft w are in the design of Pois son integrators fo r V olterra lattices. 15.2 Soft w are P ac k ages of Other Researc hers Our Mathematica co de for conserv atio n laws of PDEs has b een “translated” [108, 109, 110] in to a Maple pack age, called CONSLAW , whic h only handles PDEs in (1 + 1) dimensions. Ba sed on the concept of dila tion in v ar ia nce and the metho d of undetermined co eﬃcien ts, similar soft ware w as dev elop ed b y Deconinc k and Niv ala [27] and Y ang et a l. [10 7]. Our algor ithms [38, 42] for DDEs hav e b een adapted t o fully-discretized equations [36 , 37 ]. There are sev eral algorithms (see e.g. [101]) to sym b olically compute conse rv atio n la ws of nonlinear PDEs but few hav e b een fully implemen ted in CAS. W o lf ’s pa c k age ConLaw [101, 102, 1 06] computes ﬁrst in tegr a ls o f OD Es and conserv ation la ws of PDEs. ConLaw uses the REDUCE pack age CRACK [103, 104, 105], whic h contains to ols to solve ov erdetermined systems of PDEs. W olf ’s application pack ages heavily rely on the capabilities of CRACK , whic h to ok y ears to dev elop and p erfect. Unfortunately , no suc h pac k age is a v ailable in Mathematica . A common approac h is to use the link b etw een conserv ation la ws a nd symmetries as stated in No ether’s theorem [1 5, 66, 81]. Among the new est softw a re based on that approa c h is the Maple co de GeM [25] b y Cheviak ov [23, 24], whic h a llo ws one to compute conserv ation laws of systems of ODEs and PDEs based on the kno wledge of generalized symmetries. How ev er, the computation of such symmetries [47 ] is as diﬃcult a task as the direct computation of conserv ation laws for it requires solving systems o f ov erdetermined PDEs with, e.g., the Rif pac k age [4 7, 10 0]. Some metho ds circum v en t t he existence of a v aria tional principle (required b y No ether’s t heorem) [12, 20, 101, 106] but they still rely on softw are to solve OD Es or PDEs. The pack age DE APPLS [16, 77] also oﬀers commands for constructing conse rv atio n laws from (v ariational) symmetries by No ether’s t heorem, but the computation is not fully a uto- mated. Lik ewise, the pac k ag e Noether [43] in Maple allow s one to compute conserv ation la ws 48 from inﬁnitesimal symmetry generato r s corresp onding to (simple) Lagrangians. Based on the fo r mal symmetry approac h, Sok olov and Shabat [89], Mikhailo v et al. [74, 75], and Adler e t al. [10] classiﬁed completely in t egr a ble PDEs and DDEs in (1 + 1) dimensions. Unfortunately , the softw ar e used (see [38]) in the classiﬁcation is obsolete. F or comple teness, w e also men tion the pac k ages Jets b y Marv an [73] and Vessiot by Anderson [16, 77]. Both are general purp ose suite s of Maple pac k ages fo r computations on jet spaces. The commands within Jets and Vessiot use diﬀeren tial forms and adv a nced concepts from diﬀeren tial geometry . By av oiding diﬀeren tial forms, our co des w ere readily adaptable to nonlinear D DEs (not co vere d in Jets and Vessiot ). Finally , Deconinc k and Niv ala [26] dev elop ed Maple soft ware for the contin uous and dis- crete homotop y op erator s. Their co de is a v ailable at [28]. 16 Summary W e presen ted metho ds to sy m b olically compute conserv a tion la ws of nonlinear p olynomial and transcenden tal systems of PDEs in (1 + 1) dimensions and p o lynomial DDEs in one discrete v a riable. The ﬁrst part of this c hapter dealt with nonlinear PDEs for whic h we show ed the com- putation of densities and ﬂuxes in detail. Using the dilation in v ariance of the give n PDE, candidate p olynomial densities are constructed as linear combinations with undetermined co eﬃcien ts of scaling in v aria n t building blo c ks. F or p olynomial PDEs, the undetermined co eﬃcien ts are constants whic h m ust satisfy a linear system of algebraic equations. That system will b e parameterized b y constants app earing in the PDE, if a ny . F or tra nscende n tal PDEs, the undetermined co eﬃcien ts are functions whic h m uc h satisfy a linear system whic h is a mixture of algebraic equations and ODEs. The contin uous homotopy op erator is a p o w erful, algorithmic to ol to compute ﬂuxes ex- plicitly . Indeed, the homotopy op erator handles integration b y parts whic h allow ed us to in ve r t the total deriv ativ e op erator. The metho ds for p o lynomial PDEs a re illustrated with classical examples suc h as the K dV a nd Bo ussines q equations a nd the Drinfel’d-Sokolo v- Wilson system. The computation of conserv atio n la ws of system with transcenden tal nonlin- earities is applied to sine-Gor do n, sinh-Gordon, and Liouville equations. In the second part w e dealt with the sym b olic computation of conserv atio n la ws of non- linear DDEs. A gain, we used the scaling symmetries o f the DD E and the metho d of un- determined co eﬃcien ts to ﬁnd densities and ﬂuxes. In analog y with the con tinuous case, to compute the ﬂux one could use t he discrete ho mot op y operat or, whic h handles summation b y parts and inv erts the for ward diﬀerence op erator. How eve r, in comparison with the “splitting and shifting” tech nique, the discrete Euler a nd homotop y op erat ors ar e ineﬃcie n t to ols for t he sym b olic computat io n of conserv atio n law s of DDEs. The undetermined co eﬃcien t metho d is illustrated with classical examples suc h as the Kac-v an Mo erb eke , T o da and Ablo witz-Ladik lattices. There is a fundamen tal diﬀerence b etw een the con tinuous and discrete cases in the w ay densities (o f selected rank) are constructed. The total deriv ativ e has a w eight but the shift op erator do es not. Consequen tly , a densit y of a PDE is b ounded in order (with resp ect to x ) . Unfo rtunately , t here is no a priori b ound on the n umber of shifts in the densit y , unless 49 a leading order a nalysis is carr ied out. T o o verc ome t his diﬃcult y and other shortcomings of the undetermined co eﬃcien t metho d, w e presen ted a new metho d to compute conserv ed densities of DD Es. That metho d no longer uses dilation inv ariance and is no longer restricted to p olynomial conserv ation la ws. Instead of building a candidate densit y with undetermined co eﬃcien ts, one ﬁrst computes the leading order term in t he density and, second, generates the correction terms of lo we r order. The metho d is fast and eﬃcie n t since no unnecessary terms are computed. The new metho d w as illustrated using a mo diﬁed V olterra lattice as an example, and applied to lattices due to Bogoy a vlenskii, Belo v-Chaltikian, Bla szak-Marciniak, and Gardner. A deriv ation of the latt er lattice was also giv en. Ac kno w ledgements This material is based in part up on w ork supported by the National Researc h F oundation (NRF) of South-Africa under Gran t No. F A2007032500 003. Jan Sanders (F ree Univ ersity of Amsterdam) a nd Bernard Deconinc k (Univ ersit y of W ash- ington) are thanke d for v aluable discussions. Loren “Douglas” P o ole is thanked for pro o f reading the manus cript. Willy Hereman is gra t eful for the hospitalit y and supp ort of the Departmen t of Mathe- matics and Statistics of the Unive rsit y of Canterbury (Christc hurc h, New Zealand) and the Applied Mathematics Division of the Departmen t of Mathematical Scien ces of the Unive r sit y of Stellen b osc h (Stellen b osc h, South Africa) during his sabbatical visits in A Y 2007-2 008. References [1] M. J. Ablowitz , Nonlinear ev olution equations–con tin uous and discrete, SIAM Rev. 19 (1977) 663–684. [2] M. J. Ablow it z and P . A. Clarkson, Solitons, Nonlinear Ev olution Equations and In- v erse Scattering, London Mathematical So ciet y Lecture Note Series 149 , Cam bridge Univ ersit y Press, Cambridge, 1991. [3] M. J. Ablowitz and B. M. Herbst, On homo clinic structure and nume rically induced c haos for the nonlinear Sc hr¨ odinger equation, SIAM J. Appl. Math. 50 (199 0 ) 339–351. [4] M. J. Ablowitz and J. F. Ladik, Nonlinear diﬀeren tial-diﬀerence equations, J. Math. Ph ys. 16 (1975) 598–603 . [5] M. J. Ablowitz and J. F. L a dik, A no nlinear diﬀerenc e sche me and in v erse scattering, Stud. Appl. Math. 55 (1976 ) 2 1 3–229. [6] M. J. Ablo witz and J. F. L a dik, Nonlinear diﬀerential-diﬀerence equations and F ourier analysis, J. Math. Ph ys. 17 (1976 ) 10 11–1018. [7] M. J. Ablo witz and H. Segur, Solitons a nd The Inv erse Scattering, SIAM Studies in Applied Mathematics 4 , SIAM, Philadelphia, 1981. 50 [8] P . J. Adams, Sym b olic Computatio n of Conserv ed Densities and Fluxes fo r Systems of P artial Diﬀeren t ia l Equations with T ranscenden tal Nonlinearities, Master of Science The- sis, Dept. of Mathematical a nd Computer Sciences, Color a do School o f Mines, Go lden, Colorado, 2003. [9] P . J. Adams and W. Hereman, TransP DEDensityF lux.m : A Mathemat ica pac k a ge f o r the sym b olic computatio n of conserv ed densities and ﬂuxes for systems of part ia l diﬀer- en tial equations with transcenden tal nonlinearities (2 002). [10] V. E. Adler, A. B. Shabat, and R. I. Y amilov, Symmetry approach to the integrabilit y problem, Theor. Math. Ph ys. 125 (200 0 ) 1 6 03–1661. [11] V. E. Adler, S. I. Svinolupov, and R. I. Y amilov, Multi-comp onen t V olterra and T o da t yp e integrable equations. Ph ys. Lett. A 254 (1 9 99) 24– 36. [12] S. C. Anco and G . Bluman, Direct construction metho d for conserv ation la ws of par t ia l diﬀeren tial equations. P art I: Examples of conserv ation la w classiﬁcations, Europ. J. Appl. Math. 13 (2002) 54 5 –566. [13] S. C. Anco and G . Bluman, Direct construction metho d for conserv ation la ws of par t ia l diﬀeren tial equations. Part I I: General treat ment, Europ. J. Appl. Math. 13 (2002) 567– 585. [14] I. M. Anderson, In tro duction to the v ariational bicomplex, Contemp. Math. 132 , AMS, Pro vidence, Rho de Island (19 9 2) 51–7 3. [15] I. M. Anderson, The V ariational Bicomplex, Dept. of Mathematics, Utah State Univ er- sit y , Logan, Utah (Nov em b er, 200 4) 318 pages. Manusc ript of b o ok av ailable at http://www. math.usu.e du/ ~ fg_mp/Publi cations/VB /vb.pdf. [16] I. M. Anderson, The V essiot P a c k age, Dept. of Mathematics, Utah Sta t e Univ ersity , Logan, Utah (Nov ember, 2004) 253 pages. Av ailable at http://www.ma th.usu.edu / ~ fg_mp/Pages /Symbolics Page/VessiotDownloads.html. [17] A. Barone, F. Esp osito, C. J. Magee, and A. C. Scott, Theory and applications of the sine-Gordon equation, Riv. Nuo vo Cim. 1 (1 971) 22 7 –267. [18] A. A. Belo v and K. D. Chaltikian, Lattice analogues of W − algebras and classical inte - grable equations, Ph ys. Lett. B 309 (1993) 268–274. [19] M. Blaszak and K. Marciniak, R matrix approach to lattice integrable systems, J. Math. Ph ys. 35 (1994) 4661–46 8 2. [20] G. W. Bluman and S. C. Anco, Symmetry and In tegration Metho ds for Diﬀerential Equations, Applied Mathematical Sciences 154 , Springer V erlag, Berlin, 2002 . [21] O. I. Bo go ya vlenskii, Some constructions of in tegra ble dynamical systems , Math. USSR Izv. 31 (1988) 47–75. 51 [22] K. A. Bronnik o v, Scalar v acuum structure in general relativit y and alternate theories: conformal con tinuations, Acta Ph ys. P olonica B 11 (20 0 1) 3571 – 3592. [23] A. F. Cheviak ov, G eM softw are pac k a ge for computation of symmetries and conserv ation la ws of diﬀerential equations, Comp. Ph ys. Comm. 176 (2007) 48–61. [24] A. F. Cheviak ov, GeM sym b olic pac k age for symmetry and conserv at io n law computa- tion. Applications and recen t extens ions, Pro c. 2008 In t. Symp. on Sym b olic a nd Alge- braic Computatio n (ISSA C 2008), Hagenberg, Austria, July 20-2 3, 2008, to app ear. [25] A. F. Cheviak ov, GeM soft ware pac k age and do cumen tatio n, Av aila ble at http://www. math.ubc.ca / ~ alexch/gem/ . [26] B. D econinc k and M. Niv ala, Sym b olic in tegratio n using ho mo t op y metho ds, Preprint, Departmen t of Applied Mathematics, Univ ersity of W a shington, Seattle, W A 98195 - 2420 (2007), Math. Comput. Sim ula t. (20 0 8) in press. [27] B. Deconinc k and M. Niv ala, Maple soft w are for the sy m b olic computation of conser- v atio n la ws o f (1 + 1) − dimensional pa r tial diﬀeren tial equations. Av ailable at http: //www.amath .washingto n.edu/ ~ bernard/pap ers.html. [28] B. Deconinc k and M. Niv ala, Maple soft w are for the sym b olic integration and summa- tion of expressions that may or may not b e exact. Av ailable at http://www.ama th. washington. edu/ ~ bernard/pap ers.html. [29] R. K. Do dd and R. K. Bullough, P olynomial conserv ed densities fo r the sine-Gordon equations, Pro c. Roy . So c. Lond. A 352 (1977) 481–503 . [30] P . G. D razin and R. S. Johnson, So litons: an In tro duction, Cam bridge Univ ersit y Press, Cam bridge, 19 89. [31] V. G. Drinf el’d and V. V. Sok olov, Equations of Kortew eg-de V ries type and simple Lie algebras, So v. Math. Dokl. 23 (19 81) 457 – 462. [32] V. G. Drinfel’d and V. V. Sok olov, Lie algebras and equations of Kortew eg- de V ries t yp e, J. Sov. Math. 30 (1985) 197 5–2036. [33] T. Ergen¸ c and B. Karas¨ ozen, P oisson in tegrat o rs fo r V olterra lattice equ ations, Appl. Num. Math. 56 (2006) 879 –887. [34] H. Eklund, Sym b o lic Computation of Conserv ed Densities and Fluxes for Nonlinear Systems of Diﬀeren tia l- Diﬀerence Equations, Master of Science Thesis, D ept. of Mathe- matical a nd Computer Sciences, Color a do School of Mines, Golden, Colora do , 20 0 3. [35] H. Eklund a nd W. Hereman, DDEDensityFlux.m : A Mathematica pac k age for the sym- b olic computation of conserv ed densities and ﬂuxes for nonlinear systems of diﬀeren tial- diﬀerence equations (200 3). [36] M. Gao, Y. Kato, and M. Ito, A REDUCE pac k age for ﬁnding conserv ed densities o f systems of nonlinear diﬀerence-diﬀerence equations, Comp. Ph ys. Comm. 148 (20 0 2) 242– 255. 52 [37] M. Gao, Y. Kato, and M. Ito, A REDUCE pac k age for ﬁnding conserv ed densities o f systems of implicit diﬀerence-diﬀerence equations, Comp. Ph ys. Comm. 160 (2 004) 69 –89. [38] ¨ U. G¨ okta¸ s and W. Hereman, Sym b olic computation of conserv ed densities fo r systems of nonlinear ev olutio n equations, J. Symb. Comp. 2 4 (1997) 591–6 2 1. [39] ¨ U. G¨ okta¸ s and W. Hereman, Computatio n o f conserv ation la ws for nonlinear latt ices, Ph ysica D 132 (1998) 425–436 . [40] ¨ U. G¨ okta¸ s and W. Hereman, Algorithmic computation of higher-order symme t ries for nonlinear ev olution and lattice equations, Adv. Comput. Math. 11 (19 99) 55–80. [41] ¨ U. G¨ okta¸ s and W. Hereman, InvariantsSymme tries.m : A Mathematica in tegr a bil- it y pac k age fo r the computation of in v ariants and symmetries (1997 ) . Av ailable from MathSource (Item: 0208-9 3 2, Applications/Mathematics). [42] ¨ U. G ¨ okta¸ s, W. Hereman, and G. Erdmann, Computation of conserv ed densities for systems of nonlinear diﬀeren tial-diﬀerence equations, Ph ys. Lett. A 236 (1997 ) 3 0 –38. [43] P . D . F. Gouv eia and D. F. M. T orres, A computer algebra pac k age for determining symmetries and conserv ation la ws in the calculus of v ariations. See: http://ww2.ma t. ua.pt/delfi m/delfim/T ransparencies/Optimization2004.pdf. [44] M. H´ enon, In tegrals of the T o da lattice, Phy s. Rev. B 9 (1974) 1921–19 2 3. [45] W. Hereman, Sym b olic computation of conserv ation laws of nonlinear partial diﬀeren tial equations in m ulti- dimensions, In t. J. Quant. Chem. 106 (2006) 278–299 . [46] W. Hereman, Shallow w a t er w a ve s and solitary wa v es, In: Encyclopedia of Complexit y and Systems Science, Ed.: R. A. Mey ers, Springer V erlag, Berlin, Germany (200 8), in press. [47] W. Hereman, Sym b o lic Soft ware for L ie Symmetry Analysis. In: CRC Handb o ok of Lie Group Analysis of Diﬀeren tial Equations 3 : Ne w T rends in The oretical D eve lopmen ts and Computational Methods, Ch. 1 3, Ed.: N. H. Ibrag imo v, CR C Press, Bo ca Raton, Florida (1996) 367–41 3. [48] W. Hereman, T ransPDEDensit yFlux.m and DD ED ensit yFlux.m: Mat hematica pac k ages for t he sym b olic computation of conserv ation laws of partial diﬀeren tial equations and diﬀeren tial-diﬀerence equations (2004). Av ailable from the soft ware section at http: //www.mines .edu/fs_ho me/whereman/. [49] W. Hereman, contin uo us1DHomotop yOp erator.m and discrete1DHomotop yOp erator.m: Mathematica pac k ages for the contin uous and discrete homotop y op erators in one dimen- sion (2005). Av a ila ble from the soft w a re section a t http://www.mines. edu/fs_home/ whereman/. 53 [50] W. Hereman, GardnerLattice.m: A Math ematica pac k age f o r the deriv ation of a com- pletely integrable diﬀeren tial- diﬀerence equation with the com bined KdV-mKdV equa- tion as the con tin uous limit (2001). Av ailable from the soft w a r e se ction at http: //www.mines .edu/fs_ho me/whereman/. [51] W. Hereman, M. Colagr osso, R. Sa yers, A. Ringler, B. Deconinck , M. Niv ala, and M. S. Hic kman, Con tin uous and discrete homotopy op erators with applications in in tegrability testing. In: D iﬀ erential Equations with Sym b o lic Computation, Eds.: D. W ang and Z. Zheng, Bir kh¨ auser V erlag, Basel, Switzerland (2005) Ch. 15, 25 5–290. [52] W. Hereman, B. Deconinc k, and L. D. Poole, Con tinuous and discrete homotop y op era- tors: A theoretical approach made concrete, Math. Comp. Sim ul. 74 (20 0 7) 352– 360. [53] W. Hereman and ¨ U. G¨ okta¸ s , In tegrability tests for nonlinear ev o lut io n equations. In: Computer Algebra Systems: A Practical Guide, Ed.: M. W ester, Wiley & Sons, New Y o r k, New Y ork (1999) Ch. 12 , 21 1 –232. [54] W. Hereman, J. A. Sanders, J. Say ers, and J. P . W ang, Sym b olic computation of p olyno- mial conserv ed densities, generalized symmetries, and recursion op erators for nonlinear diﬀeren tial-diﬀerence equations, In: Group Theory a nd Numerical Analysis, CRM Pro c. & Lect. Ser. 39 , Eds.: P . Win ternitz et al. , AMS, Pro vidence, Rho de Island (200 4) 267–282. [55] M. Hic kman, Leading order in tegrability conditions for diﬀeren tial- diﬀerence equations. J. Nonl. Math. Ph ys. (2008) in press. [56] M. Hic kman, Discrete: Maple 10/11 pac k age for the sy mbolic computation of conser- v atio n la ws of diﬀeren tial-diﬀerence equations (200 7 ). Av ailable from M.Hickman@m ath. canterbury. ac.nz . [57] M. Hic kman and W. Hereman, Computat io n of densities and ﬂuxes of nonlinear diﬀeren tial-diﬀerence equations, Pro c. Roy . So c. Lond. A 459 (2003 ) 2705– 2 729. [58] M. Hic kman and W. Hereman, A direct metho d to compute conserv a tion la ws of non- linear diﬀerential-diﬀerence equations (20 08) in preparation. [59] R. Hirota, B. Grammaticos, and A. Ramani, Soliton structure of the D rinfel’d-Sok olov- Wilson equation, J. Math. Phy s. 27 (1 9 86) 149 9–1505. [60] R. Hiro ta and J. Sa tsuma, N − soliton solution of nonlinear netw o r k equations describing a V olt erra system, J. Ph ys. So c. Jpn. 40 (19 7 6) 891–900. [61] P . E. Hydon and E. L. Mansﬁeld, A v a r iational complex for diﬀerence equations, F ound. Comp. Math. 4 (2004) 187 –217. [62] M. Kac a nd P . v an Mo erb ek e, On an explicitly soluble system of nonlinear diﬀeren tial equations related to certain T o da lattices, Adv. Math. 16 (1975 ) 160- 1 69. 54 [63] A. H. Kara and F . M. Mahomed, R elat io nship b etw een symmetries a nd conserv ation la ws, Int. J. Theor. Ph ys. 39 (2000) 23–40 . [64] A. Khare, S. Habib, and A. Saxena, Exact thermo dynamics of the double sinh-Gordon theory in (1 + 1) dimensions, Phy s. R ev. Lett. 79 (1997) 3797 – 3801. [65] A. V. Kiselev, On the geometry of Liouville equation: symme tries, conserv ation la ws, and B¨ ack lund transformatio ns, Acta Appl. Math. 72 (2 002) 33 –49. [66] I. S. K r a sil’shc hik and A. M. Vinogra do v, Symmetries and Conserv ation Laws for Dif- feren tial Equations of Mathematical Ph ysics, AMS, Providenc e, Rho de Island, 1998 . [67] M. D. Krusk al, R. M. Miura, C. S. Ga rdner, and N. J. Zabusky , Kortewe g-de V ries equation and g eneralizations. V. Uniquenes s and nonexistenc e of p olynomial conserv ation la ws, J. Math. Ph ys. 11 (1970) 952–96 0 . [68] L. Lam, Ed., In tro duction t o Nonlinear Ph ysics, Springer V erlag, New Y ork, 1997. [69] G. L. Lam b, Analytical descriptions of ultra short optical pulse propagatio n in a resonan t medium, Rev. Mo d. Ph ys. 43 (1 971) 99 –124. [70] A. L. Larsen and N. Sanc hez, Sinh-G ordon, Cosh-Gordon and Liouville equations for strings and m ulti- strings in constan t curv ature spacetimes, Ph ys. Rev. D 54 (1 996) 2801–280 7. [71] E. L. Mansﬁeld and P . E. Hydon, O n a v a riational complex for diﬀerence equations. In: The Geometrical Study of Diﬀeren tial Equations, Pro c. NSF- CBMS Conf., Eds.: J. A. Leslie a nd T. P . Robart, AMS, Prov idence, Rho de Island (2 002) 121–129 . [72] E. L. Mansﬁeld a nd G. R . W. Quisp el, T o w ards a v ariational complex for the ﬁnite elemen t metho d. In: Group Theory and Numerical Ana lysis, CRM Pro c. and Lect. Series 39 , Eds.: P . Win ternitz et al. , AMS, Pro vidence, Rho de Island (2 0 04) 207–231. [73] M. Marv an, Jets : A soft ware for D iﬀeren tial Calculus on Jet Spaces and Diﬃeties, Soft ware Guide, Silesian Unive r sit y , Opav a, Czec h Republic, Jets v. 4.9 (2003). See: http://diff iety.ac.ru /curvita/marvan.html. [74] A. V. Mikhailo v, A. B. Shabat , and V. V. Sokolo v, The symmetry approach to classiﬁ- cation of in tegrable equations. In: What is In tegrability?, Ed.: V. E. Zakharo v, Springer V erlag , Berlin (1 991) 115–184 . [75] A. V. Mikhailo v, A. B. Shabat, and R. I. Y amilo v, The symmetry approach to the classiﬁcation of nonlinear equations. Complete lists of in tegra ble systems, Usp. Mat. Nauk. 24 (1987) 3–53 . Engl. T ransl.: Russ. Math. Surv eys 42 (1987) 1–63. [76] J. W. Miles, The Kortew eg-de V ries equation: a historical essa y, J. Fluid Mec h. 106 (1981) 131–147. 55 [77] C. E. Miller, V essiot: A Maple pac k age for v ariatio nal and tensor calculus in multiple co ordinate fra mes, Master of Science Thesis, Dept. of Mathematics, Utah State Uni- v ersit y , Logan, Utah, 1999. Av aila ble at http://www.mat h.usu.edu/ ~ fg_mp/Pages / ReportsAndD issertatio ns/ReportsAndDissertations.html. V essiot is av ailable a t http://www. math.usu.e du/ ~ fg_mp/Pages /Symbolics Page/Symbolics.html. [78] R. M. Miura, The Kortew eg-de V ries equation: A surv ey of results, SIAM Rev. 18 (1976) 412–459; Errata: ibid., SIAM Rev. 19 (1 977) vi. [79] R. M. Miura, C. S. Ga rdner, and M. D. Krusk al, K o rtew eg-de V ries equation and gen- eralizations. I I. Existence of conserv ation laws and constan ts of motion, J. Math. Ph ys. 9 (1 968) 1204–1 2 09. [80] A. C. New ell, Solitons in Mathematics a nd Ph ysics, SIAM, Philadelphia, Pe nnsylv ania, 1985. [81] P . J. Olver, Applications of Lie Groups to Diﬀerential Equations, 2nd. edition, Gr a duate T exts in Mathematics 107 , Springer V erlag, New Y ork, 1993. [82] Y. Ohtsuki, T. Kato, Y. F ujim ura, and S. H. Lin, Adiabatic theory of laser-induced vibrational predesorption of ph ysisorb ed molecules: Application to a CO/NaCl system, J. Chem. Ph ys. 106 (1997) 4339 –4352. [83] L. D. P o ole, Symbolic Computation of Conserv ation La ws of Nonlinear P artia l D iﬀeren- tial Equations in ( N + 1)-Dimensions Using Homotop y Op erators, Ph.D. Thesis, Dept. of Mathematical and Computer Sciences, Colorado Sc ho ol of Mines, Go lden, Colorado, 2008. [84] R. Sahadev an and S. Khousaly a, Similarity reduction, generalized symmetries and in- tegrabilit y of Belo v-Chaltikian and Blaszak-Marciniak lat tice equations, J. Math. Ph ys. 42 (2 001) 3854–0 3 870. [85] S. Y u. Sa k ov ich and T. Tsuc hida, Symmetrically coupled higher-order nonlinear Sc hr¨ odinger equations: singularity analysis and integrabilit y , J. Ph ys. A: Math. Gen. 33 (2 000) 7217–7 2 26. [86] J. Sanders and J. P . W ang, Classiﬁcation of conserv ation law s f or KdV-lik e equations, Math. Comput. Sim ula t . 44 (1997 ) 471–4 81. [87] J. M. Sanz-Serna, An explicit ﬁnite-diﬀerence sc heme with exact conserv a tion prop erties, J. Comput. Ph ys. 47 (1982 ) 199–2 10. [88] J. M. Sanz-Serna, Accuracy and conserv ation prop erties in n umerical in tegratio n: the case o f the Kortew eg-de V ries equation, Numer. Math. 75 (199 7) 4 21–445. [89] V. V. Sok olo v a nd A. B. Sha ba t, Classiﬁcation of integrable ev o lut io n equations. In: So v. Scien t. Rev. Sec. C Math. Ph ys. Rev. 4 , Harw o o d Academic Publishers, New Y ork (1984) 221–280. 56 [90] Y. B. Suris, The Problem o f Discretization: Hamiltonian Approac h, Progress in Mathe- matics 209 , Birkh¨ auser V erlag, Basel, 2003. [91] T. R. T aha, A diﬀeren tial-diﬀerence equation for a KdV- MKdV equation, Math. Comput. Sim ulat. 35 (1993) 509–51 2. [92] G. T esc hl, Jacobi Op erators and Completely In tegra ble Nonlinear Lattices, AMS Math. Surv. Monogr. 72 , AMS, Prov idence, Rho de Island, 2000. [93] M. T o da, Vibration of a chain with nonlinear interaction, J. Phy s. So c. Jpn. 22 (1967) 431–436. [94] M. T o da, Theory of nonlinear lattices, Springer V erlag, Berlin, 1 981. [95] T. Tsuc hida, U. Hideaki, and M. W a dati, Integrable semi-discretizations of the coupled nonlinear Sc hr¨ odinger equations, J. Ph ys. A: Math. Gen. 32 (1 9 99) 2239–22 6 2. [96] F. V erheest a nd W. Hereman, Conserv ation law s and solitary w av e solutions for gener- alized Sc hamel equations, Ph ys. Scripta 50 (1 994) 61 1 –614. [97] J. P . W ang, Symmetries and Conserv ation Laws of Ev olution Equations, Ph.D. Thesis, Dept. of Mathematical and Computer Sciences, F ree Univ ersit y o f Amsterdam, Amster- dam, The Netherlands, 1998 . [98] G. B. Whitham, Non- linear disp ersiv e w av es, Pro c. R oy . So c. Lond. A 283 (1965) 238– 261. [99] G. Wilson, The aﬃne Lie algebra C (1) 2 and an equation of Hirota and Satsuma, Phys . Lett. A 89 (1982 ) 332–3 34. [100] A. Wittkopf and G. Reid, The Rif Pac k age v. 1.13 (2005 ) See: Maple help fa cility and http://www. cecm.sfu.c a/ ~ wittkopf/ri f.html. [101] T. W olf , A comparison of four approac hes to the calculation of conserv ation laws, Europ. J. Appl. Math. 13 (2002 ) 129–1 52. [102] T. W olf, The REDUCE pac k age ConLaw for computing ﬁrst integrals of OD Es and con- serv ation la ws of PDEs. See: http://lie .math.broc ku.ca/crack/src/conlaw.html. [103] T. W olf, The REDUCE pac k age CRAC K fo r solving o ve rdetermined systems of ODEs and PD Es. See: http://lie. math.brock u.ca/crack/src/crack.html. [104] T. W olf and A. Brand, The computer algebra pack age CRA CK f or inv estigating PD Es. In: Pro c. ER CIM Sc ho ol on Partial Diﬀerential Equations and Group Theory , Ed.: J. F. P ommaret, G MD, Bo nn, German y (1992) 28–51 . [105] T. W olf and A. Brand, Inv estigating DEs with CRA CK and r elated programs, SIGSAM Bull. Sp ecial Issue, June 1 995, 1– 8 . 57 [106] T. W olf, A. Brand, and M. Mohammadzadeh, Computer algebra algorithms a nd rou- tines for the computation of conserv atio n laws and ﬁxing of g a uge in diﬀeren tial expres- sions, J. Sym b. Comp. 27 (199 9) 221– 2 38. [107] X.-D. Y ang, H.-Y. R uan, and S.-Y. Lou, A Maple pac k age on sym b olic computation of conserv ed densities f o r (1 + 1 ) − dimensional no nlinear ev olution systems, Comm. Theor. Ph ys. 47 (2007) 961–968 . [108] R.-X. Y ao and Z.-B. Li, An algorithm of constructing the conserv ation laws of nonlinear ev olution equations, In t. J. Mo d. Phy s. B 18 (200 4) 2 633–2639. [109] R.-X. Y ao and Z.-B. Li, CONSLA W: A Maple pac k a ge to construct the conserv ation la ws for nonlinear ev olution equations. In: Diﬀeren tial Equations with Sym b olic Com- putation, Eds.: D . W ang and Z . Zheng, Birkh¨ auser V erlag, Basel, Switzerland (200 5) 307–325. [110] R.-X. Y ao a nd Z.-B. Li, CONSLA W: A Maple pack age to construct conserv ation la ws for nonlinear ev olutio n equations, Appl. Math. Comput. 173 (2 0 06) 616 –635. [111] N. J. Zabusky , F ermi-P asta-Ulam, solitons and the fabric of nonlinear and computa- tional science: History , synergetics, and visiometrics, Chaos 15 (2005) Art. No. 015 1 02, 16 pages. [112] Z.-N. Zhu, X. W u, W. Xue, and Q. Ding , Inﬁnitely many cons erv atio n la ws for t w o Blaszak-Marciniak three-ﬁeld lattice hierarch y, Phys . Lett. A 297 (2002) 387–3 95. [113] Z.-N. Zh u, X. W u, W. Xue, and Z .- M. Zh u, Inﬁnitely man y conserv at io n law s fo r the Blaszak-Marciniak four-ﬁeld in tegrable la t tice hierar ch y, Ph ys. Lett. A 296 (200 2 ) 2 8 0– 288. [114] In Memoriam: Martin D. Krusk a l, Departmen t of Mathematics, R utgers Univ ersit y , New Br unswic k, New Jersey . See: http://www .math.rutg ers.edu/kruskal. [115] The 200 6 Steele Prizes, No t ices AMS 53 (2006 ) 464–4 7 0. 58

Direct Methods and Symbolic Software for Conservation Laws of Nonlinear Equations

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment