Computer algebra algorithms and routines for the computation of conservation laws and fixing of gauge in differential expressions

Thomas Wolf

Queen Mary & Westfield College, University of London,

Mile End Road, London E1 4NS, UK email: T.Wolf@maths.qmw.ac.uk

Andreas Brand

Institut für Informatik, Friedrich Schiller Universität Jena email: maa@cnve.rz.uni-jena.de

Majid Mohammadzadeh

Faculty of Mathematics and Computer Engineering University of Teacher Education 49 Mofateh Ave, PB 15614, Tehran

email: majid@saba.tmu.ac.ir

November 19, 2001

Abstract In this paper we present three different approaches for the determination of conservation laws. For three corresponding REDUCE computer algebra programs CONLAW1/2/3 the necessary subroutines are described. All three programs use subroutines which remove redundant functions and constants in the general solution of the conservation law conditions. The corresponding algorithm is explained. Such a program is not only applicable in the context of computing conservation laws but whenever redundacy in differential expressions is to be removed or gauge freedom to be fixed.

1 Introduction

The determination of conservation laws for single or systems of partial differential equations (PDEs) and of first integrals for ordinary differential equations (ODEs) is of interest for the exact solution of these DEs, for their understanding and classification and for supporting their numerical solution. In this paper we outline three computer algebra programs implemented in Reduce for the computation of conservation laws. We also explain in detail a procedure that fixes redundancy of arbitrary functions and constants in differential expressions. This procedure is used to extract individual conservation laws from the general solution of conservation law determining equations.

In this paper we will not make the simplifying assumption that differential equations Δ = 0 for which we want to find conservation laws, result from a variational principle. In that case any variational symmetry of the Lagrangian provides a conservation law as is known from Noether’s Theorem. Instead, we will not make any restrictive assumptions which leaves us to solve Div P = 0 either directly or to determine characteristic functions of conservation laws or to do both at once. A more application oriented comparison of these different approaches will be described elsewhere [11]. In this paper we concentrate on the computer algebra aspects.

The following section discusses the mathematical background of the three approaches. This is followed by a section on the computer algebra side of computing conservation laws. In this section the need for fixing redundancy in general solutions of conservation law conditions is motivated. A corresponding algorithm is explained in section 4 together with other procedures that are more specific for the computation of conservation laws. In section 5 examples are given that demonstrate in which sense the programs CONLAW1/2/3 are superior to other computer algebra programs that compute conservation laws. In two further sections the three methods are compared with each other and the syntax of the three programs is given.

2 The mathematical problem and the three approaches

In the following we adopt the notation of the book of Olver [7]. Independent variables will be denoted by x = (x^{1}, x^{2}, …, x^{p}). The differential equations for which we seek conservation laws are Δ(x, u^{(n)}) = 0 (i.e. Δ_{1} = 0, …, Δ_{q} = 0), for q functions u = (u^{1}, u^{2}, …, u^{q}), u^{(n)} denoting u - derivatives of order up to n. The conservation law that is to be fulfilled by solutions of Δ = 0 is Div P = ∑_{i=1}^{p} D_{i}P^{i} = 0 for a current P = (P^{1}, …, P^{p}) and total derivatives D_{i}. In this paper local conservation laws are considered where components P^{i} are expressions in x^{j}, u^{(k)} and where Div P = 0 has to be satisfied identically in x^{j}, u^{(k)} modulo Δ = 0. We will use _{J} as a multiple index denoting partial derivatives, for example, u_{J}^{α} will stand for an arbitrary partial derivative, like ∂^{k}u^{α}/(∂x^{1}∂x^{2} …).

A first approach to compute conservation laws is to solve

Div P = 0 modulo Δ = 0  (1)

directly. The components of the conserved current P^{1}, …, P^{p} that are to be calculated are allowed to depend on all independent variables x^{i}, the dependent variables u^{α} and their derivatives u_{J}^{α} up to some order.

We are not interested in trivial conservation laws P = curl V but rather in conservation laws that are obeyed by the solutions of Δ = 0. Therefore we use Δ = 0 to eliminate some of the so-called jet-variables u_{J}^{α} and substitute them in the determining conditions (1). This implies that the condition (1) has to be fulfilled identically in fewer variables. Thus the condition becomes less restrictive and so may have additional solutions apart from P = curl V. These extra non-trivial conservation laws are the ones of interest. We therefore assume Δ = 0 can be solved for leading derivatives u_{J}^{α} so that these and all their partial derivatives that occur in (1) can be substituted.

The conserved current P^{i} is defined only modulo Δ = 0, because conservation laws are to be obeyed by solutions of Δ = 0. We therefore assume, without loss of generality, that the P^{i} do not depend on the u-derivatives that we substitute. This fixes another kind of equivalence of conservation laws.

Other approaches calculate characteristic functions Q^{ν}. In ([7], p. 272) it is shown that, for a totally non-degenerate system Δ_{ν} = 0, each equivalence class of conservation laws Div P = 0 (i.e. conserved currents differing only by a curl) is determined uniquely by characteristic functions Q^{ν} satisfying

Div P = ∑_{ν} Q^{ν}Δ_{ν}  (2)

identically in all x^{i}, u^{α}, u_{J}^{α}. Equation (2) is not solved by simply eliminating Q^{1} in terms of P and Δ and the other Q^{ν} as Q^{1} would be singular for solutions of Δ = 0. To avoid this problem and because the Q^{ν} are only unique modulo Δ = 0, we can ignore dependencies of Q^{ν} on leading u-derivatives in Δ = 0 and any of their derivatives. One way to calculate the Q^{ν} is to use the basic property of the Euler operators E_{ν} = ∑_{J}(−D)_{J}∂/∂(u_{J}^{ν}), namely to give zero if and only if they act on an expression which is a divergence.

Applying the Euler operator on (2) gives

0 = ∑_{μ,J}(−D)_{J} (\frac{∂(Q^{μ}Δ_{μ})}{∂(u_{J}^{ν})}) for all ν.  (3)

Using Δ_{ν} = 0 to simplefy the condition one obtains as necessary (but not sufficient) determining conditions for the Q^{ν}:

0 = ∑_{μ,J}(−D)_{J} (Q^{μ} \frac{∂Δ_{μ}}{∂(u_{J}^{ν})}) for all ν.  (4)

The three methods and programs discussed in this paper aim at solving (1),

(2) and (4) identically in x^{i}, u_{α}, u_{J}^{α}. The program CONLAW1 tries to solve (1) for P^{i}, CONLAW2 tries to solve (4) for Q^{ν} and finally CONLAW3 tries to solve (2) for P^{i} and Q^{ν}.

Any two of the three approaches (1)-(4) differ either in the number of conditional equations to be solved, the order of these equations, the number of functions to be determined, the number of independent jet-variables, or the degree of an ansatz for P, Q that is necessary in order to obtain the same conservation laws.

To obtain solutions of (1)-(4) we assume bounds on the order of uderivatives on which the P^{i} and Q^{ν} may depend. For (1) we assume a bound for P^{i} and for (2),(4) we assume a bound for Q^{ν}. Bounds for the remaining unknown functions follow. Differentiations done in all three conditions (1)-(4) introduce jet-variables (u-derivatives) on which none of the unknown functions depend.

As the conditions have to be satisfied identically in these variables as well, over determined PDE-systems result. Those are investigated with the computer algebra package CRACK ([9],[10]).

3 The computer algebra problem

The main computer algebra problem is to solve the over determined conditions (1)-(4). Steps for solving that problem include separation, integration, application of integrability conditions (differential Gröbner Basis), solution of ODEs and other steps, all of which are described in [9],[10]. The success rate of Crack in solving over determined systems is high if the system is linear (like (1)-(4) which are linear in P^{i}, Q^{ν}) and not too big. Examples of what is currently possible are given in section 5. Any unsolved equations are returned. For example, when conservation laws of the Burgers equation are investigated in section (5), the heat equations remains as an unsolved condition and is returned. Complexity issues are discussed in [11].

From the general solution of (1),(2) or (4) a single conservation law is extracted by collecting all terms involving one of the arbitrary constants or arbitrary functions in the solution. Any redundant constants would give conservation laws which are not independent of each other.

Redundant constants and functions may result because in the process of solving the over determined system, there is no general rule for deciding whether integrations or substitutions should have a higher priority. There are examples requiring a higher priority for each of them. It therefore may happen that two equations are integrated which are not independent of each other and therefore the constants or functions of integration are not independent of each other. As a consequence, the final general solution could have redundant arbitrary constants and functions. For example, in the expression c_{1}(x)t + c_{2}xt + c_{3}, with independent variables x, t and arbitrary function c_{1}(x) and arbitrary constants c_{2}, c_{3}, the constant c_{2} is redundant as it can be absorbed by c_{1}(x) through c_{1}(x) → c_{1}(x) − c_{2}x. In this paper we will refer to redundancies also as gauge freedom.

Another application of redundancy recognition is the solution of PDE systems with some gauge freedom where the problem is to eliminate any gauge freedom from the general solution of this system. This can be accomplished by including in the solution terms representing the complete gauge freedom. For example, in the case of conditions (1)

• the general solution could be augmented by curl V and V be added to the list of free constants and functions,

• redundant constants and functions could be dropped, i.e. constants and functions that could be absorbed by V,

• finally V could be dropped from the list of free functions.

In this way trivial conservation laws could be filtered out as the free constants and functions corresponding to such conservation laws would have been absorbed by V.

In the case of computing conservation laws, one could easily drop trivial conservation laws after they have been computed by checking Div P = 0 identically in all jet-variables. Such a simple test to eliminate gauge freedom might not be available for other problems. Therefore a general algorithm is described below.

4 Subroutines

In the following subsections we describe an algorithm for fixing redundancy in the general solution of a PDE-problem (and by that extracting conservation laws from the general solution of conditions (1)-(4)), subroutines that compute Q^{ν} from P^{i} and P^{i} from Q^{ν}, and subroutines that simplify P^{i}.

4.1 Identifying redundant constants and functions in solutions

This and the following section 4.2 are self contained and could be read independently of the rest of the paper. In them we consider solutions of arbitrary differential equations, having in mind the conditions (1)-(4) (and less the equations Δ = 0, which we can not solve and which therefore are investigated for conservation laws). To avoid confusion we call these equations Ω(f_{a}, x^{i}) = 0 which are to be solved for the functions f_{a}(x^{i}). In the case of computing conservation laws, the equations Ω = 0 are the equations (1)-(4), the f_{a} are the P^{i} and/or Q^{μ}, and the x^{i} now include all jet-variables. This notation is only used in the sub-sections 4.1 and 4.2.

In this sub-section we want to discuss the problem of fixing gauge freedom, if only a PDE-solver is available that can produce solutions with possible redundancies. Before giving the steps in a summary, the algorithm is explained using simple examples.

Given is a system of DEs Ω(f_{a}, x^{i}) = 0 that is to be solved for the functions f_{a}(x^{i}). We assume

f_{b} = F_{b}(x^{i}, g_{c})  (5)

to be a general solution where F_{b} are differential expressions in x^{i}, g_{c} where g_{c} are arbitrary constants and functions of x^{i} (which, for example, resulted through integration).

The question is how to specify as many as possible of the g_{c} to fix any redundancy among them but not to loose generality of the solution. The idea behind the algorithm is the following.

If two sets of constants and functions g_{c} and ḡ_{c} in the solution (5) differ only by a gauge then for both sets the solution (5) would still be the same identically in the x^{i}. Therefore the determining conditions for the redundancy will be the invariance condition for the general solution, i.e.

F_{b}(x^{i}, g_{c}) − F_{b}(x^{i}, ḡ_{c}) = 0  (6)

identically in x^{i} where ḡ_{c} are new functions, each ḡ_{a} having the same variable dependence as the g_{a} with the same lower index.

Example: We assume the general solution of some problem for a single function f(x^{1}, x^{2}, x^{3}) being computed as

f = x^{1}g_{1}(x^{2}) + g_{2}(x^{1}, x^{2}) + x^{3}g_{3}  (7)

with arbitrary functions g_{1}(x^{2}), g_{2}(x^{1}, x^{2}) and arbitrary constant g_{3}.^{1} For the solution (7), the invariance condition (6) is

x^{1}g_{1}(x^{2}) + g_{2}(x^{1}, x^{2}) + x^{3}g_{3} − x^{1}ḡ_{1}(x^{2}) − ḡ_{2}(x^{1}, x^{2}) − x^{3}ḡ_{3} = 0. (8)

Seen as a condition for functions g_{1}(x^{2}), ḡ_{1}(x^{2}), g_{2}(x^{1}, x^{2}), ḡ_{2}(x^{1}, x^{2}) and constants g_{3}, ḡ_{3} it is solved by separation of variables (here x^{3}):

0 = g_{3} − ḡ_{3},
0 = x^{1}g_{1}(x^{2}) + g_{2}(x^{1}, x^{2}) − x^{1}ḡ_{1}(x^{2}) − ḡ_{2}(x^{1}, x^{2})

which has as a general solution

ḡ_{3} = g_{3}
ḡ_{2}(x^{1}, x^{2}) = g_{2}(x^{1}, x^{2}) + xg_{1}(x^{2}) − xḡ_{1}(x^{2}),  (9)

ḡ_{1}(x^{2}), g_{1}(x^{2}), g_{2}(x^{1}, x^{2}), g_{3} are free. If it is possible to compute the general solution of the invariance condition (6) then all g_{c} for which g_{c} itself or ḡ_{c} are computed (occur on the left hand side of this general solution) are essential free constants and functions of the general solution (5). In the solution (7) the g_{2}(x^{1}, x^{2}), g_{3} are essential.

The remaining g_{c} which occur together with their counterpart ḡ_{c} on right hand sides of the general solution of the invariance condition are redundant. They can be chosen freely (like g_{1}(x^{2}) in solution (7) of the example). Any change in them can be compensated by a change in the essential functions, as (9) shows how ḡ_{2}(x^{1}, x^{2}) has to be adapted if g_{1}(x^{2}), ḡ_{1}(x^{2}) would be changed.

Although the solvability of (6) for g_{a}, ḡ_{b} cannot be guaranteed, this should in practice not be a problem for the following reasons.

• Under normal circumstances when (5) is a solution to some PDE(s), functions g_{c} appearing on the right hand side have to be functions of less variables than f_{b} on the left hand side. The invariance condition (6) which is to be satisfied in all independent variables x^{i} identically is therefore a highly over determined problem for functions g_{c}, ḡ_{c} of only a subset of variables of (6). Even if the problem to which (5) is the solution is underdetermined, and so consequently some g_{c} have the same variable dependence as the f_{b}, still their number must be smaller than the number of f_{b}. In any case, the invariance condition (6) is a over determined set of conditions for the g_{c}, ḡ_{c} and therefore well suited to be solved by the package CRACK.

1The general solution to conservation law problems usually involves only constants, like in examples 1,3 in section 5. Example 2 in section 5 only involves functions because this is a rare case of a linearizable equation - the Burgers equation.

• Although it is over determined, the invariance condition (6) has always at least one trivial solution ḡ_{c} = g_{c}, for all c. This is the general solution of the invariance condition if no redundancy exists.

• If the equations Ω = 0 are linear in f_{a}, then their solution is linear in the arbitrary functions g_{c}, i.e. equation (6) is linear in the g_{c}, ḡ_{c}. This is the case for the computation of conservation laws.^{2}

• If the equations Ω = 0 are non-linear in f_{a}, then solving (6) should still be easier than the solving Ω = 0, an equation which we assumed was possible to solve.

In the following example we show how new redundancy can result from solving the invariance condition. The example also illustrates that it may happen that one is able to compute gauge freedom but that it is impossible to remove it.

Example: As mentioned above, the more underdetermined a problem is, the more difficult it is in general to fix redundancy in its general solution. In this respect the following well known simple problem (one condition for 3 unknowns) is comparatively hard. Equation

0 = Divf = ∂_{1}f_{1} + ∂_{2}f_{2} + ∂_{3}f_{3}  (10)

for functions f_{1}, f_{2}, f_{3} of variables x^{1}, x^{2}, x^{3} has the general solution

f_{1} = ∂_{3}g_{2} − ∂_{2}g_{3}
f_{2} = ∂_{1}g_{3} − ∂_{3}g_{1}  (11)
f_{3} = ∂_{2}g_{1} − ∂_{1}g_{2}

where g_{1}, g_{2}, g_{3} are arbitrary functions of x^{1}, x^{2}, x^{3}. (Equation (10) is not a conservation law problem. It is not to be solved modulo some equation, but describes a static source free field.)

2Conditions become non-linear if Δ = 0 does contain parameters and we want to calculate their value such that conservation laws exist.

The general solution (11) must have gauge freedom as it has the same number of free functions g_{i} of all variables as there are functions f_{i} to be determined initially. To solve the invariance conditions

0 = ∂_{3}g_{2} − ∂_{2}g_{3} − ∂_{3}ḡ_{2} + ∂_{2}ḡ_{3}  (12)
0 = ∂_{1}g_{3} − ∂_{3}g_{1} − ∂_{1}ḡ_{3} + ∂_{3}ḡ_{1}  (13)
0 = ∂_{2}g_{1} − ∂_{1}g_{2} − ∂_{2}ḡ_{1} + ∂_{1}ḡ_{2}  (14)

a computer algebra program could do the following steps:

• Solve (12) by introducing a new function k_{1}(x^{1}, x^{2}, x^{3}) with

ḡ_{2} − g_{2} = ∂_{2}k_{1}, ḡ_{3} − g_{3} = ∂_{3}k_{1}.  (15)

• Solve (13) by substituting g_{3} from (15) and integrate it with respect to x^{3} with a new function k_{2}(x^{1}, x^{2}) to get

ḡ_{1} − g_{1} = ∂_{1}k_{1} + k_{2}.

• equation (14) finally requires ∂_{2}k_{2} = 0, giving the general solution

ḡ_{1} = g_{1} + ∂_{1}k_{1} + k_{2}  (16)
ḡ_{2} = g_{2} + ∂_{2}k_{1}  (17)
ḡ_{3} = g_{3} + ∂_{3}k_{1}  (18)

with k_{1} = k_{1}(x^{1}, x^{2}, x^{3}), k_{2} = k_{2}(x^{1}).

This example shows how solving the invariance condition may require solving differential equations, therefore the solution may involve extra redundant functions, here k_{2}. The free function k_{1} is a function of all variables. Therefore, if k_{1} occurs in one of equation (16)-(18) only algebraically then at least locally this equation could be solved for k_{1}, and k_{1} could be substituted so that only two of the g_{c} are essential.

The situation in solution (16)-(18) is different. k_{1} can not be solved for and the redundancy can not be fixed. This deficiency is not only a technical aspect. If a derivative of k_{1} could be renamed to k̃_{1} in order to have k̃_{1} appearing algebraically in at least one equation and everywhere else only k̃_{1} (not k_{1}) or derivatives of k̃_{1} then this would be done in the corresponding routine in CRACK. However, in (16)-(18) k_{1} comes with different derivatives, therefore it is not a deficiency of the algorithm or computer program that gauge freedom can not be fixed explicitly in a general solution like (11).

The following summary is a list of steps as they are performed in the program CRACK to fix gauge freedom in a given solution (5)

f_{b} = F_{b}(x^{i}, g_{c})

to some differential problem:

• Formulate a set of conditions (6)

F_{b}(x^{i}, g_{c}) − F_{b}(x^{i}, ḡ_{c}) = 0

where ḡ_{c} are new functions, each ḡ_{c} having the same variable dependence as g_{c}. Regard equation (6) as a system of equations for the set of unknown functions {g_{c}, ḡ_{c}}, to be satisfied identically in the x^{i}.

• Find the general solution of the system (6) as

g̃_{c} = G_{c}(x^{i}, h_{d})  (19)

where g̃_{c} is a subset of {g_{a}, ḡ_{b}}, and G_{c} are algebraic or differential expressions in functions h_{d} which are the remaining {g_{a}, ḡ_{b}} and extra constants and functions of integration. (In the example above the G_{c} are the right hand sides of (16)-(18), the g̃_{c} are ḡ_{1}, ḡ_{2}, ḡ_{3} and the h_{d} are g_{1}, g_{2}, g_{3}, k_{1}, k_{2}.) The h_{d} are arbitrary. Any function g_{a} or ḡ_{a} appears only once on a left hand side of (19) or only on right hand sides.

• If for any index c, both g_{c} and ḡ_{c} appear only on right hand sides of (19) then g_{c} is redundant and can be set to zero in all F_{b} in (5) and all G_{c} in (19).

• The system (19) will now be solved for the ḡ_{c}: if for any index c, both g_{c} and ḡ_{c} appear only on left hand sides of (19) in the equations g_{c} = G_{c} and ḡ_{c} = Ḡ_{c} then these two equations are replaced by ḡ_{c} = g_{c} −G_{c} +Ḡ_{c} in (19).

• If for any index c, g_{c} appears on a left hand side of (19) and ḡ_{c} appears only on right hand sides then the equation with left hand side g_{c} is solved for ḡ_{c} and replaced by the new equation ḡ_{c} = Ḡ_{c}(g_{c}, …). With this new equation ḡ_{c} is substituted on any right hand side of (19).

• There remains only the case of ḡ_{c} being on the left hand side of an equation and g_{c} being on right hand sides. The system (19) now has the form

ḡ_{c} = Ḡ_{c}(x^{i}, g_{a}, h̄_{b})  (20)

where h̄_{b} are arbitrary constants and functions of integration which arose during the solution of (6). ḡ_{c} do not occur on right hand sides as they would be redundant and would have been set to zero otherwise.

• Free constants and functions h̄_{b} that occur only algebraically on the right hand side of at least one equation and that depend on all independent variables in that equation, will be chosen to make as many Ḡ_{c} as possible zero. In such a case also the redundant g_{c} is set to zero in (5) and (20). (In equation (16) k_{2} could not be used as k_{2} does not depend on all variables and k_{1} could not be used as it does not occur algebraically in at least one equation, also not through the substitution of a derivative of k_{1} by k_{1} everywhere.) As we do not have to know h̄_{b} explicitly, it is enough to find equations in (20) which include such a h̄_{b}. Assuming local solvability of Ḡ_{c} = 0 for h̄_{b} we conclude redundancy of g_{c}.

• Finally, all remaining h̄_{b} which cannot be used to make a right hand side zero are set to zero themselves (in example (16)-(17) these are k_{1}, k_{2}). The final form of (20)

ḡ_{c} = Ḡ_{c}(x^{i}, g_{a})  (21)

provides substitutions which turn the old solution (5) (f_{b} = F_{b}(x^{i}, ḡ_{c}) after g_{c} have been renamed ḡ_{c}) into a new solution f_{b} = F̃_{b}(x^{i}, g_{c}) for which as much redundancy as possible is fixed.

4.2 Identifying redundant constants and functions in solutions and remaining unsolved equations

The above steps for fixing redundancy are not only applicable once a general solution of a PDE(-system) Ω(f_{a}, x^{i}) = 0 has already been found. Also, if a preliminary solution is known but some equations remain to be solved, then fixing at least some redundancy is possible. For example, the computation of conservation laws for the Burgers equation in section 5 returns a solution together with the heat equation which remains unsolved.

What could the fixing of redundancies mean under these circumstances? Substitutions that find redundancies in g_{c} in a preliminary solution f_{b} = F_{b}(x^{i}, g_{c}) and in remaining differential equations D(x^{i}, g_{c}) = 0 which have to be satisfied by g_{c}, should leave the solution completely unchanged but should be allowed to change the remaining equations D = 0 into a set of equivalent equations. This is in general a much harder problem than we had before as it involves equivalence of PDE-systems. A simpler question is to find redundancy fixing substitutions that leave the solution and equations unchanged. This is done by extending redundancy conditions (6) by D(x^{i}, g_{c}) − D(x^{i}, ḡ_{c}) = 0.

Although this simplified version of fixing redundancies seems to potentially miss redundancies, we have not encountered such a case yet, given that our remaining equations are at most polynomially non-linear (for conservation law calculations linear), unknown functions are not nested, or depending on expressions (like f(x+y)), and in each equation there is at least one function of all variables.

Example: Using the above steps to delete redundancies in the preliminary solution

f = g_{1}(x^{1}, x^{3}) − g_{2}(x^{1}, x^{2}) − g_{3}(x^{2})  (22)

and the remaining condition

g_{2}(x^{1}, x^{2}) + g_{3}(x^{2}) + g_{4}(x^{3}) = 0  (23)

would show that g_{3} is redundant but g_{4} not. The reason for not recognizing that after substituting g_{2} from (23) in (22), g_{4} can be absorbed by g_{1}, is that equation (23) does not involve a function of all variables x^{1}, x^{2}, x^{3}.

A situation like this can not occur in practice, as the program CRACK has routines to handle equations in which there is no function of all variables. An equation like (23) can therefore not remain unsolved. CRACK would find g_{4} =const. (and g_{2} + g_{3} =const.) and then g_{4} would be absorbed.

To look at the problem in terms of redundancy fixing substitutions that transform (23) into an equivalent condition, a computer program would have to be able to see that equation (23) is equivalent to

g_{2}(x^{1}, x^{2}) + g_{3}(x^{2}) + g_{5} = 0, and g_{4}(x^{3}) = g_{5} = const.

and then g_{5} would be found to be redundant.

The possibility of fixing gauge freedom even in the presence of yet unsolved equations opens the possibility of running a gauge-fixing step during the process of solving over determined PDE-systems (which is currently not done). By that the number of unknown functions could be reduced and the remaining equations simplified.

4.3 Computing characteristic functions from conserved currents

The remainder of the paper is concerned with the computation of conservation laws. In this sub-section we explain how to compute the characteristic functions Q from the conserved current P.

The first approach (1) to compute conservation laws is attractive compared with (2),(4) as it generates only one PDE to be solved which is of first order and involves fewer jet-variables than approach (2) because it is computed modulo Δ = 0. Also, it has fewer functions to compute than approach (2). A negative aspect is that it provides only the conserved current P and not the characteristic functions Q.

If expressions Q^{νJ} in a relation (24) below are known then partial integrations (25) yield the characteristic functions Q^{ν} and the corresponding conserved current P − R:

Div P = 0 modulo Δ_{ν} = 0 ↔
∃ Q^{νJ} : Div P = ∑_{ν,J} Q^{νJ}Δ_{ν}^{(J)}  (identically in all x, u_{J}^{α}) (24)
= ∑_{ν,J} D_{J}(Q^{νJ}Δ_{ν}) − D_{J}(Q^{νJ})Δ_{ν}  (repeatedly)(25)
= Div R + ∑_{ν} Q^{ν}Δ_{ν}

Equation (24) cannot be regarded as a linear algebraic equation to determine Q^{νJ} as there is the additional requirement that the Q^{νJ} are non-singular for solutions of Δ = 0. Instead, Div P is calculated and substitutions of a different form from before are made. For example, if conservation laws for the Harry Dym equation 0 = Δ = u_{t} − u^{3}u_{xxx} are investigated and if for the derivation of (1) there had been done substitutions u_{t} = u^{3}u_{xxx}, u_{tx} = (u^{3}u_{xxx})_{x}, … before then now the substitutions would be u_{t} = Δ + u^{3}u_{xxx}, u_{tx} = Δ_{x} + (u^{3}u_{xxx})_{x}, … which provide the right hand side of (24). The computation of Q^{ν} and P^{i} − R^{i} from P^{i} is part of CONLAW1.

4.4 Computing conserved currents from characteristic functions

The inverse computation is necessary in CONLAW2 where the conserved current P^{i} has to be computed from Q^{μ} by integrating Div P = ∑_{ν} Q^{ν}Δ_{ν}.

A direct way is based on a formula given by Anco & Bluman in [1]:

P^{i} = ∫_{0}^{1} \frac{dλ}{λ} (S^{i}(u) + N_{μ}^{i}(u)u^{μ} + N_{μ}^{ij}(u)D_{j}u^{μ} + …) |_{u→λu} (26)
S^{i}(u) = Q^{ν} \frac{∂Δ_{ν}}{∂u_{i}^{μ}} u^{μ} + Q^{ν} \frac{∂Δ_{ν}}{∂u_{ij}^{μ}} u_{j}^{μ} − (Q^{ν} \frac{∂Δ_{ν}}{∂u_{ij}^{μ}} )_{j} u^{μ} + …  (27)

N_{μ}^{i}(u) = \frac{∂Q^{ν}}{∂u_{i}^{μ}} Δ_{ν} − (\frac{∂Q^{ν}}{∂u_{ij}^{μ}} Δ_{ν})_{j} + (\frac{∂Q^{ν}}{∂u_{ijk}^{μ}} Δ_{ν})_{jk} − …  (28)

N_{μ}^{ij}(u) = \frac{∂Q^{ν}}{∂u_{ij}^{μ}} Δ_{ν} − (\frac{∂Q^{ν}}{∂u_{ijk}^{μ}} Δ_{ν})_{k} + (\frac{∂Q^{ν}}{∂u_{ijkl}^{μ}} Δ_{ν})_{kl} − …  (29)

where summation is done over double indices.

These formulas are programmed in the more compact form

V = Q^{ν}Δ_{ν},
W^{i} = n(i)u^{μ} \frac{∂V}{∂u_{i}^{μ}} +

n(ij) (u_{j}^{μ} − u^{μ}D_{j}) \frac{∂V}{∂u_{ij}^{μ}} +

n(ijk) (u_{jk}^{μ} − u_{j}^{μ}D_{k} + u^{μ}D_{j}D_{k}) \frac{∂V}{∂u_{ijk}^{μ}} +

⋮

T^{i} = x^{i} ∫_{0}^{1} dλ λ^{p−1}V |_{u→0,x→λx}
P^{i} = T^{i} + ∫_{0}^{1} \frac{dλ}{λ}W^{i}|_{u→λu}  (30)

where in W^{i} it is summed over equal indices (not the i, j, k, … in n(i, j, …)) with i ≤ j ≤ k ≤ … and n(i, j, …) = ∏_{a} r_{a}!/(∑_{b} r_{b})! with r_{a} being the multiplicities of different arguments i, j, … of n (e.g. n(i) = 1, n(i, i) = 1, n(1, 2) = 1/2) which also occur in (27) - (29). p is the number of variables x^{i} and T^{i} are non-zero only if u ≡ 0 does not solve Δ = 0 (T^{i} have to enter (26) in that case as well).

Despite being an elegant formula there may be problems in computing the integral analytically. More seriously, the integral may be singular for λ = 0, 1. That is the case, for example, for the non-polynomial characteristic functions of the Harry-Dym equations in the next section.

If taking the limit u → 0 in V |_{u→0,x→λx} or taking the limit u → λu in W^{i}|_{u→λu} produces a singularity then this can be circumvented by choosing a ũ(x) (e.g. ũ = 1 or ũ = x) such that after replacing u → λu + (1 − λ)ũ the limits V |_{u→ũ,x→λx} and W^{i}|_{u→λu+(1−λ)ũ} become non-singular [2].

Because of these potential difficulties, the default procedure to compute P^{i} is to use the integration module of CRACK to x^{1}-integrate ∑_{ν} Q^{ν}Δ_{ν}, to x^{2}- integrate the remaining unintegrated terms, and so on. In case terms remain after the last x^{p}-integration, the process is restarted on the remaining terms until all terms are integrated or at most a fixed number of times. If this method does not work because not all determining conditions can be solved as, for example, for the Burgers equation below then (30) is used.

If both methods fail then the program concludes that solutions Q of (4) represent adjoint symmetries and not characteristic functions of conservation laws.

4.5 The simplification of P in two variables

After deleting trivial conservation laws and identifying equivalent conservation laws through the computation of characteristic functions Q it remains to simplify the conserved current P through the addition of some curl: P → P + curl V. This is done if there are only two independent variables, say x^{1}, x^{2}. The aim is to lower the order of x^{2}-derivatives in P^{1} through changes P^{1} → P^{1} − D_{2}R, P^{2} → P^{2} + D_{1}R. R is found by repeated partial integration of terms in P^{1} with highest x^{2}-derivatives of u. For that, partial integration routines of CRACK are used which are limited in applicability to expressions which are at most polynomially non-linear in u and derivatives of u.

5 Examples

Computation times refer to a 24 MB REDUCE 3.6 session under LINUX on a 133 MHz Pentium PC with the Jan. 1998 version of CRACK.

Example 1: The advantage of using the package CRACK for solving determining equations is that they can be PDEs and do not have to be restricted to algebraic equations for coefficients of a polynomial ansatz for the conservation law. By that it is possible to find non-polynomial conservation laws and conservation laws that have an explicit x^{i} dependence. An example is the Harry Dym equation

Δ = u_{t} − u^{3}u_{xxx}, u = u(t, x)

which was used below to substitute u_{t} and derivatives of u_{t}. These calculations were done with CONLAW1. In the following the calculation to find conservation laws with P^{t} of order 0 is described. The results for this case and the higher order cases are given afterwards.

If P^{t} is of order 0, i.e. P^{t} = P^{t}(x, t, u), the condition (1): Div P = 0 will involve u_{t} from D_{t}P^{t}. After substitution of u_{t} from the Harry Dym equation it will have u_{xxx} as the highest u-derivative which has to match the highest u-derivative in D_{x}P^{x}. Therefore, CONLAW1 makes the ansatz P^{t} = P^{t}(x, t, u), P^{x} = P^{x}(x, t, u, u_{x}, u_{xx}) for conservation laws with P^{t} of order 0, then formulates the condition Div P = 0 as

u_{xxx}(u^{3}P_{u}^{t} + P_{u_{xx}}^{x}) + P_{t}^{t} + P_{x}^{x} + P_{u}^{x}u_{x} + P_{u_{x}}^{x}u_{xx} = 0  (31)

and calls CRACK to solve it. As P^{t}, P^{x} do not depend on u_{xxx}, a separation performed by CRACK gives

0 = u^{3}P_{u}^{t} + P_{u_{xx}}^{x}  (32)
0 = P_{t}^{t} + P_{x}^{x} + P_{u}^{x}u_{x} + P_{u_{x}}^{x}u_{xx}.  (33)

P^{t} does not depend on u_{xx}, so integration of (32) with respect to u_{xx} introduces a function of integration g_{1}(x, t, u, u_{x}) and substitution of P^{x} in (33) gives an equation that can be separated with respect to u_{xx}. Two equations for P^{t}(x, t, u) and g_{1}(x, t, u, u_{x}) result. One of the equations can be integrated with respect to u_{x} with a function of integration g_{2}(x, t, u) such that g_{1} can be substituted in the other equation which then becomes separable with respect to u_{x}, giving 4 equations. One of the four equations 3P_{xu}^{t} + uP_{xuu}^{t} = 0 can be integrated 3 times to get u^{2}P^{t} + g_{3}(t, u) + g_{4}(x, t) + u^{2}g_{5}(x, t)/2 = 0. Substituting P^{t} in the other 3 equations yields among 2 other equations g_{3},_{uuu} = 0 which integrated and g_{3} being substituted provides 2 equations. One of these equations can be u-integrated such that g_{2} can be substituted in the other equation, which in turn becomes u-separable, giving 4 equations.

The computation so far needed 14 steps in CRACK. Continuing this process of integrations, substitutions, separations for another 14 steps gives the solution

P^{t} = −(g_{17},_{x} u^{2} + 2ug_{10} + 2xg_{14} + x^{2}g_{15} + 2g_{16} − g_{7}u^{2})/(2u^{2}) (34)
P^{x} = (g_{17},_{t} −xg_{7},_{t} −2uu_{xx}g_{10} + u_{x}^{2}g_{10} + 4u_{x}g_{14} − 4xu_{xx}g_{14}
−4ug_{15} + 4xu_{x}g_{15} − 2x^{2}u_{xx}g_{15} − 4u_{xx}g_{16})/2  (35)

with free g_{17}(x, t), g_{7}(t), g_{10}, g_{14}, g_{15}, g_{16}. The algorithm of section 4.1 finds g_{7} to be redundant and sets it to zero which concludes the run of CRACK. In CONLAW1, the terms in the general solution returned by CRACK which involve g_{7} are found to fulfill Div P = 0 identically. Therefore this conservation law is trivial and is dropped. The coefficients of the remaining 4 constants give 4 non-trivial conservation laws. As described in section 4.3 the characteristic functions Q are computed which finally gives: P^{t} of order 0: time to formulate (1): 0.32 sec, to solve (1): 1.34 sec, 28 steps:

\begin{aligned}2u^{−2} · Δ & = & D_{t}(−2u^{−1}) & + & D_{x}(u_{x}^{2} − 2uu_{xx}) \\ 2u^{−3} · Δ & = & D_{t}(−u^{−2}) & + & D_{x}(−2u_{xx}) \\ 2xu^{−3} · Δ & = & D_{t}(−xu^{−2}) & + & D_{x}(2u_{x} − 2xu_{xx}) \\ 2x^{2}u^{−3} · Δ & = & D_{t}(−x^{2}u^{−2}) & + & D_{x}(4xu_{x} − 2x^{2}u_{xx} − 4u) \end{aligned}

Similar computations for higher orders give: P^{t} of order 1: time to formulate (1): 0.32 sec, to solve (1): 2.6 sec, 49 steps:

(2uu_{xx} − u_{x}^{2})u^{−2} · Δ =
D_{t}(−u_{x}^{2}u^{−1}) + D_{x}((2u_{t}u_{x} − u_{xx}^{2}u^{3} + u_{xx}u_{x}^{2}u^{2} − u_{x}^{4}u/4)u^{−1})

P^{t} of order 2: time to formulate (1): 1.7 sec, to solve (1): 110 sec, 193 steps^{3}:

(−8u_{xxxx}u^{3} − 16u_{xxx}u_{x}u^{2} − 12u_{xx}^{2}u^{2} + 12u_{xx}u_{x}^{2}u − 3u_{x}^{4})u^{−2} · Δ =

3The times and the number of steps depend strongly on the relative priorities of the individual steps done (integrations, substitutions,...) which is especially the case for longer computations and therefore these measures characterize the complexity only roughly.

D_{t}((−4u_{xx}^{2}u^{2} − 3u_{xx}u_{x}^{5}tu − u_{x}^{4})u^{−1})+
D_{x}((8u_{tx}u_{xx}u^{2} + 3u_{tx}u_{x}^{5}tu − 8u_{t}u_{xxx}u^{2} − 8u_{t}u_{xx}u_{x}u + 4u_{t}u_{x}^{3}+
4u_{xxx}^{2}u^{5} + 4u_{xx}^{3}u^{4} − 6u_{xx}^{2}u_{x}^{2}u^{3} + 3u_{xx}u_{x}^{4}u^{2})u^{−1})

Computations for conservation laws of higher order of equations Δ = 0 of higher order in more variables x^{i} and more functions u^{α} are similar to the one described above. However, the application of integrability conditions may be necessary (like in a computation for a differential Gröbner basis) and a more general separation method may be necessary if no function depends on all variables and each variable occurs in at least one function. More difficult computations may involve several thousand steps.

Example 2: The Burgers equation in the form

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

is an example for the case that the determining equations cannot be solved completely. It has 0th order conservation laws

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

with f = f(t, x) satisfying the linear reverse heat equation f_{t} + f_{xx} = 0. The above law is also an example that CONLAW allows the computation of conservation laws with non-rational terms which is not possible with approaches based on a polynomial ansatz. A remaining linear PDE and the occurrence of free functions in the conservation law indicates linearizability of Δ = 0, which is the case with the Burgers equation.

Example 3: The MVDNLS equations (Modified Vector Derivative Nonlinear Schrödinger equations) describe oblique propagation of magneto-hydrodynamic waves in warm plasmas [8]. For functions u = u(t, x), v = v(t, x) and b = const. they are

Δ_{1} = u_{t} + [u(u^{2} + v^{2}) + bu − v_{x}]_{x}  (38)
Δ_{2} = v_{t} + [v(u^{2} + v^{2}) + u_{x}]_{x}.  (39)

Both equations have the form of conservation laws. Using the abbreviations (introduced by hand afterwards)

E = −v_{x} + u(u^{2} + v^{2})
F = u_{x} + v(u^{2} + v^{2} − b)
G = 2u_{xx} + 6v_{x}(u^{2} + v^{2}) − 3u(u^{2} + v^{2})^{2} − 2bu^{3}
H = 2v_{xx} − 6u_{x}(u^{2} + v^{2}) − 3v(u^{2} + v^{2})^{2} + 2bv^{3}
I = b(u^{4} − v^{4}) + (u^{2} + v^{2})^{3} − 2u_{x}^{2} − 2v_{x}^{2}

and using equations (38), (39) to substitute for u_{t}, v_{t}, further conservation laws calculated by CONLAW2/3 have the characteristics {Q^{1}, Q^{2}} :

{u, v}, {E, F }, {G, H},  (40)
{(bt − 2x)E − 2tG + b(bt − x)u + v, (bt − 2x)F − 2tH + b(bt − x)v − u}, (41)
{−H_{x} + 2uvH + (b + 2u^{2})G + uI, G_{x} + 2uvG + 2v^{2}H + vI}.  (42)

CONLAW2 can compute one more conservation law with Q^{1}, Q^{2} of 4th order and 36 terms each. Run times are listed in table 1.^{4}

Apart from (41) these conservation laws are given in [8] where also a bi-Hamiltonian structure is provided. Although an infinite sequence of conserved densities can be calculated from the resulting recursion operator, the conservation law (41) is not contained in that sequence and is new - it has an explicit t, x-dependence.

In the scope of CONLAW1 to find conservation laws with P^{1} of order 1 are conservation laws (40),(41).

6 Comparison of the three methods

The determining equations (1)-(4) differ in the number of functions, number of variables and their order. For example, for the MVDNLS equations (38),(39) the condition (1) for conservation laws with P^{1} of order 2 and the conditions (2),(4) for conservation laws with Q^{μ} of order 3 have the following characteristics:

4With the improved (Sep. 1998) version of CRACK times are reduced. For example, the 1977 sec are now 76 sec and the CONLAW1 second order case can now be solved in 880 sec.

order of P^{t} for CONLAW1, order of Q for CONLAW2/3

  CONLAW  0  1  2  3  4

t_{1} t_{2}  t_{1}  t_{2} t_{1}  t_{2} t_{1} t_{2} t_{1}  t_{2}
1  0.15 2.9 0.15 1977
2 1.7 2.0 2.7 16 4.5 194 8.5 722 17 2784
3  0.17 4.5 0.18 11.7 0.3 28.5 0.6 377 1.9 low memory

Table 1: Run times t_{1} to formulate and t_{2} to solve determining conditions of conservation laws of the MVDNLS equations

(1): 1 condition in 12 variables (t, x, u, v, u_{x}, v_{x},…, u_{4x}, v_{4x}), 2 of which occur only explicitly (u_{4x}, v_{4x}), with 55 terms linear in functions P^{t} of 8 variables (t, x, u, v, u_{x}, v_{x}, u_{xx}, v_{xx}) and P^{x} of 10 variables (t, x, u, v,…, u_{xxx}, v_{xxx}) and their 1st order derivatives. The difference in the dependencies of P^{t}, P^{x} at the beginning of CONLAW1 is a consequence of using (38), (39) for substituting a first order t-derivative of u by a second order x-derivative.

(2): 1 condition in 22 variables (t, x, u, v, …, u^{(3)}, v^{(3)}), 6 of which occur only explicitly (2nd order derivatives of u_{t}, v_{t}), with 37 terms linear in functions P^{t}, P^{x} of 14 variables (t, x, u, v, …, u^{(2)}, v^{(2)}) and their 1st order derivatives, and furthermore functions Q^{1}, Q^{2} of 10 variables (t, x, u, v, …, u_{xxx}, v_{xxx}).

(4): 2 coupled conditions in 14 variables (t, x, u, v, u_{x}, v_{x}, …, u_{5x}, v_{5x}), 4 of which occur only explicitly (u_{4x}, v_{4x}, u_{5x}, v_{5x}), with 131 and 132 terms linear in functions Q^{1}, Q^{2} of 10 variables (t, x, u, v, …, u_{xxx}, v_{xxx}) and their 1st and 2nd order derivatives.

The following are general features of equations (1)-(4). Equation (1) is of first order and therefore only highest order u-derivatives which are not substituted due to Δ = 0 are not variables of the P^{i} and can be used for direct separation. Equation (1) therefore is only weakly over determined with the application of integrability conditions playing an important role. A general problem with computing a differential Gröbner Basis is that the complexity of these calculations depends heavily on the total ordering of derivatives of functions P, Q chosen for which there is currently no complete theory available. Choices made by the program can be particularly good or bad for the problem at hand.

In contrast, equations (4) are of higher order with more jet-variables that occur only explicitly and that can be used for direct separation. Although these equations are of higher order they are highly over determined and simpler to solve in general. An efficient way of doing direct separations and handling large equations is of importance for this approach.

Finally, in equations (2) the P^{i} depend initially on all jet-variables (apart from highest order u-derivatives), and also those substituted through Δ = 0 on which the Q^{μ} do not depend. On the other hand, the Q^{μ} do depend on highest order u-derivatives initially. The efficiency in solving (2) therefore depends on the efficiency of a module for indirect separation, i.e. on a module for handling equations which have no function depending on all variables but which have also no variable occurring only explicitly so that no direct separation with respect to any variable is possible. Such a module is described in [9].

All techniques are used to solve the over determined system of all three approaches, only some are used more often in one approach than in the other.

There is another issue. If the order of derivatives with respect to different variables differs, such as for example the Harry Dym equation u_{t}−u^{3}u_{xxx} = 0, then it matters whether this equation is used to do substitutions u_{t} = u^{3}u_{xxx} or u_{xxx} = u_{t}/u^{3}. Substituting u_{t} gives a lower increase in complexity when successively higher order ansätze for P or Q are made. On the other hand, one has to go to higher orders of P and Q to cover the same equivalence classes of conservation laws compared to substituting u_{xxx}. As equations (4) involve already higher order u-derivatives, a further increase could explode the size of (4) even more.

Another relation between (2) and (4) is that one could look at (4) as resulting from a differential-Gröbner-Basis calculation done with (2), with the aim to eliminate the P^{i} first. It is of course more efficient to exploit knowledge of the structure of (2) and to apply the Euler operator to write down (4) directly rather than to do the differential Gröbner Basis calculation step by step with (2). On the other hand, CRACK includes a number of modules to take advantage of special situations (e.g. to integrate exact PDEs or to recognize and solve PDEs that are ODEs for some partial derivatives and to solve them using ODESOLVE [6]). For a concrete problem it is very likely that there exists a quicker way to solve (2) than to first eliminate all P^{i}. The question of which of the CONLAW programs is best depends on the effectiveness of different submodules of the program CRACK which is used to solve conditions (1)-(4). With the current version of CRACK (Jan. 1998), programs CONLAW1/3 are better for simpler conservation law problems while CONLAW2 is better for larger problems.

7 Syntax of CONLAW

Example: The input to find conservation laws with Q of order 0-4 for the MVDNLS equations (38),(39) is

depend u,x,t;
depend v,x,t;
conlaw2({{df(u,t) = - df( u*(u**2+v**2) + b*u - df(v,x) ,x),
          df(v,t) = - df( v*(u**2+v**2) + df(u,x)       ,x) },
         {u,v}, {t,x}
        },
        {0, 4, t, {}, {}});

In REDUCE lists are enclosed in { }. The input of CONLAWi (i=1,2,3) consists of two lists, of which the first encodes the PDE problem. It contains a list of equations with the derivative to be substituted on the left hand side, a list of functions and a list of independent variables. The second parameter to CONLAWi is a list that specifies the conservation laws to be computed. Its first two elements are the minimum and maximum order of P^{1} in the case of CONLAW1 and the order of Q^{μ} in the case of CONLAW2/3. The third element is t or nil and specifies whether the conservation law may depend explicitly on the x^{i} or not. The fourth element is a list of functions to be determined in an ansatz made for P^{i} or Q^{μ} and the last element is a list of inequalities to be satisfied.

More details about investigating an ansatz is given in a manual file that comes with the three CONLAWi.RED files.

8 Summary

Supplied with subroutines to fix gauge freedom in differential expressions, the programs CONLAW1/2/3 proved to be an efficient tool for the computation of conservation laws of differential equations. Compared with other programs, (see for instance [3]) the programs CONLAWi show the following new features:

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

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

• There is no limit to 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.

• It is possible to determine values of parameters in the DE such that conservation laws exist (examples in [11]).

• For each of the programs CONLAWi, an ansatz for P^{i} and/or Q^{μ} can be input to specify conservation laws to be calculated.

Compared with the program of Göktaş and Hereman, CONLAW is able to find more general conservation laws and to make a definitive statement if local conservation laws do not exist and if the order is not too high to complete the computations.

The strength of the program described in [3] is to sometimes get a higher order that still can be handled by concentrating on polynomial conservation laws because the program only has to solve algebraic systems for coefficients of a polynomial ansatz. They were also able to extend applicability to differential-difference systems [4], [5].

The comparison of the three approaches (1)-(4) showed that each of them has advantages in special circumstances. It also serves as a comparison between using a general purpose program to find the quickest way of solving over determined PDE systems directly (CONLAW1/3) and an approach to derive integrability conditions by applying extra information about the structure of the PDE system (CONLAW2).

The programs including a manual and a test file are available via ftp

from ftp://ftp.maths.qmw.ac.uk/pub/tw. The package will be submitted

to the REDUCE network library.

9 Acknowledgement

The authors want to thank Malcolm MacCallum, George Labahn and the referees for comments on an earlier version of the manuscript. Also, Jan Sanders and Willy Hereman are thanked for discussions along a research visit of TW at CAN Netherlands/Amsterdam. Discussions with Winfried Neun during multiple visits at the Konrad Zuse Institut/Berlin are gratefully acknowledged. Finally, Frank Verheest, Willy Sarlet, Micheline Musette (Belgium) and the Relativity group at Hall University (England) provided PDEs to test the code.

References

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

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

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

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

[5] Göktaş, Ü. and Hereman, W. (1998) Computation of conserved densities for nonlinear lattices, Physica D (1998) in press.

[6] MacCallum, MAH. (1988) An Ordinary Differential Equation Solver for REDUCE Proc. ISSAC‘88, Springer Lect. Notes in Comp. Sc., 358, 196-205.

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

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

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

[10] 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.

[11] Wolf, T. A comparison of four approaches to the calculation of conservation laws, in preparation.