reduce-historical/r37/packages/crack/conca.tex

% This is a LaTeX file
\documentstyle[12pt]{article}
%Sets size of page and margins
%\oddsidemargin -2mm  \evensidemargin -2mm
\topmargin -25mm %   \headheight 0pt   \headsep 0pt
%\footheight 14pt  \footskip 40pt
\textheight 25.5cm  
%\textwidth 15.4cm

%spaces lines at one and a half spacing
%\def\baselinestretch{1.5}

%\parskip = \baselineskip

\title{\mbox{\quad} \\ \mbox{\quad} \\  \mbox{\quad} \\  \mbox{\quad} \\ 
       Computer algebra algorithms and routines for the
       computation of conservation laws and fixing of gauge
       in differential expressions}

\author{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\\ Institute for Informatik, Friedrich Schiller
     Universit\"{a}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
     }

\begin{document}
\maketitle
\begin{abstract}
Three different approaches for the determination of conservation
laws of differential equations are presented.
For three corresponding {\tt REDUCE} computer algebra programs 
{\tt CONLAW1/2/3} the necessary subroutines
are described. One of them simplifies general solutions of overdetermined
PDE systems so that all remaining free functions and constants
correspond to independent conservation laws.
It determines redundant functions and constants in differential expressions
and is equally useful for the determination of symmetries or the
fixing of gauge freedom in differential expressions.
\end{abstract}
%\tableofcontents
%\listoftables
%\listoffigures

\section{Introduction}
The determination of conservation laws (CLs) 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, classification and for supporting 
their numerical solution. In this paper we outline three computer
algebra programs for the computation of CLs and explain
in more detail subroutines to fix gauge freedom in differential
expressions which in this context is used to extract individual CLs
from the general solution of CL-determining equations.

In the following we adopt the notation of the book of Olver \cite{PO}.
Independent variables will be denoted by $x = (x^1, x^2, \ldots , x^p)$.
The differential equations are $\Delta(x,u^{(n)}) = 0
\;\;\;(\mbox{i.e.}\;\, \Delta_1=0, \ldots , \Delta_q=0)$,
for $q$ functions $u = (u^1, u^2, \ldots , u^q),\;\; u^{(n)}$ denoting
$u$-derivatives of order up to $n.$ The conservation law that is to be
fulfilled by solutions of $\Delta = 0$ is $\mbox{Div}\,P = 0$ with conserved
current $P = (P^1, \ldots , P^p).$ We will use $_J$ as a multiple index
denoting partial derivatives, for example, $u^{\alpha}_{J}$ will stand for
an arbitrary partial derivative, like
$\partial^k u^{\alpha}/(\partial x^1\partial x^2\ldots)$.
%, which also may be denoted as $D_J u^{\alpha}$.

If the differential equations $\Delta = 0$ result from a variational
principle then any Lie-symmetry of $\Delta = 0$ provides a conservation
law as is known from Noether's Theorem. Instead, we will not make
any restrictive assumptions which leaves us to solve $\mbox{Div}\,P = 0$
either directly or to determine characteristic functions of conservation
laws or to do both at once. A comparison of these different approaches
with respect to their complexity, and an extension to find non-local
conservation laws and applications to PDEs with parameters will be
described elsewhere \cite{TW}; here we concentrate on the computer 
algebra aspects.

\section{The mathematical problem and the three approaches}
In this section we describe three ways to formulate
determining conditions for conservation laws.

The first and most direct approach is to solve 
\begin{equation} 
\mbox{Div}\,P = 0                                     \label{a1} 
\end{equation}
modulo $\Delta = 0$ directly. The corresponding program is
{\tt CONLAW1}.
The components of the conserved current $P^1,\ldots,P^p$ that are 
to be calculated
are functions of all independent variables $x^i$, the dependent variables
$u^{\alpha}$ and their derivatives $u^{\alpha}_{J}$ up to some order.

Because we are not interested in trivial CLs $P = {\tt curl}\; V$
but in CLs that solutions of $\Delta = 0$ obey,
we use $\Delta = 0$ to eliminate some 
of the so-called jet-variables
$u^{\alpha}_{J}$ and substitute them
in the determining conditions (\ref{a1}). 
By that, the conditions (\ref{a1}) have to be fulfilled identically
in less variables, they become 
less restrictive and they may have additional solutions apart from 
$P = {\tt curl}\; V$. These extra non-trivial CLs are the
ones of interest.
We therefore assume $\Delta = 0$ can be solved for 
leading derivatives $u^{\alpha}_{J}$ so that they and 
all their partial derivatives that occur in (\ref{a1}) can be substituted.
We also, w.l.o.g., assume that the $P^i$ do not depend on 
$u$-derivatives we substitute, which fixes a kind of
equivalence of CLs.

Other approaches calculate characteristic functions $Q^{\nu}$.
A theorem can be proven (\cite{PO}, p.\ 272) that for a totally
non-degenerate system $\Delta_{\nu}=0$, 
each equivalence class of CLs
$\mbox{Div}\,P = 0$ (i.e.\ conserved currents differing only by a curl)
is determined uniquely by characteristic functions
$Q^{\nu}$ satisfying
\begin{equation}
\mbox{Div}\,P = \sum_{\nu} Q^{\nu} \Delta_{\nu}  \label{a2}
\end{equation}
identically in {\it all} $x^i,u^{\alpha},u^{\alpha}_{J}$.
Equ.\ (\ref{a2}) is not solved by simply eliminating $Q^1$ in terms
of $P$ and $\Delta$ and other $Q^{\nu}$ as it would be singular for
solutions of $\Delta=0$. To avoid this and because the $Q^{\nu}$ are 
unique only modulo $\Delta=0$,
we w.l.o.g.\ ignore dependencies of $Q^\nu$ on
leading $u$-derivatives in $\Delta=0$ and any of their derivatives.
A way to calculate the $Q^{\nu}$ is to use the property of the
Euler operators $E_{\nu} = \sum_J (-D)_J \partial/\partial(u^{\nu}_{J})$ 
which acting
on an expression gives identically zero iff this expression is a divergence.
The $D$ are total derivatives.
Applying this operator on (\ref{a2}) and using $\Delta_{\nu}=0$ one
obtains as determining conditions for the $Q^{\nu}$:
\begin{equation}
0 = \sum_{\mu,J} (-D)_J
   \left( Q^{\mu} \frac{\partial \Delta_{\mu}}
                       {\partial(u^{\nu}_{J})}
   \right) \;\;\; \forall \nu.                           \label{a3}
\end{equation}
The second and third approach are to solve identically in
$x^i,u_{\alpha},u^{\alpha}_{J}$ either 
(\ref{a3}) for $Q^{\nu}$ or (\ref{a2}) for $P^i, Q^{\nu}$.
The corresponding programs are
{\tt CONLAW2} for (\ref{a3}) and {\tt CONLAW3} for (\ref{a2}).

The three approaches (\ref{a1})-(\ref{a3})
differ in the number of equations to be solved or their order or the
number of functions to be determined or the number of independent
jet-variables or the degree of an ansatz for $P,Q$ in order to
obtain the same conservation law. 

To obtain solutions of (\ref{a1})-(\ref{a3}) we assume
bounds on the order of $u$-derivatives on which the $P^i$ and $Q^{\nu}$
may depend. For (\ref{a1}) we assume a bound for $P^1$ and for
(\ref{a2}),(\ref{a3}) we assume a bound for $Q^{\nu}$.
Bounds for the remaining unknown functions follow.
Differentiations done in all three conditions
(\ref{a1})-(\ref{a3}) introduce jet-variables ($u$-derivatives)
on which the $P^i$ resp.\ $Q^{\nu}$ do not depend so that overdetermined
conditions result in which there is no unknown function $P^i,Q^{\nu}$
of all jet-variables $u^{\alpha}_{J}$
in which the conditions have to be satisfied identically. 
The resulting overdetermined
PDE-systems are investigated with the computer algebra
package {\tt CRACK}.

\section{The computer algebra problem}
The main computer algebra problem is to solve the overdetermined
conditions (\ref{a1})-(\ref{a3}). Steps undertaken  
include the separation, integration, application of integrability
conditions (differential Gr\"{o}bner Basis), solution of ODEs and other steps
which are described in \cite{CRACK1},\cite{CRACK2}.

If the overdetermined system is linear ((\ref{a1})-(\ref{a3}) are
linear in $P^i,Q^{\nu}$)  and not too big - we give an example below
for what is currently possible - then {\tt CRACK} will solve the
system either completely or partially and return unsolved equations 
e.g.\ return the heat equation when investigating conservation laws of the
Burgers equation.

In the general solution of the CL-condition(s) a CL is extracted by
collecting all terms involving one of the arbitrary constants or
arbitrary functions in the solution. If some of them were redundant
then CLs extracted would not be independent of each other.

Redundant constants and functions may result because
in the process of solving the overdetermined system there is no general
rule for what should have a higher priority, integrations or 
the application of integrability conditions,
as 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 x t + 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) \rightarrow 
c_1(x) - c_2 x$.

Recognizing redundancy can become cumbersome in the case of many
independent variables or if arbitrary constants/functions appear 
non-linearly.

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 (\ref{a1}) the general
solution could be augmented by ${\tt curl}\; V$ and $V$ be added
to the list of free constants and functions. 
In this way trivial CLs could be filtered
out as the free constants and functions corresponding to them
would be redundant to $V$.

Although in the case of computing CLs, one easily could drop
trivial CLs after they have been computed by checking
$\mbox{Div}\,P = 0$ identically in all jet-variables,
such a simple test to eliminate gauge might not be available
for other problems.

\section{Subroutines}
In the following subsections we describe subroutines which
extract CLs from the general solution of conditions 
(\ref{a1})-(\ref{a3}), subroutines that compute $Q^{\nu}$ from
$P^i$ and $P^i$ from $Q^{\nu}$ and subroutines that simplify $P^i$.

\subsection{Identifying redundant constants and functions}
The problem of finding the general solution of a PDE system with
some existing gauge freedom fixed can be reduced to the problem 
of finding the general solution of a PDE system without fixing
gauge in the following way.

Given a system of DEs $0 = \Omega(f_a,x^i)$ to be solved for the
functions $f_a(x^i)$, we assume 
\begin{equation}
f_b = F_b(x^i,g_c) \label{c0} 
\end{equation}
to be a general solution where $F_b$ are differential expressions in 
$x^i,g_c$ where $g_c$ are arbitrary constants and functions.
They may include functions from the original set $f_a$  
and constants and functions of integration.

The question is to specify the $g_c$ to fix any redundancy but
not to lose generality of the solution. The steps are:
\begin{itemize}
\item Formulate a set of conditions 
\begin{equation}
 0 = F_b(x^i,g_c) - F_b(x^i,\bar{g}_c)     \label{c1}
\end{equation}
where $\bar{g}_c$ are new functions, each $\bar{g}_c$ having the same variable
dependence as $g_c$. Regard equ.\ (\ref{c1}) as a system
of equations for the set of unknown functions $\{g_c, \bar{g}_c\}$,
to be satisfied identically in the $x^i$.
\item Find the general solution of the system (\ref{c1}) as
\begin{equation}
\tilde{g}_c = G_c(x^i,h_d)  \label{c2}
\end{equation}
where $\tilde{g}_c$ is a subset of $\{g_a,\bar{g}_b\}$, and
$G_c$ are algebraic or differential expressions in functions
$h_d$ which are the remaining $\{g_a,\bar{g}_b\}$ and extra constants and functions
of integration. The $h_d$ are arbitrary. Any function $g_a$ or $\bar{g}_a$
appears only once on a left-hand-side (lhs) of (\ref{c2}) or only on 
right-hand-sides (rhs's).
\item
If for any index $c$ both, $g_c$ and $\bar{g}_c$ appear only on rhs's
of (\ref{c2}) then $g_c$ is redundant and can be set to zero in all $F_b$
in (\ref{c0}) and all $G_c$ in (\ref{c2}).
\item
If for any index $c$ both, $g_c$ and $\bar{g}_c$ appear only on 
lhs's of (\ref{c2}) in the equations $g_c=G_c$ and $\bar{g}_c=\bar{G}_c$
then these two equations are replaced by $\bar{g}_c = g_c-G_c+\bar{G}_c$ in (\ref{c2}).
\item
If for any index $c$, $g_c$ appears on a lhs of (\ref{c2}) 
and $\bar{g}_c$ appears only on rhs's then the equation with lhs $g_c$ is solved
for $\bar{g}_c$ in terms of $g_c$ and other functions and replaced by the
new equation $\bar{g}_c=\bar{G}_c(g_c,\ldots)$. With this new equation $\bar{g}_c$ is
substituted on any rhs of (\ref{c2}).
\item
There remains only the case of $\bar{g}_c$ being on the lhs of an equation
and $g_c$ being on rhs's such that the system (\ref{c2}) now has the form 
\begin{equation}
\bar{g}_c = \bar{G}_c(x^i,g_a,\bar{h}_b)    \label{c3}
\end{equation}
where $\bar{h}_b$ are arbitrary constants and functions of integration which
arose during the solution of (\ref{c1}). $\bar{g}_c$ do not occur on rhs's
as they would be redundant and would have been set to zero otherwise.
\item
Finally, free constants and functions $\bar{h}_b$ on rhs's will be chosen
to make as many $\bar{G}_c$ as possible zero 
and to set the redundant $g_c$ to zero in (\ref{c0}) and (\ref{c3}).
As we do not have to know $\bar{h}_b$ explicitly, it is enough to find
equations in (\ref{c3}) which include an arbitrary function $\bar{h}_b$
of all variables $x^i$ in this equation. Assuming local solvability
of $0 = \bar{G}_c$ for $\bar{h}_b$ we conclude redundancy of $g_c$.
\item 
All remaining $\bar{h}_b$ which cannot be used to make a rhs zero are
set to zero themselves and the final form of (\ref{c3}) 
$\bar{g}_c = \bar{G}_c(x^i,g_a)$ provides substitutions which turn 
$F_b(x^i,\bar{g}_c)$ into the gauge fixed final solution 
$f_b = F_b(x^i,g_c)$.
\end{itemize}
\noindent Two comments: 

Although the solvability of (\ref{c1}) for $g_a,\bar{g}_b$ and the
solvability of $0 = g_c - G_c(x^i,\bar{g}_c,\ldots)$ for $\bar{g}_c$ cannot
be guaranteed, this should in practice not be a problem for the
following reasons.
\begin{itemize}
\item
Usually there is no arbitrary function $g_a,\bar{g}_b$ depending
on all (jet-) variables
of (\ref{c1}) such that (\ref{c1}) is very overdetermined and
therefore easy to solve. 
\item
If the equ.s $0=\Omega$ are linear in $f_a$ then their solution is
linear in the arbitrary functions $g_c$ which is the case for the
computation of CLs \footnote{Conditions become non-linear
if we want to calculate parameter values such that CLs exist.}.
\item
If equ.s $0=\Omega$ are non-linear in $f_a$ then 
solving (\ref{c1}) should still be simpler than the solution of $0=\Omega$
which we assume was possible to derive.
\item
Equ.s (\ref{c1}) have the special 
solution $\bar{g}_c = g_c, \;\;\forall c$.
\end{itemize}


The above steps for fixing gauge freedom are not only 
applicable once a general solution of a PDE(-system) 
$0 = \Omega(f_a,x^i)$ has already been found. For example, the
computation of conservation laws for the Burgers equation below
returns the heat equation which remains unsolved.
In order to find redundancies in constants and functions which 
turn up in a preliminary solution $f_b = F_b(x^i,g_c)$ and 
which additionally have to satisfy remaining differential equations
$0 = D(x^i,g_c)$, one can extend redundancy conditions (\ref{c1})
by $ 0 = D(x^i,g_c) - D(x^i,\bar{g}_c)$. These conditions are sufficient
but not necessary as only equivalence of $0=D(x^i,g_c)$ and 
$0=D(x^i,\bar{g}_c)$ is required, not equality.

The possibility to fix at least some gauge freedom even in the 
presence of yet unsolved equations opens the possibility to
run a gauge-fixing step during the process of solving
overdetermined PDE-systems. By that the number of unknown 
functions could be reduced and the remaining equations be simplified.

\subsection{Computing characteristic functions from \\ conserved currents}
The first approach (\ref{a1}) is attractive compared with (\ref{a2}),(\ref{a3})
as it generates only one PDE to be solved which is of first
order and involves less jet-variables than approach (\ref{a2})
because it is computed modulo $\Delta = 0$. Also, it has less functions
to compute than approach (\ref{a2}). A negative aspect is that
it provides only the conserved current $P$ and not the characteristic
functions $Q$. 

If expressions $Q^{\nu J}$ in a relation (\ref{trafo1}) below are known
then partial integrations (\ref{trafo2})
yield the characteristic functions $Q^{\nu}$
and the corresponding conserved current $P-R$:
\begin{eqnarray}
\mbox{Div}\,P = 0 & &\!\!\!\!\!\!\!\mbox{mod}\;\;\Delta_{\nu}=0 \;\;
\leftrightarrow \nonumber \\
\exists\, Q^{\nu J}: \mbox{Div}\,P & = & \sum_{\nu,J} Q^{\nu J}
\Delta_{\nu}^{(J)} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\,
(\mbox{identically in {\it all}} \;\; x, u^{\alpha}_J ) \label{trafo1} \\
& = & \sum_{\nu,J} D_J(Q^{\nu J} \Delta_{\nu}) - D_J(Q^{\nu J}) \Delta_{\nu}
\;\;\;\;\;\;\;\;\;\;\;\;(\mbox{repeatedly}) \label{trafo2}   \\
& = & \mbox{Div}\,R + \sum_{\nu} Q^{\nu} \Delta_{\nu}  \nonumber 
\end{eqnarray}
Equation (\ref{trafo1}) cannot be regarded as a linear algebraic
equation to determine $Q^{\nu J}$ as
there is the additional requirement that the $Q^{\nu J}$ are non-singular
for solutions of $\Delta=0$. Instead,
$\mbox{Div}\,P$ is calculated 
and substitutions of a different form than before
are made. For example, if CLs for the Harry Dym equation 
$0 = \Delta = u_{t} - u^3u_{xxx}$ are investigated and if for the
derivation of (\ref{a1}) there had been done substitutions
$u_{t}=u^3u_{xxx},\;\; u_{tx}=(u^3u_{xxx})_{x},\ldots$ before
then now the substitutions would be
$u_{t}=\Delta+u^3u_{xxx},\,$ 
$u_{tx}=\Delta_{x}+(u^3u_{xxx})_{x},\ldots$
which provide the rhs of (\ref{trafo1}).
The computation of $Q^{\nu}$ and $P^i-R^i$ from $P^i$ is part of {\tt CONLAW1}.

\subsection{Computing conserved currents from \\ 
            characteristic functions}
The inverse computation is necessary in {\tt CONLAW2} where the
conserved current $P^i$ has to be computed from $Q^{\mu}$ by
integrating $\mbox{Div}\,P = \sum_{\nu} Q^{\nu} \Delta_{\nu}$.

A direct way is based on a formula % to compute $P^i$ from $Q^{\mu}$
given by Anco \& Bluman in \cite{AB}:
\begin{eqnarray}
P^i & = & \int^1_0 \frac{d \lambda}{\lambda} 
       \left(S^i(u) + N^i_{\mu}(u)u^{\mu} + N^{ij}_{\mu}(u)D_ju^{\mu} 
        + \ldots\right)|_{u\rightarrow \lambda u}  \label{ab0} \\
S^i(u) & = & Q^{\nu}\frac{\partial \Delta_{\nu}}{\partial u^{\mu}_i}u^{\mu} + 
          Q^{\nu}\frac{\partial \Delta_{\nu}}{\partial u^{\mu}_{ij}}u^{\mu}_j -
          \left(Q^{\nu}\frac{\partial \Delta_{\nu}}{\partial 
                                  u^{\mu}_{ij}}\right)_{j}u^{\mu} + \ldots 
         \label{ab1} \\
N^i_{\mu}(u) & = & \frac{\partial Q^{\nu}}{\partial u^{\mu}_i}\Delta_{\nu} - 
          \left(\frac{\partial Q^{\nu}}{\partial u^{\mu}_{ij}}\Delta_{\nu}\right)_{j} +
          \left(\frac{\partial Q^{\nu}}{\partial u^{\mu}_{ijk}}\Delta_{\nu}\right)_{jk} - \ldots   \label{ab2} \\
N^{ij}_{\mu}(u) & = & \frac{\partial Q^{\nu}}{\partial u^{\mu}_{ij}}\Delta_{\nu} - 
          \left(\frac{\partial Q^{\nu}}{\partial u^{\mu}_{ijk}}\Delta_{\nu}\right)_{k} +
          \left(\frac{\partial Q^{\nu}}{\partial u^{\mu}_{ijkl}}\Delta_{\nu}\right)_{kl} - \ldots  \label{ab3}
\end{eqnarray}
where summation is done over double indices.
% in each term including a not shown combinatorical factor .

A slightly more compact formulation (and way to compute $P^i$) is
\begin{eqnarray} 
V&=&Q^{\nu} \Delta_{\nu}, \nonumber \\
W^i
&=& n(i) u^{\mu} \frac{\partial V}{\partial u^{\mu}_i} +  \nonumber \\
& & n(ij) \left(u^{\mu}_j -u^{\mu}D_j\right) \frac{\partial V}{\partial u^{\mu}_{ij}} +  \nonumber \\
& & n(ijk) \left( u^{\mu}_{jk} - u^{\mu}_jD_k + u^{\mu}D_jD_k \right)
    \frac{\partial V}{\partial u^{\mu}_{ijk}} +  \nonumber \\     
& & \vdots  \nonumber \\
   T^i&=& x^i \int^1_0 d \lambda \lambda^{p-1} V|_{u\rightarrow 0, 
x\rightarrow \lambda x}  \nonumber \\
   P^i&=& T^i + \int^1_0 \frac{d \lambda}{\lambda} W^i|_{u\rightarrow \lambda u} 
\label{ab}
\end{eqnarray}
where in $W^i$ it is summed over equal indices (not the $i,j,k,\ldots$ in $n$) and 
$n(i,j,\ldots)=\prod_a r_a!/ (\sum_b r_b)! $ with $r_a$ being the multiplicities 
of different
arguments $i,j,\ldots$ of $n$ (e.g.\ $n(i)=1, n(i,i)=1, n(1,2)=1/2$)
which also occur in (\ref{ab1}) - (\ref{ab3}). 
$p$ is the number of variables $x^i$ and $T^i$ are non-zero only if
$u\equiv 0$ does not solve $\Delta=0$ ($T^i$ have to enter (\ref{ab0})
in that case as well).

Although being an elegant formula there may be problems in computing the
integral analytically. More seriously, the integral may be singular for
$\lambda=0,1$. That is the case, for example, for the non-polynomial
characteristic functions of the Harry-Dym equations in the next section.
Although in some cases it might help do take $P^i = \int^1 
\frac{d \lambda}{\lambda} W^i|_{u\rightarrow \lambda u}$ 
this need not always be the case.

Because of these potential difficulties
the default procedure to compute $P^i$ is to use
the integration module of {\tt CRACK} to
$x^1$-integrate $\sum_{\nu} Q^{\nu} \Delta_{\nu}$, 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 had been
solved as, for example, for the Burgers equation
below then (\ref{ab}) is used.

\subsection{The simplification of $P$ in two variables}
After deleting trivial CLs and identifying equivalent CLs
through the computation of characteristic functions $Q$ it remains
to simplify the conserved current $P$ through the addition of
some curl: $P \rightarrow P+\mbox{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 \rightarrow P^1 - D_2 R,\,
P^2 \rightarrow 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 {\tt CRACK} are
used which are limited in applicability
to expressions at most polynomially non-linear
in $u$ and derivatives of $u$.

\section{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 {\tt CRACK}.

{\it Example 1}: \newline 
The advantage of using the package {\tt 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 CL. By that it is possible to find non-polynomial
CLs and CLs that have an explicit $x^i$ dependence.
An example is the Harry Dym equation
\[ \Delta = u_{t} - u^3u_{xxx}, \;\;\;\; u = u(t,x) \]
which was used below to substitute $u_{t}$ and derivatives of 
$u_{t}$. These calculations were done with {\tt CONLAW1}.\newline
$P^t$ of order 0: time to formulate (\ref{a1}): 0.32 sec, to solve (\ref{a1}):
1.34 sec, CLs:
\[ \begin{array}{rclcl}
2u^{-2} \cdot \Delta & = & D_t(-2u^{-1}) & + & D_x(u_{x}^{\;\,2}-2uu_{xx}) \vspace{1.5mm} \\
2u^{-3} \cdot \Delta & = & D_t(-u^{-2}) & + & D_x(-2u_{xx}) \vspace{1.5mm} \\
2xu^{-3} \cdot \Delta & = & D_t(-xu^{-2}) & + & D_x(2u_{x}-2xu_{xx}) \vspace{1.5mm} \\
2x^2u^{-3} \cdot \Delta & = & D_t(-x^2u^{-2}) & + & D_x(4xu_{x}-2x^2u_{xx}-4u)
\end{array} \]
$P^t$ of order 1: time to formulate (\ref{a1}): 0.32 sec, to solve (\ref{a1}):
2.6 sec, CLs: 
\[(2uu_{xx}-u_{x}^{\;\,2})u^{-2} \cdot \Delta = \]
\[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 (\ref{a1}): 0.7 sec, to solve (\ref{a1}):
158 sec, CL: 
\[
(-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} \cdot \Delta = \]
\[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})
\]

{\it Example 2}: \newline 
The Burgers equation in the form
\begin{equation}
 \Delta = u_t - u_{xx} - \frac{1}{2}u_x^{\,2} = 0, \;\;\;\; u = u(t,x)
 \label{BE1}
\end{equation} 
is an example for the case that the determining equations cannot be
solved completely. It has zeroth order CLs
\begin{equation}
 fe^{u/2}\Delta = D_t(2fe^{u/2}) + D_x(e^{u/2}(2f_x-fu_x)) \label{BE1cl}
\end{equation}
with $f = f(t,x)$ satisfying the linear reverse heat equation
$0 = f_t + f_{xx}.$ 
This CL is also an example that {\tt CONLAW} allows the computation
of CLs 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 CL indicates linearizability of $\Delta=0$ which
is the case with the Burgers equation. 

{\it Example 3}: \newline
The MVDNLS equations (Modified Vector Derivative Nonlinear Schr\"{o}dinger
equations) describe oblique propagation of magnetohydrodynamic waves in
warm plasmas \cite{NS}. For
functions $u = u(t,x), \; v = v(t,x)$ and $b = $ const.\ they are
\begin{eqnarray}
\Delta_1 & = & u_t + [u(u^2+v^2) + bu - v_x]_x \label{mvdnls1} \\
\Delta_2 & = & v_t + [v(u^2+v^2) + u_x]_x.      \label{mvdnls2}
\end{eqnarray}
Both equations have the form of CLs. Using the abbreviations (introduced
by hand afterwards)
\begin{eqnarray*}
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
\end{eqnarray*}
and using equ.s (\ref{mvdnls1}), (\ref{mvdnls2}) to substitute
for $u_t, v_t$,
further CLs calculated by {\tt CONLAW2/3} have the characteristics
$\{Q^1,Q^2\}\;\,$:
\begin{equation} \{u, v\},\;\; \{E, F\},\;\; \{G, H\}, \label{mvdnls3} 
\end{equation}
\vspace{-5mm} 
\begin{equation} \{(bt-2x)E-2tG+b(bt-x)u+v,\;\, (bt-2x)F-2tH+b(bt-x)v-u\}, 
\label{mvdnls4} \end{equation}
\vspace{-3mm} 
\begin{equation}
\{-H_x+2uvH+(b+2u^2)G+uI,\;\, G_x+2uvG+2v^2H+vI \}. 
\label{mvdnls5} \end{equation}
{\tt CONLAW2} can compute one more CL with $Q^1, Q^2$ of 4'th order
and 36 terms each. Run times are listed in table 1.

Apart from (\ref{mvdnls4}) these CLs are given in \cite{NS} 
where also a bi-Hamiltonian structure is provided.
Although from the resulting recursion operator, an infinite sequence
of conserved densities can be calculated, the CL (\ref{mvdnls4}) is not
contained in that sequence and is new - it has an explicit
$t,x$-dependence. 

\small
\begin{table}
\begin{center}
\begin{tabular}{|c|r|r|r|r|r|r|r|r|r|r|}  \hline
             & \multicolumn{10}{c|}{order of $P^t$ for {\tt CONLAW1}, 
                 order of $Q$ for {\tt CONLAW2/3}} \\ \cline{2-11}
{\tt CONLAW} & 
\multicolumn{2}{c|}{0} &
\multicolumn{2}{c|}{1} &
\multicolumn{2}{c|}{2} &
\multicolumn{2}{c|}{3} &
\multicolumn{2}{c|}{4} \\  \cline{2-11}
 & $t_1$ & $t_2$ & $t_1$ & $t_2$ & $t_1$ & $t_2$ & $t_1$ & $t_2$ & $t_1$ & $t_2$     \\ \hline
 1 & 0.15 & 2.9 & 0.15 & 1977 & & & & & &     \\ \hline
 2 & 1.7 & 2.0 & 2.7 & 16 & 4.5 & 194 & 8.5 & 722 & 17 & 2784 \\ \hline
 3 & 0.17 & 4.5 & 0.18 & 11.7 & 0.3 & 28.5 & 0.6 & 377 & 1.9 & low memory  \\ \hline
\end{tabular}
\caption{Run times $t_1$ to formulate and $t_2$ to solve 
determining conditions of CLs of the MVDNLS equations}
\end{center}
\end{table}
\normalsize
In the scope of
{\tt CONLAW1} to find CLs with $P^1$ of order 1 are CLs 
(\ref{mvdnls3}),(\ref{mvdnls4}) and if equations 
(\ref{mvdnls1}),(\ref{mvdnls2}) are used to substitute $u_{xx}, v_{xx}$
then also (\ref{mvdnls5}) is included.
Such a run of {\tt CONLAW1} returns a differential Gr\"{o}bner
Basis of 2 equations for one function in 3 variables and 2 equations
for one function in 2 variables,
which could not be solved completely because one of the ODEs is a 
second order ODE that could not be solved automatically.

\section{Comparison of the three methods}
The determining equations (\ref{a1})-(\ref{a3}) differ in the
number of functions, number of variables and their order. 

For example, for the MVDNLS equations (\ref{mvdnls1}),(\ref{mvdnls2})
the condition (\ref{a1}) for CLs with $P^1$ of order 2 and
the conditions (\ref{a2}),(\ref{a3}) for CLs with $Q^{\mu}$ of
order 3 have the following characteristics:

(\ref{a1}): 1 condition in 12 variables 
$(t,x,u,v,u_x,v_x,\!\ldots\!,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,\!\ldots\!,u_{xxx},v_{xxx})$
and their 1st order derivatives. The unsymmetry in the dependencies of
$P^t,P^x$ at the beginning of {\tt CONLAW1}
is necessary because of the unsymmetry in using (\ref{mvdnls1}),
(\ref{mvdnls2}) to substitute a first order $t$-derivative of $u$
by a second order $x$-derivative.

(\ref{a2}): 1 condition in 22 variables $(t,x,u,v,\ldots,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,\ldots,u^{(2)},v^{(2)})$
and their 1st order derivatives, and furthermore
functions $Q^1,Q^2$ of 10 variables $(t,x,u,v,\ldots,u_{xxx},$ $v_{xxx})$.

(\ref{a3}): 2 coupled conditions in 14 variables $(t,x,u,v,u_x,v_x,\ldots,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,\ldots,u_{xxx},v_{xxx})$
and their 1st and 2nd order derivatives.

The following are general features of equations (\ref{a1})-(\ref{a3}).\\
Equ.\ (\ref{a1}) is of first order and therefore only highest
order $u$-derivatives which are not substituted due to $0=\Delta$
are not variables to the $P^i$ and can be used for direct separation.
Equ.\ (\ref{a1}) therefore is only weakly overdetermined with
the application of integrability conditions playing an important role.
A general problem with computing a differential Gr\"{o}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, equ.s (\ref{a3}) 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 overdetermined 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 equ.s (\ref{a2})
the $P^i$ depend initially on all jet-variables
(apart from highest order $u$-derivatives), also those
substituted through $0=\Delta$ on which the $Q^{\mu}$
do not depend. On the other hand the $Q^{\mu}$ do 
depend on highest order $u$-derivatives initially.
The efficiency in solving (\ref{a2}) 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 \cite{CRACK1}.

To solve the overdetermined system of all three approaches, all
techniques are used, only some are used more often in one approach
than in the other.

There is another issue.
If the order of derivatives w.r.t.\ different variables differs,
like, for the Harry Dym equation $0 = u_{t} - u^3u_{xxx}$,
then it matters whether this equation is used to do substitutions
$u_{t}=u^3u_{xxx}$ or $u_{xxx}=u_{t}/u^3$. Substituting $u_{t}$
gives a lower increase in complexity when successively higher
order ans\"{a}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 CLs compared to substituting $u_{xxx}$.
As equ.s (\ref{a3}) involve already higher order $u$-derivatives,
a further increase could explode the size of (\ref{a3}) even more.

Another relation between (\ref{a2}) and (\ref{a3}) is that one could
look at (\ref{a3}) as resulting from a differential-Gr\"{o}bner-Basis
calculation done with (\ref{a2}), with 
the aim to eliminate the $P^i$ first. It is of course more efficient to
exploit knowledge of the structure of (\ref{a2}) and to apply the Euler
operator to write down (\ref{a3}) directly rather than to do the differential
Gr\"{o}bner Basis calculation step by step with (\ref{a2}). On the other
hand {\tt 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 {\tt ODESOLVE} \cite{MM}).
For a concrete problem
it is very likely that there exists a quicker way to solve (\ref{a2})
than to eliminate at first all $P^i$. The question which of the
{\tt CONLAW} programs is more effective depends on the effectiveness of different
submodules of the program {\tt CRACK} which solves (\ref{a1})-(\ref{a3}).
With the current version of {\tt CRACK} (Jan.\ 1998),
programs {\tt CONLAW1/3} are better for simpler CL problems
and {\tt CONLAW2} is better for larger problems.

\section{Syntax of {\tt CONLAW}}
\noindent
{\it Example:} The input to find CLs with $Q$ of order 0-4 for the MVDNLS
equations (\ref{mvdnls1}),(\ref{mvdnls2}) is
\begin{verbatim}
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, {}, {}});
\end{verbatim}
In {\tt REDUCE} lists are enclosed in \verb+{ }+.
The input of {\tt CONLAWi} (i=1,2,3) consists of two lists, 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 {\tt CONLAWi} is a list
that specifies the CLs to be computed. Its first two elements are the 
minimum and maximum order of $P^1$ in the case of {\tt CONLAW1} and the
order of $Q^{\mu}$ in the case of {\tt CONLAW2/3}. The third element is
{\tt t} or {\tt nil} and specifies whether the CL 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^{\mu}$ 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 {\tt CONLAW} files.

\section{Summary}
Supplied with subroutines to fix gauge freedom in differential expressions
the programs {\tt CONLAW1/2/3} proved to be a efficient tool for the
computation of CLs of differential equations. Compared with other
programs, a list of which and a short description is given in \cite{GH1},
the programs {\tt CONLAWi} show the following new features:
\begin{itemize}
\item
By solving systems of overdetermined differential equations
it is possible to find CLs with non-polynomial, even non-rational
$P, Q$.
\item
It is possible to find CLs with an explicit dependence of $P, Q$ on
the independent variables.
\item
There is no limit on the number of DEs nor the number of independent
variables to be investigated for CLs other than a limit through the
complexity of computations.
\item
It is possible to determine values of parameters in the DE such that 
CLs exist (examples in \cite{TW}). 
\item
For each of the programs {\tt CONLAWi} 
an ansatz for $P^i$ and/or $Q^{\mu}$ can be input to specify
CLs to be calculated. 
\end{itemize}
Compared with the program of G\"{o}kta\c{s} and Hereman, {\tt CONLAW}
is able to find more general CLs and to make a definitive statement
if local CLs do not exist and the order is not too high to complete
the computations. 

The strength of the program described in \cite{GH1} is to get 
sometimes higher in the order that still can be handled
by concentrating on polynomial CLs having to solve
algebraic systems for coefficients of a polynomial ansatz.
They were also able to extend applicability to differential-difference
systems \cite{GH2}.

The comparison of the three approaches (\ref{a1})-(\ref{a3}) 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 overdetermined PDE systems directly
({\tt CONLAW1/3}) and an approach to derive integrability 
conditions by applying extra information about the structure of
the PDE system ({\tt CONLAW2}).

The programs including a manual and a test file
are available via ftp from {\tt lie.maths.qmw.ac.uk},
directory {\tt pub/compalg}.
%Programs are available from TW and 
The package will be submitted to the REDUCE network library.

\section{Acknowledgement}
The authors want to thank Malcolm MacCallum for comments on a first
version of the manuscript. Further, the support
at a research visit of TW at CAN Netherlands/Amsterdam
and discussions with Jan Sanders and Willy Hereman
and multiple visits at the Konrad Zuse Institute/Berlin 
and discussions with Winfried Neun are gratefully acknowledged.
Frank Verheest, Willy Sarlet, Micheline Musette and the Relativity group 
at Hall University are thanked for suggesting PDEs to test the code.

\begin{thebibliography}{99}
\bibitem{AB} Anco, S.C.\ and Bluman, G.\ (1997). 
        {\it Direct Construction of Conservation Laws from Field Equations}
        Phys.\ Rev.\ Let.\ {\bf 78}, no 15, 2869-2873.
\bibitem{GH1} G\"{o}kta\c{s}, \"{U}. and Hereman, W.\ (1996) {\it Symbolic
        Computation of Conserved Densities for Systems of Nonlinear
        Evolution Equations.} Journal of Symbolic Computation, 
        {vol. 24}, pp. 591-621 (1997).
\bibitem{GH2} G\"{o}kta\c{s}, \"{U}., Hereman, W.\ and Erdmann, G.\ (1997) 
        {\it Computation of conserved densities for systems of nonlinear
        differential-difference equations.}, 
        Physics Letters A, {vol. 236}, pp. 30-38 (1997).
%\bibitem{GH1} G\"{o}kta\c{s}, \"{U}. and Hereman, W.\ (1998)
%{\it Computation of conserved densities for nonlinear lattices}, 
%Physica D (1998) in press.
\bibitem{MM} MacCallum, MAH.\ (1988) {\it An Ordinary Differential Equation 
        Solver for REDUCE} Proc.\ ISSAC`88, Springer Lect.\ Notes in 
        Comp.\ Sc., 358, 196-205.
\bibitem{PO} Olver, P.J.\ (1986). {\it Applications of Lie Groups to
        Differential Equations.} Grad.\ Texts in Math.\ Berlin:
        Springer-Verlag.
\bibitem{NS} Willox, R., Hereman, W.\ and Verheest, F.\ (1995).
        {\it Complete Integrability of a Modified Vector Derivative
        Nonlinear Schr\"{o}dinger Equation.} Physica Scripta {\bf 52}, 21-26.
\bibitem{CRACK1} Wolf, T.\ and Brand, A.\ (1992)
        {\it The Computer Algebra Package CRACK for Investigating
        PDEs}, manual + software in the {\tt REDUCE} network library.
\bibitem{CRACK2} Wolf, T.\ (1996) {\it The program CRACK for solving PDEs
        in General Relativity} in {\it Relativity and Scientific Computing}
        Eds.\ F.W.Hehl, R.A.Puntigam, H.Ruder, Springer, 241-257.
\bibitem{TW} Wolf, T. To be published elsewhere.

\end{thebibliography}

\end{document}