A comparison of four approaches to the calculation of conservation laws

Thomas Wolf Department of Mathematics, Brock University St. Catharines, ON L2S 3A1, Canada

April 11, 2001

Abstract The paper compares computational aspects of four approaches to compute conservation laws of single differential equations (DEs) or systems of them, ODEs and PDEs. The only restriction, required by two of the four corresponding computer algebra programs, is that each DE has to be solvable for a leading derivative. Extra constraints for the conservation laws can be specified. Examples include new conservation laws that are non-polynomial in the functions, that have an explicit variable dependence and families of conservation laws involving arbitrary functions. The following equations are investigated in examples: Ito, Liouville, Burgers, Kadomtsev-Petviashvili, Karney-Sen-Chu-Verheest, Boussinesq, Tzetzeica, Benney.

1 Introduction

As is well known, conservation laws play an important role in Mathematical Physics. The knowledge of conservation laws is useful in the numerical integration of partial differential equations (PDEs) [20], for example, to control numerical errors. Also, the investigation of conservation laws of the Korteweg de Vries equation was the starting point of the discovery of a number of techniques to solve evolutionary equations [25] (Miura transformation, Lax pair, inverse scattering technique, bi-Hamiltonian structures). The existence of a large number of conservation laws of a PDE (system) is a strong indication of its integrability. Conservation laws play an important role in the theory of non-classical transformations [22],[23] and in the theory of normal forms and asymptotic integrability [24]. Programs described below are able to find conservation laws involving the independent variables explicitly. Finding such conservation laws is a good challenge for the inverse scattering technique.

The purpose of the methods described below is to pose as few restrictions as possible on the differential equations (DEs) to be investigated. For example, it is not assumed that any Lie-symmetries are known, nor that the equations are equivalent to the Euler-Lagrange equations of a variational problem. Instead we attempt to solve the conservation law condition directly. The strategy will be to make a local ansatz involving only the dependent variables and their derivatives. Further, the order of the derivatives is bounded in order to obtain an over determined PDE problem which subsequently is solved with the computer algebra package Crack [33], [34].

In an earlier paper [35] three of the methods were discussed with emphasis put on the computer algebra algorithms involved. In this paper we present an additional fourth method and compare these methods in terms of complexity and functionality.

The rest of the paper is organised as follows. After specifying the notation that is used, in section 3 a reminder on issues of the equivalence of conservation laws will provide the motivation for the four approaches which are explained in section 4 followed by an overview. In section 5 related computer algebra programs are shortly described and examples are given. Extensions of the basic usage of these programs are discussed in section 6.

2 Notation and setup

We adopt the notation of the book of Olver [27] where the question of equivalence of conservation laws is described in more detail in chapter 4.3. This compact notation will be explained using the sine-Gordon equation which will be used as an example throughout the paper.

• In general derivations or generally applicable formulas the independent variables are x = (x^{1}, x^{2}, …, x^{p}). In examples independent variables are t, x and y.

• Functions are denoted in general by u = (u^{1}, u^{2}, …, u^{q}). In examples the function is u and in a few examples an additional function is v.

• Partial derivatives are written as lower index like in u_{tx}−sin u = 0. If partial derivatives are repeated n times (n > 2) then this may be indicated by writing n in front of the variable in the index like in ∂^{7}u/(∂t^{4}∂x^{3}) = u_{4t3x} (which is essentially used in the appendix due to the high derivatives occuring there).

• An upper index in brackets like u^{(n)} denotes the set of all derivatives of all components u^{1}, …, u^{q} of order up to n including order zero, i.e. the functions themselves. In our example with u = u(t, x) this would, for example, mean u^{(2)} = {u, u_{t}, u_{x}, u_{tt}, u_{tx}, u_{xx}}.

• The differential equations that are to be investigated concerning conservation laws are denoted 0 = Δ(x, u^{(n)}), in our example Δ(x, u^{(2)}) = u_{tx} − sin u. If we have q functions u^{1}, …, u^{q} then it is assumed that a system of q equations 0 = Δ_{1}, 0 = Δ_{2}, …, 0 = Δ_{q} is given.

• Latin indices i, j have a range 1, …, p used for x^{i} and Greek indices α, μ have a range 1, …, q used for u^{α} and Δ_{μ}.

• Whereas u^{(n)} denotes the set of all possible indices up to order n we will also need a way to specify any particular partial derivative of an unspecified order. For that we will use the so-called multiple index _{J}. Then any u_{x}, u_{t}, u_{2t5x} are all examples for u_{J}. An example for the use of _{J} is the notation of total derivatives:

• For the total derivative we use the symbol D:

D_{x^{i}} = ∂_{x^{i}} + ∑_{α=1}^{q} (u_{x^{i}}^{α}∂_{u^{α}} + ∑_{j=1}^{p} u_{x^{j}x^{i}}^{α}∂_{u_{x^{j}}^{α}} + …) .

Using the multi index _{J} a more compact notation is:

D_{x^{i}} = ∂_{x^{i}} + ∑_{α=1}^{q} ∑_{J} u_{Jx^{i}}^{α}∂_{u_{J}^{α}}.

• If a relation 0 = W (x, u^{(n)}) is said to be fulfilled identically for any solutions of the equation(s) 0 = Δ then this means that W can be written as a linear combination of Δ and any total derivatives of Δ with arbitrary x− and u^{(n)}−dependent coefficients which are non-singular for 0 = Δ. The adopted notation for that relation is

0 = W |_{Δ=0}

which may be described as “W vanishes modulo Δ = 0”. The way to check this circumstance in practise is

– to solve each Δ_{μ} for one leading derivative (such that no other derivative in Δ_{μ} is a derivative of its leading derivative and that no two leading derivatives from two different Δ_{μ}, Δ_{ν} coincide and that none is a derivative of another one),

– to substitute in W (x, u^{(n)}) all these eliminated leading derivatives as well as all derivatives of the leading derivatives.

– We have 0 = W |_{Δ=0} if and only if W vanishes identically after all possible substitutions have been performed.

Example: For the sine-Gordon equation 0 = Δ = u_{tx} − sin u we will take u_{tx} as leading derivative and do substitutions u_{tx} = sin u. Let us assume an expression

W = −8u_{xx}u_{xxt} + 4u_{x}^{3}u_{tx} + 8u_{x}u_{xx} cos u − 4u_{x}^{3} sin u  (1)

is given and we want to compute W |_{Δ=0}. The highest derivative of u_{tx} is u_{xxt} = ∂_{x}u_{tx}. We therefore substitute u_{txx} = ∂_{x}(sin u) = u_{x} cos u giving

W = 4u_{x}^{3}u_{tx} − 4u_{x}^{3} sin u.  (2)

Substituting now u_{tx} = sin u gives W = 0. We therefore found that W as given in (1) satisfies W |_{Δ=0} = 0.

• A conservation law of Δ will be given in form of a so-called conserved current P^{i} where each P^{i} is a differential expression in x and u, i.e. P^{i} = P^{i}(x, u^{(m)}) which has a vanishing divergence due to 0 = Δ. In other words

0 = ∑_{i} D_{x^{i}}P^{i}|_{Δ=0} = Div P |_{Δ=0} .

Example: For the current

P^{t} = −4u_{xx}^{2} + u_{x}^{4}, P^{x} = 4u_{x}^{2} cos u  (3)

we find

DivP = D_{t}P^{t} + D_{x}P^{x} = −8u_{xx}u_{txx} + 4u_{x}^{3}u_{tx} + 8u_{x}u_{xx} cos u − 4u_{x}^{3} sin u  (4)

which is equal W in (1). Because of Div P |_{Δ=0} = W |_{Δ=0} = 0 the vector P in (3) is a conserved current of the sine-Gordon equation and represents a conservation law.

• About the conserved quantity: By using Gauss law and integrating 0 = Div P over a region of the p-dimensional x^{i}-space we obtain the vanishing of a surface integral over that region:

0 = ∮ P^{i}dS_{i}.

B choosing a region with a pancake like shape with the bottom and top side lying in the p − 1 dimensional surfaces x^{1} = a = constant and x^{1} = b = constant the surface integral takes the form

0 = ∫_{x^{1}=b} P^{1}dx^{2} … dx^{p} − ∫_{x^{1}=a} P^{1}dx^{2} … dx^{p}

if we assume that P^{1} falls off sufficiently quickly for any x^{i} → ∞ and therefore ∫_{mantle} P^{i}dS_{i} → 0. The minus sign comes in because the normal vectors of the top and bottom surface point into opposite directions.

If one of the coordinates plays the role of time, say x^{1}, then ∫ P^{1}dx^{2} … dx^{p} is called a constant of motion as it does not change from any time a to any time b.

Example: Continuing our example we have as a constant of motion

∫_{−∞}^{∞}(−4u_{xx}^{2} + u_{x}^{4}) dx

if u is a solution of the sine-Gordon equation.

3 The equivalence of conservation laws

Although two conservation laws may look rather different, i.e. their two conserved currents, say P and P̃, may be different, nevertheless their information content may be the same. To have a method of counting conservation laws and of comparing them, we need a unique way to characterise them. To do that we first look at ways conservation laws can look different but be equivalent. Two conservation laws Div P = 0 and Div P̃ = 0 are equivalent if 0 = Div(P − P̃) = Div R is a trivial conservation law.

i) The first kind of equivalence of two conservation laws is the case that R = 0 for all solutions of Δ = 0, i.e. if P and P̃ differ only by multiples of Δ and by total derivatives of Δ (i.e. by D_{J}Δ). To test whether this is the case for two given conserved currents P and P̃ one has to check whether

0 = (P − P̃)|_{Δ=0}

holds. As described above one solves the equation 0 = Δ (or system of equations 0 = Δ_{μ}) for the leading derivative(s) and substitutes that (them) in P − P̃. If P − P̃ becomes identically zero then the two conservation laws based on P and P̃ are equivalent.

If the conservation laws are not yet calculated and one wants to ensure that the computation of P gives a unique result, without arbitrariness due to terms vanishing because of Δ = 0 then there is no need to solve Δ = 0 for some leading derivative(s) u_{J}. In that case one just drops from the beginning of the calculation the dependency of P on the leading derivative(s) u_{J} and all derivatives of u_{J}.

Example: When computing conserved currents P for the sine-Gordon equation 0 = u_{tx} − sin u then in the ansatz for P the components P^{t}, P^{x} are assumed to be independent of u_{tx} and derivatives of u_{tx}, like u_{ttx}, u_{txx}, … (see appendix).

ii) The second kind of equivalence of two conservation laws occurs if R^{i} = P^{i} − P̃^{i} = ∑_{j} D_{j}V^{ij} for some expressions V^{ij}(x, u^{(m)}) = −V^{ji}, anti-symmetric in i, j, because in that case Div R = ∑_{i,j} D_{i}D_{j}V^{ij} = 0 (due to the symmetry D_{i}D_{j} = D_{j}D_{i} and the anti-symmetry V^{ij} = −V^{ji}) for any functions u, not necessarily solutions u(x) of Δ = 0. The existence of V^{ij}(x, u^{(m)}) satisfying P^{i} = P̃^{i}+∑_{j} D_{j}V^{ij} may not be obvious and may require a computation checking Div P = Div P̃. For ODEs this problems does not occur as there is only one independent variable and no antisymmetric V^{ij}.

The solution to this problem is not to compare conservation laws by comparing their conserved currents P and P̃ but by comparing them by their integrating factors, for PDEs they are called characteristic functions, in the following way. For a conservation law to satisfy Div P |_{Δ=0} = 0 means that Div P is identical to a linear combination of Δ_{μ} and total derivatives D_{J}Δ_{μ}. Partial integration can rewrite that as a divergence plus a linear combination of the Δ_{μ} alone:

Div P |_{Δ=0} = 0  (5)
⟺ ∃Q_{ν}^{J} : Div P = ∑_{ν,J} Q_{ν}^{J}D_{J}Δ_{ν}  (identically in all x, u_{J}^{α})  (6)
= ∑_{ν,J} D_{J}(Q_{ν}^{J}Δ_{ν}) − D_{J}(Q_{ν}^{J})Δ_{ν} (repeated partial integration)
= Div R + ∑_{ν} Q_{ν}Δ_{ν}
⟺ Div P = ∑_{ν} Q_{ν}Δ_{ν}  (after renaming (P −R) → P) (7)

The integrating factors Q_{ν} are called characteristic functions as it is known ([27], p. 272) that for a totally non degenerate system Δ = 0, the equivalence class of conservation laws Div P |_{Δ=0} = 0 is characterised uniquely by the functions Q_{ν} up to equivalence of type i).

One can look at equation (7) as a determining equation for P and Q_{ν} as it is done by the method described in section 4.2 below. Alternatively one can formulate a system of conditions that is equivalent to (7) but which involves only functions Q_{ν}. That is achieved using the property of Euler operators (also called variational derivatives) E_{ν} = ∑_{J}(−D)_{J}∂/∂u_{J}^{ν} when acting on an expression they give identically zero iff this expression is a divergence. Conditions for the Q_{ν} are therefore

∀ν : 0 = E_{ν} (∑_{μ} Q_{μ}Δ_{μ}) = ∑_{J}(−D)_{J} (\frac{∂}{∂u_{J}^{ν}} ∑_{μ} Q_{μ}Δ_{μ}) .  (8)

Example: Allowing Q to be of 2^{nd} order, i.e. Q = Q(t, x, u, u_{t}, u_{x}, u_{tt}, u_{xx}) where we dropped the dependence on u_{tx} which is equal sin u (see remark above) we find that Q · Δ = Q · (u_{tx} − sin u) does depend on u, u_{t}, u_{x}, u_{tt}, u_{tx}, u_{xx} and the Euler operator therefore reads

E = ∂_{u} − D_{t}∂_{u_{t}} − D_{x}∂_{u_{x}} + (−D_{t})(−D_{t})∂_{u_{tt}} + (−D_{t})(−D_{x})∂_{u_{tx}} + (−D_{x})(−D_{x})∂_{u_{xx}}
= ∂_{u} − D_{t}∂_{u_{t}} − D_{x}∂_{u_{x}} + D_{t}^{2}∂_{u_{tt}} + D_{t}D_{x}∂_{u_{tx}} + D_{x}^{2}∂_{u_{xx}}.

Requiring condition(s) (8) to be satisfied identically in all x^{i}, u^{α} and derivatives of u^{α} (i.e. u_{J}^{α}) is equivalent to the condition (7). System (8) is often very large. It can be considerably shortened if it is projected onto the space of solutions |_{Δ=0} (as described above):

0 = ∑_{μ,J}(−D)_{J} (Q_{μ} \frac{∂Δ_{μ}}{∂u_{J}^{ν}} )|_{Δ=0} ∀ν.  (9)

Conditions (9) are known as adjoint symmetry conditions which are necessary but not sufficient for the Q_{μ} to be characteristic functions of first integrals.

iii) For any two conservation laws 0 = Div P and 0 = Div P̃, 0 = Div(P + P̃) is also a conservation law. By determining conservation laws with characteristic functions of successively increasing order, constant multiples of characteristic functions of lower order can be dropped.

iv) In the case of (systems of) ODEs the characteristic functions are called integrating factors, and P is a scalar, called a first integral. Any arbitrary function of first integrals is a first integral as well.

The four approaches described in the following four sections are to solve conditions (5), (7), (8) and (9).

4 The four approaches

4.1 A first approach

The first approach is to solve

Div P |_{Δ=0} = 0  (10)

directly.

The condition (10) is made over determined by restricting the P^{i} to be differential expressions in the u of at most some order k, i.e. P^{i} = P^{i}(x, u^{(k)}). Characteristic features of this approach are

(+) A single, first order PDE involving only few terms is to be solved.

(0) Characteristic functions have to be computed from P in a straight forward calculation (described in [35]). This is done within the computer algebra program ConLaw1 which implements the first approach including the computation of related characteristic functions Q_{μ}.

(−) It would be computationally expensive for a corresponding computer program to drop during the process of solving (10) any free functions V^{ij} = −V^{ji} (see the discussion of the second kind of equivalence of conservation laws in the previous section) which correspond to trivial conservation laws^{1}. Hence, the condition (10) has to be solved first in full generality and trivial conservation laws (i.e. V^{ij}) have to be identified and dropped afterwards. That means that the task for the computer program is unnecessarily hard through the presence of the trivial conservation laws in the general solution of (10). A rule of thumb says that the difficulty in solving a linear over determined PDE system depends less on the order or size of the PDE but more on the complexity of the result^{2}. That means the trivial conservation laws will complicate the solution of (10), the more so the more independent variables are present.

1An algorithm for that is given in [35].

2For example, if an over determined PDE (system) has no solution then a differential Gröber Basis calculation will quickly produce PDEs of lower and lower order until a contradiction is reached. On the other hand, if a PDE system has arbitrary functions in its general solution (as is the case with the PDE (10)) then computing a differential Gröbner Basis will not produce a system that is solvable by only integrating ODEs, it will involve PDEs

(−) In most cases the expressions for the P^{i} are more complicated than the expressions for the characteristic functions Q_{μ} which by the above rule of thumb indicates a more difficult computation than the solution of equations involving only Q_{μ}.

To illustrate and compare all four approaches we will apply each to finding conservation laws of the sine-Gordon equation

u_{tx} − sin(u) = 0.  (11)

If the program ConLaw1 is called to find conservation laws with conserved current P^{t}, P^{x} of order 0, then it will reply that it is not applicable. This is because Div P would be of first order in u, so equation (11) could not be used to substitute u_{tx} (when computing |_{Δ=0} in (10)) and therefore any conservation laws found would be valid for any function u(t, x), not necessarily only for solutions of the sine-Gordon equation. These conservation laws would therefore be trivial, falling into category ii) in section 3.

Details of higher order investigations are given in table 1 below. u^{(n)} stands for all derivatives of u of order 0 to n. u_{tx}^{(n)} stands for all derivatives of u_{tx} up to order n, for example, u_{tx}^{(1)} would be the derivatives u_{tx}, u_{ttx}, u_{txx}. Finally, u^{(n)}/u_{tx}^{(k)} stands for all derivatives of u up to order n apart from u_{tx} and all its derivatives up to order k. The conservation laws are given in the appendix (in the table only the equation number is cited). For each conservation law in the appendix (apart from the first) there exists another one resulting from the exchange t ↔ x.

The times given in the table are measured on a 266 MHz Pentium PC running a 80 MByte Reduce 3.6 session under Linux using the Sep. 1998 version of the program Crack for solving the over determined conditions. The 80 MByte were not necessary. For example, it is possible (using ConLaw2 which implements the 4^{th} method described below) to investigate up to 4^{th} order laws with 4 MByte and up to 7^{th} order laws with 8 MByte. To get this high in order with relatively low memory consumption, one has to give in Crack the study of integrability conditions a higher priority than the integration of equations. The price is a higher computing time. The times in the last column are to be understood only as very rough indicators^{3}. They depend sensitively on the order of priorities with which modules are to be used within the program Crack (see the manual [33] and about its availability the end of the section 7).

When condition (10) is solved, the P^{i} that are computed initially do not contain u_{tx} nor its derivatives. Afterwards a computation as outlined in (5) - (7) is performed such that finally ConLaw1 is able to return the conservation law in the form (7). In the process of computing this form (7) the new P^{i} may now involve u_{tx} (through R in (6) - (7)).

3For example, the computing times reported in [35] are now (July 1999) reduced by a factor of more than ten for higher orders.

order no of  independent  functions to  cons.  time to
of P^{i} terms variables, [no of var.]  compute, [no of arg.]  laws solve (10)

found

1  8  t, x, u^{(2)}/u_{tx}, [7]  P^{t}, P^{x}(t, x, u^{(1)}), [5]  (29),(30)  9 sec

2  12  t, x, u^{(3)}/u_{tx}^{(1)}, [9]  P^{t}, P^{x}(t, x, u^{(2)}/u_{tx}), [7]  (31)  38 sec

3  18  t, x, u^{(4)}/u_{tx}^{(2)}, [11] P^{t}, P^{x}(t, x, u^{(3)}/u_{tx}^{(1)}), [9] none^{4}  -

4  26  t, x, u^{(5)}/u_{tx}^{(3)}, [13] P^{t}, P^{x}(t, x, u^{(4)}/u_{tx}^{(2)}), [11] none^{5}  -

Table 1: The program ConLaw1 applied to compute conservation laws of the sine-Gordon equation.

4.2 A second approach

The next approach consists in solving

Div P = ∑_{ν} Q_{ν}Δ_{ν}  (12)

directly, i.e. finding P^{i}, Q_{μ} that satisfy (12) identically in x^{i}, u_{J}^{α}. Equations Δ = D_{J}Δ = 0 are not used for substitutions in (12) but they are used to reduce dependencies of the Q_{μ}.

The problem becomes over determined by restricting the order of the Q_{μ}, i.e. Q_{μ} = Q_{μ}(x, u^{(k)}) for some k and by taking Q_{μ} ← Q_{μ}|_{Δ=0}, i.e. having Q_{μ} independent of one leading u-derivative (and their derivatives) from each one of the equations Δ_{ν}. If Q_{μ} would be allowed to depend on all u^{(n)} which occur in (12) then this equation could simply be solved algebraically, by eliminating one of the Q_{μ}. But that would mean division through one Δ_{μ} and therefore Q_{μ} being singular for solutions of Δ_{μ} = 0.

The second approach has the following characteristics:

(+) The conservation law condition (12) is a single first order PDE as in the first approach.

(+) By calculating characteristic functions Q_{μ} and furthermore characteristic functions Q_{μ}|_{Δ=0}, conservation laws are uniquely characterized.

(+) The effort in formulating conditions is as low as in the first approach.

(0) The P^{i} and Q_{μ} are computed in one computation.

4Crack was not able to solve all equations completely because the general solution of (10) involves free functions (related to trivial conservation laws) which complicates the problem considerably for the computer program.

5The computer memory was not sufficient to complete the computation.

(0) The number of functions to compute is higher than in the first approach and also the number of derivatives of u on which these functions depend on because no substitutions are done in (12). The resulting complication is not too big as more variables means a higher over determination and simplification.

(−) If the order of Δ is n and the order of Q_{μ} is chosen to be k then the order of P^{i} at the start of the computation can be assumed without loss of generality to be max(k, n) (see [26]). ‘Without loss of generality’ means that a trivially conserved current P̃ can be subtracted from P such that P − P̃ is of order max(k, n). If the right hand side of equation (12) is known to be linear in the highest derivatives of order max(k, n) then P^{i} at the start of the computation can even be assumed without loss of generality to be of order max(k, n) − 1.

In this approach the investigations with k < n are not much simpler than the case k = n. This matters when the order n of Δ and the number p of variables x are high. Therefore this approach is not very efficient for low order conservation laws of high order equations.

For example, for zeroth order conservation laws (k = 0) of the Kadomtsev-Petviashvili equation (19) the P^{i} are taken initially as functions of the 23 variables t, x, y, u, u_{t}, u_{x}, u_{y}, u_{tt}, …, u_{yy}, u_{ttt}, …, u_{yyy} and the conservation law condition (12) is a condition in 38 variables (including the 4^{th} order u-derivatives). That is a much harder problem than the corresponding conditions (8),(9). For example, in this case condition (8) is a single 4^{th} order PDE in also 38 variables but for only one function Q of only four variables!

(−) When looking for conservation laws with the first method, gradually increasing the order of the conserved current P gives each conservation law in its lowest order form, i.e. a form where P is of minimal order. This is not necessarily the case using the 2^{nd} method. The transformation (6)-(7) adding R to P may increase the order of P. This implies an increase of complexity having to go up in order to get the equivalent conservation law. To give an example, the Tzetzeica equation u_{xt} = e^{u} −e^{−2u} (analysed in [30],[21]) has the conservation law

0 = D_{t} [3u_{xxx}^{2} − 5u_{xx}^{3} + 15u_{xx}^{3}u_{x}^{2} + u_{x}^{6}] +
D_{x} [−3e^{u} (u_{xx}^{2} + 2u_{xx}u_{x}^{2} + 2u_{x}^{4}) − 3e^{−2u} (2u_{xx}^{2} − 8u_{xx}u_{x}^{2} + u_{x}^{4})]

with a third order conserved current. (In [21] an infinite list of conservation laws is given.) Bringing the above conservation law to the form (12) as it would be found with the second method, it becomes

6 (u_{xxxxx} + 5u_{xxx}u_{xx} − 5u_{xxx}u_{x}^{2} − 5u_{xx}^{2}u_{x} + u_{x}^{5}) (u_{tx} − e^{u} + e^{−2u})
= D_{t} [3u_{xxx}^{2} − 5u_{xx}^{3} + 15u_{xx}^{2}u_{x}^{2} + u_{x}^{6}] +
D_{x} 3 [2u_{tx}u_{xxxx} − 2u_{txx}u_{xxx} + 5u_{tx}u_{xx}^{2} − 10u_{tx}u_{xx}u_{x}^{2}
+e^{u} (2u_{xxx}u_{x} − 2u_{xxxx} − 6u_{xx}^{2} + 8u_{xx}u_{x}^{2} − 2u_{x}^{4})
+e^{−2u} (2u_{xxxx} + 4u_{xxx}u_{x} + 3u_{xx}^{2} − 2u_{xx}u_{x}^{2} − u_{x}^{4})]

with a 4^{th} order conserved current.

Applying the program ConLaw3 that corresponds to the above method to the sine-Gordon equation (11) gives the following table.

order no of  independent  functions to  cons.  time to
of Q terms variables, [no of var.] compute, [no of arg.]  laws solve (12)

found

0  10  t, x, u^{(2)}, [8]  P^{t}, P^{x}(t, x, u^{(1)}), [5]  none  3 sec
Q(t, x, u), [3]
1  10  t, x, u^{(2)}, [8]  P^{t}, P^{x}, Q(t, x, u^{(1)}), [5] (29),(30) 8.3 sec

2  10  t, x, u^{(2)}, [8]  P^{t}, P^{x}(t, x, u^{(2)}), [8]  none  3.5 sec
Q(t, x, u^{(2)}/u_{tx}), [7]
3  16  t, x, u^{(3)}, [12]  P^{t}, P^{x}(t, x, u^{(3)}), [12]  none^{6}  -
Q(t, x, u^{(3)}/u_{tx}^{(1)}), [9]

Table 2: The program ConLaw3 applied to compute conservation laws of the sine-Gordon equation.

4.3 A third approach

Instead of calculating the conserved current P^{i} directly, the third approach is to calculate characteristic functions Q_{μ} first (see, e.g. Proposition 5.33 in [27]) and from them P^{i} afterwards using formulas of Anco & Bluman [2],[3],[4] in a form described in [35] or using repeatedly the Crack routine for integrating exact DEs (see also the section on Homotopy Operators in [27]). The condition (8) (as derived in [27],[2]) is:

0 = ∑_{J}(−D)_{J} (\frac{∂}{∂u_{J}^{ν}} ∑_{μ} Q_{μ}Δ_{μ}) ∀ν.  (13)

Typical features are:

(+) Equations (13) are equivalent to (12) and therefore necessary and sufficient.

(+) The usually more complicated P^{i} are eliminated and as in the 2^{nd} method, no trivial conservation laws are calculated which otherwise unnecessarily complicate the calculation.

(+) The highest u-derivatives in conditions (13) are of the order 2n where n is the order of the u-derivatives in ∑_{μ} Q_{μ}Δ_{μ}. The harder the problem, i.e. the higher n and the higher the number of variables, the more u-derivatives occur only explicitly in (13) and can be used for a direct separation (splitting). Higher over determination simplifies the solution of (13).

(−) Equations (13) consist of as many equations as there are dependent variables u^{μ}. The unknown functions Q_{μ} appear with n^{th} order derivatives.

6The computer memory was not sufficient to complete the computation.

(−) For an increasing order of the Q_{μ}, increasing number of u^{ν} and increasing number of x^{i}, the size of (13) can soon become unmanageable.

Applying the program ConLaw4 that corresponds to the above method to the sine-Gordon equation (11) gives the following table 3. The striking feature of this approach is the quick increase of the size of conditions. Apart from the order 0 case they increase by a factor of about 7 which itself is increasing slightly with the order. The size of conditions prevents going higher in the order. On the other hand, the completeness of the generated conditions simplifies the solution in difficult cases and speeds up the solution of the over determined system as long as it is not already too large at the beginning.

order no of  independent  functions to  cons.  time to
terms variables, [no of var.] compute, [no of arg.] laws solve (13)
found h:min:sec
0  7  t, x, u^{(1)}, u_{tx}, [6]  Q(t, x, u), [3]  none  0.7 sec

1  22  t, x, u^{(2)}, [8]  Q(t, x, u^{(1)}), [5]  (29)  2.8 sec

2  154  t, x, u^{(3)}, [17]  Q(t, x, u^{(2)}/u_{tx}), [7] none  4.7 sec

3  1116  t, x, u^{(4)}, [24]  Q(t, x, u^{(3)}/u_{tx}^{(1)}), [9] (31) 5 min 17 sec

4  8402  t, x, u^{(5)}, [34]  Q(t, x, u^{(4)}/u_{tx}^{(2)}), [11] none 10 h 49 min^{7}

5 64064  t, x, u^{(6)}, [41]  Q(t, x, u^{(5)}/u_{tx}^{(3)}), [13]  -  > 2 days

Table 3: The program ConLaw4 applied to compute conservation laws of the sine-Gordon equation.

4.4 A fourth approach

Projecting conditions (13) into the space of solutions of Δ = 0 we obtain

0 = ∑_{μ,J}(−D)_{J} (Q_{μ} \frac{∂Δ_{μ}}{∂u_{J}^{ν}} )|_{Δ=0} ∀ν.  (14)

Characteristic features of this method are similar to those of the third method with the following modifications:

7This time was nearly completely spent to formulate the condition and to separate it into 823 individual equations for Q. Then already the 3^{rd} step gave that Q can not depend on 4^{th} order derivatives.

(+) The conditions usually involve fewer terms than in the third approach which can be decisive. But as the conditions (14) are not sufficient, they are less over determined and may be harder to solve than those in the third approach.

(−) Because the conditions (14) investigated by this fourth method are not equivalent to (13) and sometimes less restrictive than conditions (13), the solutions of (14) need not represent conservation laws. Therefore after computing the Q_{μ} from (14), it has to be checked whether P^{i} exist that satisfy Div P = ∑_{ν} Q_{ν}Δ_{ν} ([2],[3],[4],[35]). If they do not exist then the Q_{μ} correspond to an adjoined symmetry but not to a conservation law.

(−) If the 4^{th} method finds adjoined symmetries which are not conservation laws, then it is still possible that these adjoined symmetries can be combined linearly to give conservation laws. But how to combine adjoined symmetries to give conservation laws is not answered by solving (14), it has to be investigated separately. This shall be illustrated with the following simple example.

Applying the third method through the program ConLaw4 to the ODE u′′ + u = 0 and restricting the search to integrating factors Q = a(x)u + b(x) that are linear in u, the program finds the following 5 integrating factors

cos(x)^{2}u′ + cos(x) sin(x)u, 2 cos(x)^{2}u − 2 cos(x)u′ sin(x) − u, cos(x), sin(x), u′,

and therefore it finds 5 first integrals. In comparison, using the same ansatz for Q but now using the fourth method with the less restrictive condition (14) the corresponding program ConLaw2 finds 8 solutions for Q:

−2 cos(x)^{2}u′ − 2 cos(x) sin(x)u + u′, cos(x), sin(x), u′,
cos(x)u′u + sin(x)u^{2}, cos(x)u^{2} − u′ sin(x)u, − cos(x)^{2}u + cos(x)u′ sin(x), u.

Only for the first 4 of these 8 Q-values does a P exist such that D_{x}P = Q · (u′′ + u), i.e. only 4 first integrals are found, the remaining 4 solutions represent only adjoined symmetries. This inability to find 5 first integrals is not a weakness of the computer program but the price to pay for the fourth method to have shorter conditions (14) compared with the conditions (13) of the third method. If the conditions (14) of the fourth method are less restrictive than the conditions (13) of the third method, why then does it find fewer first integrals, 4 instead of 5? The answer is that the 5^{th} conservation law is contained in the second half of the 8 solutions: the 8^{th} solution for Q plus 2 times the 7^{th} solution gives a value for Q that is not only an adjoined symmetry but also an integrating factor for an additional first integral.

To summarise, the fourth method gives shorter, more manageable conditions which, strictly speaking, do not compute conservation laws but adjoined symmetries. Conservation laws can be derived through appropriate linear combinations of adjoined symmetries which has to be investigated separately. This theoretical weakness of the fourth method does usually play no role in practical applications.

Applying the program ConLaw2 that corresponds to the above method to the sine-Gordon equation (11) gives the following table. The typical feature of this approach is the slower increase of the size of conditions. Apart from the order 0 case they increase by a factor of about 2 which itself is increasing slightly with the order. Compared with the previous method the size of conditions grows slower which allows to go higher in the order. Because the conditions that are generated are only necessary, not sufficient, they are slightly more difficult and expensive to solve. This causes longer running times for low order investigations. Time limitations could be overcome to some extend by faster computers.

order no of  independent  functions to  cons.  time to
of Q terms variables, [no of var.] compute, [no of arg.] laws  solve (14)
found h:min:sec
0  6  t, x, u^{(1)}, [5]  Q(t, x, u), [3]  none  1 sec

1  21  t, x, u^{(2)}/u_{tx}, [7]  Q(t, x, u^{(1)}), [5]  (29)  4.3 sec

2  45  t, x, u^{(3)}/u_{tx}^{(1)}, [9]  Q(t, x, u^{(2)}/u_{tx}), [7] none  12 sec

3  99  t, x, u^{(4)}/u_{tx}^{(2)}, [11] Q(t, x, u^{(3)}/u_{tx}^{(1)}), [9] (31)  50 sec

4  202  t, x, u^{(5)}/u_{tx}^{(3)}, [13] Q(t, x, u^{(4)}/u_{tx}^{(2)}), [11] none 2 min 43 sec

5  435  t, x, u^{(6)}/u_{tx}^{(4)}, [15] Q(t, x, u^{(5)}/u_{tx}^{(3)}), [13] (32) 16 min 10 sec

6  870  t, x, u^{(7)}/u_{tx}^{(5)}, [17] Q(t, x, u^{(6)}/u_{tx}^{(4)}), [15] none 49 min 20 sec

7  1836 t, x, u^{(8)}/u_{tx}^{(6)}, [19] Q(t, x, u^{(7)}/u_{tx}^{(5)}), [17] (33)  8 h

8  3643 t, x, u^{(9)}/u_{tx}^{(7)}, [21] Q(t, x, u^{(8)}/u_{tx}^{(6)}), [19] none 5 h 22 min

9  7434 t, x, u^{(10)}/u_{tx}^{(8)}, [23] Q(t, x, u^{(9)}/u_{tx}^{(7)}), [21] (34)  25 h

Table 4: The program ConLaw2 applied to compute conservation laws of the sine-Gordon equation.

4.5 Overview

The circumstance that the number of conservation laws for the sine-Gordon equation that were found by the different methods varies is due to the varying computational complexity of the determining equations they generate. The first three methods will find all conservation laws if memory and time requirements would not matter. The fourth approach is different in that it generates only necessary conditions which are often sufficient (if they have the same solution set as the other methods) but sometimes not. In that case the conditions (14) have additional adjoined symmetries as solutions. It may be that only specific linear combinations of them give a conservation law, as demonstrated with the example in the previous section. The following comments concentrate on complexity issues and other characteristic differences between the four approaches.

Arranging the methods as in the table below one can compare rows I,II and columns A,B.

A  B
I  Div P |_{Δ=0} = 0  ∑_{μ,J}(−D)_{J} (Q_{μ} \frac{∂Δ_{μ}}{∂u_{J}^{ν}} )|_{Δ=0} = 0 ∀ν

II  Div P = ∑_{ν} Q_{ν}Δ_{ν} ∑_{J}(−D)_{J} (\frac{∂}{∂u_{J}^{ν}} ∑_{μ} Q_{μ}Δ_{μ}) = 0 ∀ν

Table 5: The four approaches arranged in a table.

I-II: The conditions in row I are to be solved in the space of solutions (|_{Δ=0}), in row II they are not. This means that methods of row I can not be applied if equations or constraints Δ_{μ} = 0 can not be solved for a leading derivative but methods of row II can. Due to these substitutions the conditions in row I have fewer terms and involve fewer different derivatives of u than conditions in row II. The complexity of conditions and the number of conservation laws up to some order obtained in row I depend on whether Δ_{μ} = 0 is used to substitute lower order u-derivatives by higher ones or higher ones by lower ones. There are two reasons for this.

1) Substitutions based on 0 = Δ in Q may give extra restrictions for Q. For example, determining Q for conservation laws of the Korteweg de Vries equation 0 = Δ = u_{t}−u_{xxx}−uu_{x} and restricting Q to be of 2^{nd} order, then a substitution u_{t} = u_{xxx} + uu_{x} would imply Q = Q(t, x, u, u_{x}, u_{xx}), whereas a substitution u_{xxx} = u_{t} − uu_{x} would not restrict Q = Q(t, x, u, u_{t}, u_{x}, u_{tt}, u_{tx}, u_{xx}).

2) If a lower u-derivative is substituted by higher ones using 0 = Δ in the conservation law conditions in row I then such substitutions may increase the order of u-derivatives in which the conservation law conditions have to be satisfied identically. By that the desired effect of lowering the number of u-derivatives in which the conditions have to be fulfilled identically is lost. For example, condition IB for Δ = u_{tt} − u_{xxt}^{2} is 0 = Q_{tt} + 2(Qu_{xxt})_{xxt} which includes up to 6^{th} order u-derivatives (if Q is not of higher than 3^{rd} order). By substituting u_{tt} = u_{xxt}^{2} the order would increase to seven.

Hence, substituting lower order u-derivatives by higher order u-derivatives gives more over determined conditions for a less general ansatz. Such conditions are easier to solve, which may allow higher orders of Q to be investigated. However, one then may miss conservation laws of some order in P or Q.

These aspects are not an issue in row II as no substitutions are made there. A-B: In column A the single first order conservation law condition itself is to be solved, and in column B the integrability conditions of column A, which result when the conserved current P is eliminated are to be solved. Conditions in column B involve as many equations as there are functions u^{μ} and they are of the same order as the highest derivatives of u^{μ} in Δ_{ν} = 0.

The methods in column A compute P and therefore also trivial conservation laws when determining the general solution of the determining equations. This may complicate the solution of the determining equations to some extend. After the determining equations are solved the trivial conservation laws are easily dropped. For method IA one simply checks whether Div P = 0 holds identically and in method IIA trivial conservation laws have zero integrating factors Q_{μ}. The general solution of the determining conditions in column B do not generate trivial conservation laws. Also, conditions in column B are more straight forward to solve, they can be separated with respect to many high order u-derivatives and yield highly over determined systems. The disadvantage of methods in column B is that already their formulation may exceed available computational resources. Another potential problem with using methods in column B is the following. If one or more linear PDEs from the over determined conditions remain unsolved (for example, when investigating the Burgers equation (17) then the heat equation remained unsolved (see equations (17),(18)) then the program will usually not be able to compute P^{i} from the Q_{μ}. A way out is to use methods IB or IIB to get Q_{μ} and to use that as input to get P^{i} from method IIA.

Differences between the approaches are amplified with problems that involve an increasing number of PDEs and an increasing number of independent variables. A recommendation for tackling a new single PDE / system of PDEs would be:

• If any ansatz for some or all of the P^{i} is to be made then the method IA has to be used. (An example is the question whether an easily integrable conservation law D_{x}P^{x} = 0 with P^{y} = 0 exists.)

• Otherwise try IIB first.

• If the conditions become too large to handle then try IB or even IIA or IA.

• If IIB or IB provide the Q_{μ} but not the P^{i} (for example because arbitrary functions appear or some conditions remain unsolved) then try IIA using the computed Q_{μ} as input.

5 The computer algebra programs

The names of computer algebra programs for the four approaches are: IA: ConLaw1,IB: ConLaw2, IIA: ConLaw3, IIB: ConLaw4. They and the program Crack for solving the over determined conditions are written in the computer algebra system Reduce (ConLaw1 through ConLaw4 by the author, Crack by the author and Andreas Brand). From the general solution individual conservation laws are extracted by picking one arbitrary constant or function and setting all other arbitrary constants or functions to zero. The problem is to find, whether all arbitrary constants and functions are independent or whether some can be dropped without loss of generality. This problem itself leads to an over determined system of conditions which in general is very over determined and easy to solve. A description of that method is given in [35] where also the computation of Q from P and P from Q is explained.

Compared to other computer algebra programs, the package Crack has a wide variety of techniques for solving over determined PDE-systems. This allows the following new features as compared with other computer programs, a list of which and a short description is given in [12]:

• In all four computer programs P as well as Q are computed.

• By solving systems of over determined differential equations it is possible to find conservation laws with non-polynomial, even non-rational P, Q.

• If memory requirements are not too high then the program will make a definite statement about the existence of conservation laws of a given order. In the majority of these cases the program will find the explicit form of the conservation law, otherwise it will return unsolved equations.

• It is possible to find conservation laws with an explicit dependence of P, Q on the independent variables.

• There is no limit on the number of DEs nor the number of independent variables to be investigated for conservation laws other than a limit through the complexity of computations. Although not demonstrated in this paper, the program is able to handle ordinary differential equations (ODEs) as well.

• It is possible to determine values of parameters in the DE such that conservation laws exist.

• For each of the four programs ConLaw1 through ConLaw4 an ansatz for P^{i} and/or Q^{μ} can be input to specify to some extend conservation laws to be calculated.

A program written by Göktaş and Hereman [8] for computing conservation laws of PDEs makes a polynomial ansatz for conservation laws and finds the coefficients in this ansatz by solving a linear algebraic system of equations. Compared with that, the programs ConLaw1 through ConLaw4 are able to find more general conservation laws and to make a definitive statement in case the order is not too high to complete the computations. On the other hand, the program of Göktaş and Hereman was later extended to handle differential-difference systems [9],[10],[11].

Before showing examples which highlight the special abilities of ConLaw1 through Con-Law4 a comment to the treatment of ODEs shall be made. Although all methods and programs are applicable equally well to ODEs, the form of the ansatz for the integrating factor Q or for the first integral P to be made will usually be different. An nth order ODE has always first integrals of order n − 1 and any arbitrary function of first integrals is a first integral as well. In order to obtain an over determined system of conditions, the ansatz for a first integral must not contain functions of all variables x, u, u′, …, d^{n−1}u/dx^{n−1} but, for example, a polynomial in d^{n−1}u/dx^{n−1} with arbitrary functions of x, u, u′, …, d^{n−2}u/dx^{n−2} as coefficients or any other combination of functions of less than n + 1 variables, see also [4] for more details. Special features of the ConLaw programs that are not available with other programs are demonstrated with the following examples.

Example: The Ito equations for two functions u = u(t, x), v = v(t, x) read [15]

u_{t} = u_{xxx} + 6uu_{x} + 2vv_{x}
v_{t} = 2(uv)_{x}.

The first 7 conservation laws calculated by the program ConLaw1 which in turn calls Crack to solve condition (10), have the following values of P^{t}:

u, v, u^{2} + v^{2}, u_{x}^{2} − 2u^{3} − 2uv^{2}, (4uv^{2} − v_{x}^{2})/v^{3},

u_{xx}^{2} − 10uu_{x}^{2} − 4vv_{x}u_{x} + 5u^{4} + 6u^{2}v^{2} + v^{4}, ((2vv_{xx} − 4uv^{2} − 3v_{x}^{2})^{2} + 16v^{6})/v^{7}.

(P^{x} is not shown due to its length. It could be computed easily from P^{t}.) What is interesting is that 2 of the 7 conservation laws have a non-polynomial expression for P^{t} and as far as the author knows these conservation laws have not been known so far.

Example: The following equations [16] describe low - frequency Alfvén waves propagating parallel to an external magnetic field in a relativistic electron-positron plasma [31]. Typical for them is the symmetry with respect to interchanging the two functions u = u(t, x), v = v(t, x) due to the same charge-to-mass ratio for both kinds of particles. The equations are

Δ_{1} = u_{t} + r_{x} = 0, with r = u(u^{2} + v^{2}) + u_{xx},
Δ_{2} = v_{t} + s_{x} = 0, with s = v(u^{2} + v^{2}) + v_{xx}.  (15)

The equations themselves have the form of conservation laws. We find the following additional ones:

4uΔ_{1} + 4vΔ_{2} = D_{t}[2(u^{2} + v^{2})] +
D_{x}[4uu_{xx} − 2u_{x}^{2} + 4vv_{xx} − 2v_{x}^{2} + 3(u^{2} + v^{2})^{2}]
4rΔ_{1} + 4sΔ_{2} = D_{t}[(u^{2} + v^{2})^{2} − 2u_{x}^{2} − 2v_{x}^{2}] +
D_{x} [4u_{t}u_{x} + 4v_{t}v_{x} + 2u_{xx}^{2} + 2v_{xx}^{2} + 4(u^{2} + v^{2})×
((3(u^{2} + v^{2})t − x)(uu_{t} + vv_{t}) + uu_{xx} + vv_{xx})]
4(xu − 3tr)Δ_{1} + 4(xv − 3ts)Δ_{2} =
D_{t} [3t((2u_{x}^{2} + 2v_{x}^{2} − (u^{2} + v^{2})^{2}) + 2x(u^{2} + v^{2})] +
D_{x}2 [(uu_{t} + vv_{t}) (−x^{2} + (u^{2} + v^{2}) (6tx − 9t^{2}(u^{2} + v^{2})))
−3t(u^{2} + v^{2})^{3} + 3x(u^{2} + v^{2})^{2} + 2x(uu_{xx} + vv_{xx}) − 3tr^{2}
−3ts^{2} − 2uu_{x} − 2vv_{x} − xu_{x}^{2} − xv_{x}^{2} − 6tu_{t}u_{x} − 6tv_{t}v_{x}]

Whereas the first two are known [31], the last one shows an explicit x, t-dependence and is new. Further investigation provides that no conservation laws exist with the characteristic functions Q_{μ} of 3^{rd} or 4^{th} order (if u_{t}, v_{t} are substituted using (15)).

Example: The following equation of Gibbons and Tsarev [7]

0 = u_{xx} + u_{y}u_{xy} − u_{x}u_{yy} + 1  (16)

for u = u(x, y) is unusual in that it has already 5 conservation laws of first order. Characteristic functions contain x, y explicitly. Up to first order they are:

1, u_{y}, 3u_{y}^{2} + 2u_{x} + 3x, 2u_{y}^{3} + 3u_{x}u_{y} + 4u_{y}x + y,
10u_{y}^{4} + 6u_{x}^{2} + 24u_{x}u_{y}^{2} + 20u_{x}x + 30u_{y}^{2}x + 12u_{y}y + 2u + 15x^{2},

3u_{y}^{5} + 6u_{x}^{2}u_{y} + 10u_{x}u_{y}^{3} + 18u_{x}u_{y}x + 4u_{x}y + 12u_{y}^{3}x + 6u_{y}^{2}y + 12u_{y}x^{2} + 2u_{y}u + 6xy.

Example: The Liouville equation for a function u = u(x, y) reads

Δ = u_{xy} − e^{u}.

Conservation laws of order zero found by ConLaw2 are

(f_{x} + fu_{x})Δ = D_{x}(−e^{u}f) + D_{y}(f_{x}u_{x} + fu_{x}^{2}/2), f = f(x) arbitrary
(g_{y} + gu_{y})Δ = D_{y}(−e^{u}g) + D_{x}(g_{y}u_{y} + gu_{y}^{2}/2), g = g(y) arbitrary.

Because the ansatz made is investigated in full generality, any free functions in the conservation law will be found if the conditions can be solved completely by Crack. Otherwise the remaining conditions are returned as in the following example.

Example: The Burgers equation in the form

Δ = u_{t} − u_{xx} − \frac{1}{2}u_{x}^{2} = 0, u = u(t, x)  (17)

has zeroth order conservation laws

fe^{u/2}Δ = D_{t}(2fe^{u/2}) + D_{x}(e^{u/2}(2f_{x} − fu_{x}))  (18)

with f = f(t, x) satisfying the linear reverse heat equation 0 = f_{t} + f_{xx}.^{8}

The occurrence of free functions in the conservation law indicates linearizability of Δ = 0, which is the case for both previous examples. The following example involves more than 2 variables.

Example: The Kadomtsev-Petviashvili equation for u = u(t, x, y) with the abbreviation

w = u_{t} + 2uu_{x} + u_{xxx}

is

0 = Δ = w_{x} − u_{yy}.  (19)

Its zeroth order conservation laws include an arbitrary function c = c(t):

cΔ = D_{x}(cw) + D_{y}(−cu_{y})  (20)
cyΔ = D_{x}(cyw) + D_{y}(cu − cyu_{y})  (21)
(2cx + c_{t}y^{2})Δ = D_{t}(−2cu) +
D_{x} ((2cx + c_{t}y^{2})w − 2cu_{xx} − 2cu^{2}) +  (22)
D_{y} (−(2cx + c_{t}y^{2})u_{y} + 2c_{t}uy)
(6cxy + c_{t}y^{3})Δ = D_{t}(−6cyu) +
D_{x} ((6cxy + c_{t}y^{3})w − 6cyu_{xx} − 6cyu^{2}) +  (23)
D_{y} (−(6cxy + c_{t}y^{3})u_{y} + 3c_{t}uy^{2} + 6cxu) .

It is somewhat remarkable that although equation (19) does not involve u_{t} but only u_{xt} nevertheless the conserved density P^{t} in the last two conservation laws involves u and not u_{x}.

In the following section we give examples for an extension of our method to compute nonlocal conservation laws and report on the possibility to determine parameters in the equation such that conservation laws exist.

8Although already used in [2], [35] this example is shown again as it also serves to demonstrate an extension to non-local conservation laws in section 6.

6 Extending applicability

6.1 Non-local conservation laws

The implementations of the four methods have a common limitation: the characteristic functions Q and the conserved current P must depend functionally only on a finite number of derivatives of the u^{α}. No dependencies on integrals are possible. The same restriction is usually made when generators of Lie-symmetries are determined for differential equations. Whereas this restriction is less severe when calculating symmetries of PDEs, it is a serious restriction for the determination of conservation laws. To give an example, the Burgers equation in the form

Δ = u_{t} − u_{xx} − uu_{x} = 0, u = u(t, x)  (24)

has as low order conservation law only the trivial one D_{t}u − D_{x}(u_{x} + u^{2}/2) = 0. In order to include dependencies on ∫ u dx one could set u = v_{x} for some function v(x, t) and investigate conservation laws depending on v and derivatives of v. For the Burgers equation such a substitution alone is not enough. In addition one has to realize that (24) can be integrated with respect to x to f(t)_{t} = v_{t} − v_{xx} − v_{x}^{2} for some function f = f(t). Renaming v − f → u gives (17) and its conservation laws (18).

To give a further example, we consider the Boussinesq equation describing surface water waves whose horizontal scale is much larger than the depth of the water [1],[13]

u_{tt} − u_{xx} + 3uu_{xx} + 3u_{x}^{2} + αu_{xxxx} = 0.  (25)

Calculating conservation laws, using (25) to substitute u_{xxxx}, the only characteristic functions Q up to 4^{th} order are 1, x, t, xt. On the other hand, substituting u = v_{x}, integrating (25) with respect to x and renaming v − f → v gives

v_{tt} − v_{xx} + 3v_{x}v_{xx} + αv_{xxxx} = 0  (26)

having 2 conservation laws with characteristic functions 1, t which x-differentiated give the conservation laws above with characteristic functions 1, t. In addition two new conservation laws with characteristic functions v_{x}, v_{t} result. Repeating this step again: v = w_{x}, x-integration of (26), w − f → w gives

w_{tt} − w_{xx} + 3/2w_{xx}^{2} + αw_{xxxx} = 0  (27)

with three third order conservation laws. Two of them have characteristic functions w_{xxx}, w_{xxt} which correspond to the above conservation laws with characteristic functions v_{x}, v_{t}. In addition one extra conservation law with Q = w_{txx} − w_{txx}w_{xx} + w_{tx}w_{xxx} − \frac{2}{3}w_{ttt} exists.

A third example is the Kadomtsev-Petviashvili equation already discussed above^{9}. After a substitution u = v_{x}, x-integration of (19) and v − f → v the equation is

0 = [v_{t} + v_{xxx} + v_{x}^{2}]_{x} − v_{yy}.

Apart from conservation laws equivalent to (20),(21) three new conservation laws result with characteristic functions

−c_{tt}y^{2} − 2c_{t}x + 4cv_{x}
−c_{3t}y^{3} − 6c_{tt}xy + 12c_{t}yv_{x} + 24cv_{y},

9The hint to try KP for this extension was given by Alan Fordy.

−c_{4t}y^{4} − 12c_{3t}xy^{2} + 24c_{tt}y^{2}v_{x} − 12c_{tt}x^{2} + 48c_{t}xv_{x} + 96c_{t}yv_{y} + 48c_{t}v + 144cv_{t}.

Conserved currents are omitted due to their length. Repeating this transformation again does not yield conservation laws with characteristic functions of order less than three.

The purpose of this paragraph was to show that even if computer algebra programs ConLaw, Crack do only allow the investigation of local conservation laws depending on a finite number of derivatives of the unknown functions, we still may be able to enlarge the range of search by a contact transformation and integration of the PDE.

In the next section we extend the computation of conservation laws to the computation of parameters such that conservation laws exist.

6.2 Differential equations with parameters

In applications it is common that the DEs contain parameters and usually it would be desirable to know conservation laws which are valid for all possible values of these parameters. But as the example below shows, often conservation laws exist only for special values of parameters. Even if these parameter values are not of interest from the application side of view, the conservation laws valid for these values can at least be used, for example, to test numerical code. Another purpose for determining parameters together with conservation laws could be to find integrable equations from a more general class of equations.

The problem to determine parameters such that conservation laws exist is potentially much harder than determining conservation laws which are valid for any values of these parameters. This is because the conservation law determining equations become non-linear. Expressions may become unmanageably large and many sub cases may have to be considered. To use ConLaw1 through ConLaw4 for such calculations one only has to specify in its call the names of parameters to be computed (more details in the ConLaw manual).

Example: The 5^{th} order Korteweg - de Vries equation

u_{t} + αu^{2}u_{x} + βu_{x}u_{xx} + γuu_{3x} + u_{5x} = 0  (28)

with constant parameters α, β, γ includes well known special cases [6], [8], [14], [18], [28]: for α = 30, β = 20, γ = 10 the Lax equation [19], for α = 5, β = 5, γ = 5 an equation due to Sawata, Kotera [29] and Dodd and Gibbon [5], for α = 20, β = 25, γ = 10 an equation due to Kaup [17] and Kupershmidt, for α = 2, β = 6, γ = 3 an equation due to Ito [15].

The following zeroth and first order conservation laws are calculated with ConLaw1 (omitting P^{x} due to its length in the last two of these conservation laws):

• Q = 1, P^{t} = u

• α = β = γ = 0 : As (28) becomes linear, a conservation law is obtained with a characteristic function Q = Q(x, t) satisfying the adjoint PDE Q_{t} + Q_{5x} = 0 with P^{t} = Qu.

• α = 0, γ = β/3 : Q = x^{2}, P^{t} = x^{2}u

• α = 0, γ = β/3 : Q = x, P^{t} = xu

• γ = β/2 : Q = 2u, P^{t} = u^{2} • α = \frac{1}{10}(−2β^{2} + 7βγ − 3γ^{2}) :

Q = 60u_{xx}t(β − 3γ) + 6u^{2}t(2β^{2} − 7βγ + 3γ^{2}) + 60x P^{t} = 30u_{x}^{2}t(−β + 3γ) + u^{3}t(4β^{2} − 14βγ + 6γ^{2}) + 60ux

• α = \frac{1}{10}(−2β^{2} + 7βγ − 3γ^{2}) : Q = 30u_{xx} + 3u^{2}(2β − γ), P^{t} = −15u_{x}^{2} + u^{3}(2β − γ)

We find the same conservation laws as found by the program of Göktaş and Hereman and in addition a few conservation laws with explicit x, t-dependence.

7 Summary

Four approaches to find conservation laws have been compared with respect to their complexity and other characteristic features.

In a number of examples, conservation laws have been given, some of them new, which show that the programs ConLaw1,...,ConLaw4 and Crack can be used to find local, not necessarily polynomial, conservation laws with explicit variable dependence and free functions. The programs are, in principle, applicable to problems with arbitrarily many equations, functions and variables.

The programs including a manual and a test file are available via ftp from ftp.maths.qmw.ac.uk, directory pub/tw/crack. A demo web page which allows the use of ConLaw for problems of restricted size, is accessible from

http://cathode.maths.qmw.ac.uk/ConLaw demo.html. The package is part of the Re-

duce network library.

8 Acknowledgement

The author wants to thank Stephen Anco and especially Alexander Mikhailov for many comments. Willy Hereman and George Labahn are thanked for their comments on the manuscript. The Symbolic Computation Group at the University of Waterloo is thanked for its hospitality during the authors sabbatical in the fall of 1998.

References

[1] Ablowitz, M.J., Clarkson, P.A. (1991) Solitons, Nonlinear Evolution Equations and Inverse Scattering. London Math. Soc. Lec. Note Ser. 149. London: Cambridge University Press.

[2] Anco, S.C. and Bluman, G. (1997). Direct Construction of Conservation Laws from Field Equations Phys. Rev. Let. 78, no 15, 2869-2873.

[3] Anco, S.C. and Bluman, G. (1998). Direct construction method for conservation laws of differential equations, preprint.

[4] Anco, S.C. and Bluman, G. (1998). Integrating factors and first integrals for ordinary differential equations, Eur. J. Applied Math 9, p245-259.

[5] Dodd, R.K., Gibbon, J.D. (1977) The prolongation structure of a higher order Korteweg - de Vries equation, Proc. R. Soc. Lond. A 358, 287-296.

[6] Fordy, A., Gibbons, J. (1980) Some remarkable nonlinear transformations Phys. Lett. 75A 325.

[7] Gibbons, J. and Tsarev, S.,P. (1996) Reductions of the Benney Equations Phys. Lett. A 211, no 1, 19-24.

[8] Göktaş, Ü. and Hereman, W. (1997) Symbolic Computation of Conserved Densities for Systems of Nonlinear Evolution Equations. Journal of Symbolic Computation, 24, pp. 591-621.

[9] Göktaş, Ü., Hereman, W. and Erdmann, G. (1997) Computation of conserved densities for systems of nonlinear differential-difference equations, Physics Letters A, 236, pp. 30-38.

[10] Göktaş, Ü. and Hereman, W. (1998) Computation of conserved densities for nonlinear lattices, Physica D, 123, pp. 425-436.

[11] Hereman, W., Göktaş, Ü., Colagrosso, M. and Miller, A. (1998) it Algorithmic integrability tests for nonlinear differential and lattice equations, Computer Physics Communications, 115, pp. 428-446.

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

[13] Hickernell, F.J. (1983) The evolution of large - horizontal - scale disturbances in marginally stable, inviscid, shear flows. Stud. Appl. Math. 69, 23-49.

[14] Hirota R., Ito, M. (1983). Resonance of solitons in one dimension J. Phys. Soc. Jpn. 52 744-748.

[15] Ito, M. (1980). An extension of nonlinear evolution equations of the K-dV (mK-dV) type to higher orders. J. Phys. Soc. Jpn. 49, 771-778.

[16] Karney, C.F.F., Sen, A., Chu, F.Y.F. (1978). The complex modified Korteweg - de Vries equation, a non-integrable evolution equation in A.R. Bishop & T. Schneider (Eds), Solitons and condensed matter physics, Springer, p 71-75.

[17] Kaup, D.J. (1980) On the inverse scattering problem for cubic eigenvalue problems of the class Ψ_{xxx} + 6QΨ_{x} + 6RΨ = λΨ Stud. Appl. Math. 62 189-216.

[18] Kupershmidt, B.A., Wilson, G. (1981) Modifying Lax equations and the second Hamiltonian Structure Invent. Math. 62 403-436.

[19] Lax, P.D. (1968) Integrals of nonlinear equations of evolution and solitary waves Commun. Pure Appl. Math. 21,467-490.

[20] LeVeque, R.J. (1992). Numerical methods for Conservation laws. Lec. in Math., Basel: Birkhäuser Verlag.

[21] Mikhailov, A.,V. (1979) On the Integrability of the two-dimensional Generalization of the Toda Lattice Prisma v ZhETF, 30, 443.

[22] Mikhailov, A.V., Shabat, A.B. and Yamilov, R.I. (1987) On an Extension of the Module of Invertible Transformations Soviet Mathematical Doklady (Soviet Mathematical Survey) , 295, 288-291.

[23] Mikhailov, A.V., Shabat, A.B. and Yamilov, R.I. (1988) Extension of the Module of Invertible Transformations and Classification of Integrable Systems Comm. in Math. Phys. 115, 1-19.

[24] Kodama, Y. and Mikhailov, A.V. (1996) Obstacles to Asymptotic Integrability in Algebraic aspects of Integrability, ed. I.M. Gelfand and A. Fokas, Birkhauser, Boston, Basel, Berlin, 173-204.

[25] Newell, A.C. (1983). The history of the soliton. J. Appl. Mech. 50, 1127-1137.

[26] Olver, P.J. (1983). Conservation laws and null divergences. Math. Proc. Camb. Phil. Soc., 94, 529-540.

[27] Olver, P.J. (1986). Applications of Lie Groups to Differential Equations. Grad. Texts in Math. Berlin: Springer-Verlag.

[28] Satsuma, J., Kaup, D.J. (1977). A Bäcklund transformation for a higher order Korteweg - De Vries equation J. Phys. Soc. Jpn. 43 692-697.

[29] Sawada, K., Kotera, T. (1974). A method for finding N-Soliton solutions of the K.d.V. equation and K.d.V.-like equation Prog. Theor. Phys. 51 1355-1367.

[30] Shabat, A., Zhiber, A. (1979) Doclady Acad. Nauk SSSR, 247.

[31] Verheest, F. (1996). Solitary Alfvén modes in relativistic electron-positron plasmas. Physics Letters A 213, 177-182.

[32] Willox, R., Hereman, W. and Verheest, F. (1995). Complete Integrability of a Modified Vector Derivative Nonlinear Schrödinger Equation. Physica Scripta 52, 21-26.

[33] Wolf, T. and Brand, A. (1992) The Computer Algebra Package CRACK for Investigating PDEs, manual + software in the REDUCE network library.

[34] Wolf, T. (1996) The program CRACK for solving PDEs in General Relativity in Relativity and Scientific Computing Eds. F.W.Hehl, R.A.Puntigam, H.Ruder, Springer, 241-257.

[35] Wolf, T., Brand, A., Mohammadzadeh, M., Computer algebra algorithms and routines for the computation of conservation laws and fixing of gauge in differential expressions, preprint (1998), accepted for publication in J. of Symb. Comp., 18 pages.

9 Appendix:

Conservation Laws of the sine-Gordon equation

In this appendix conservation laws for the sine-Gordon equation

u_{tx} − sin(u) = 0

are shown as they have been computed by ConLaw1-4 and as they are refered to in tables above. They are not new, we provide them only to illustrate computer results. Except for the first conservation law, for all the following there is an additional conservation law due to the x ↔ t symmetry. These results are further examples of the ability of the programs to compute non-polynomial conservation laws.

2(tu_{t} − xu_{x})(u_{tx} − sin(u)) = D_{t} [2 cos(u)t − u_{x}^{2}x] + D_{x} [−2 cos(u)x + u_{t}^{2}t]  (29)
2u_{t}(u_{tx} − sin(u)) = D_{t} [2 cos(u)] + D_{x} [u_{t}^{2}]  (30)
(8u_{3t} + 4u_{t}^{3})(u_{tx} − sin(u)) = D_{t} [4 cos(u)u_{t}^{2} + 8u_{tx}u_{tt} − 8u_{tt} sin(u)] + D_{x} [−4u_{tt}^{2} + u_{t}^{4}] (31)

2(−8u_{5t} − 20u_{3t}u_{t}^{2} − 20u_{tt}^{2}u_{t} − 3u_{t}^{5})(u_{tx} − sin(u))  (32)
= D_{t}2 [−8 cos(u)u_{3t}u_{t} + 4 cos(u)u_{tt}^{2} − 3 cos(u)u_{t}^{4} − 8u_{tx}u_{4t} − 20u_{tx}u_{tt}u_{t}^{2}
+8u_{4t} sin(u) + 8u_{3t}u_{ttx} + 12u_{tt}u_{t}^{2} sin(u)]
+ D_{x} [−8u_{3t}^{2} + 20u_{tt}^{2}u_{t}^{2} − u_{t}^{6}]

8(−16u_{7t} − 56u_{5t}u_{t}^{2} − 224u_{4t}u_{tt}u_{t} − 168u_{3t}^{2}u_{t} − 280u_{3t}u_{tt}^{2} − 70u_{3t}u_{t}^{4} − 140u_{tt}^{2}u_{t}^{3} − 5u_{t}^{7})

×(u_{tx} − sin(u))  (33)
= D_{t}8 [16u_{6t} sin(u) − 16u_{tx}u_{6t} − 16 cos(u)u_{5t}u_{t} + 16 cos(u)u_{4t}u_{tt} − 8 cos(u)u_{3t}^{2}

−40 cos(u)u_{3t}u_{t}^{3} − 20 cos(u)u_{tt}^{2}u_{t}^{2} − 5 cos(u)u_{t}^{6} − 56u_{tx}u_{4t}u_{t}^{2} − 112u_{tx}u_{3t}u_{tt}u_{t} −56u_{tx}u_{tt}^{3} − 70u_{tx}u_{tt}u_{t}^{4} + 16u_{5t}u_{ttx} − 16u_{4t}u_{3tx} + 40u_{4t}u_{t}^{2} sin(u) + 56u_{3t}u_{ttx}u_{t}^{2}

+160u_{3t}u_{tt}u_{t} sin(u) + 40u_{tt}^{3} sin(u) + 30u_{tt}u_{t}^{4} sin(u)]

+ D_{x} [64u_{4t}^{2} − 224u_{3t}^{2}u_{t}^{2} + 112u_{tt}^{4} + 280u_{tt}^{2}u_{t}^{4} − 5u_{t}^{8}]

2(−128u_{9t} − 576u_{7t}u_{t}^{2} − 3456u_{6t}u_{tt}u_{t} − 7296u_{5t}u_{3t}u_{t} − 6720u_{5t}u_{tt}^{2} − 1008u_{5t}u_{t}^{4} −4416u_{4t}^{2}u_{t} − 24192u_{4t}u_{3t}u_{tt} − 8064u_{4t}u_{tt}u_{t}^{3} − 5824u_{3t}^{3} − 6048u_{3t}^{2}u_{t}^{3} − 24864u_{3t}u_{tt}^{2}u_{t}^{2} −840u_{3t}u_{t}^{6} − 6384u_{tt}^{4}u_{t} − 2520u_{tt}^{2}u_{t}^{5} − 35u_{t}^{9})

×(u_{tx} − sin(u))  (34)
= D_{t}2 [128u_{8t} sin(u) − 128 cos(u)u_{7t}u_{t} + 128 cos(u)u_{6t}u_{tt} − 128 cos(u)u_{5t}u_{3t} + 64 cos(u)u_{4t}^{2}

−448 cos(u)u_{5t}u_{t}^{3} − 1344 cos(u)u_{4t}u_{tt}u_{t}^{2} − 1568 cos(u)u_{3t}^{2}u_{t}^{2} − 1344 cos(u)u_{3t}u_{tt}^{2}u_{t} −560 cos(u)u_{3t}u_{t}^{5} + 336 cos(u)u_{tt}^{4} − 840 cos(u)u_{tt}^{2}u_{t}^{4} − 35 cos(u)u_{t}^{8} − 1008u_{t,x}u_{4t}u_{t}^{4} −576u_{t,x}u_{6t}u_{t}^{2} − 2304u_{t,x}u_{5t}u_{tt}u_{t} − 4992u_{t,x}u_{4t}u_{3t}u_{t} − 4416u_{t,x}u_{4t}u_{tt}^{2} − 128u_{t,x}u_{8t} −5184u_{t,x}u_{3t}^{2}u_{tt} − 4032u_{t,x}u_{3t}u_{tt}u_{t}^{3} − 4256u_{t,x}u_{tt}^{3}u_{t}^{2} − 840u_{t,x}u_{tt}u_{t}^{6} + 128u_{7t}u_{tt,x} −128u_{6t}u_{3t,x} + 448u_{6t}u_{t}^{2} sin(u) + 128u_{5t}u_{4t,x} + 576u_{5t}u_{tt,x}u_{t}^{2} + 2688u_{5t}u_{tt}u_{t} sin(u)

−576u_{4t}u_{3t,x}u_{t}^{2} + 4480u_{4t}u_{3t}u_{t} sin(u) + 1152u_{4t}u_{tt,x}u_{tt}u_{t} + 4032u_{4t}u_{tt}^{2} sin(u) +560u_{4t}u_{t}^{4} sin(u) + 1920u_{3t}^{2}u_{tt,x}u_{t} + 5824u_{3t}^{2}u_{tt} sin(u) + 3264u_{3t}u_{tt,x}u_{tt}^{2}

+1008u_{3t}u_{tt,x}u_{t}^{4} + 4480u_{3t}u_{tt}u_{t}^{3} sin(u) + 3360u_{tt}^{3}u_{t}^{2} sin(u) + 280u_{tt}u_{t}^{6} sin(u)]

+ D_{x} [−128u_{5t}^{2} + 576u_{4t}^{2}u_{t}^{2} − 1280u_{3t}^{3}u_{t} − 3264u_{3t}^{2}u_{tt}^{2} − 1008u_{3t}^{2}u_{t}^{4} + 2128u_{tt}^{4}u_{t}^{2}

+840u_{tt}^{2}u_{t}^{6} − 7u_{t}^{10}]