# Long guide to the STANDARD FORM Package

LONG GUIDE TO THE STANDARD FORM PACKAGE
RELEASE 2:  AUGUST 27, 1993

Gregory J. Reid and Allan D. Wittkopf

REQUIREMENTS:
Access to Maple: Version 5 (releases 1 or 2) or later version.
The programs in the Standard Form Package
which is available by anonymous ftp to:
math.ubc.ca
(i.e. after user type anonymous).  Then change directory to
pub/reid/standard
in that directory for instructions.

ACKNOWLEDGEMENTS:

We wish to especially thank: Cameron Connell who substantially contributed
to the writing of this manual (especially to sections 3.6 and 6);
David Weih (who provided the preversions of the coordinate
transformation and interface routines); Alan Boulton who
was an indispensable participant in many lively discussions about the package
and who contributed the initial interface to the Liesymm package;
and Ian Lisle whose beautiful moving frame algorithm is an inspiration
for future developments.

GJR gratefully acknowledges the generous support of a grant from
the Natural Sciences and Engineering Research Council of Canada.

We also wish to thank:  Satish Reddy who
performed the first computer algebra experiments back in 1987;
Carl Wulfman for asking interesting questions which first turned
GJR's mind to the direction of applying formal operations to differential equations to expose their properties; Sukeyuki Kumei for initial
discussions; Beth Filliman for a test translation
of the Macsyma prototype program into Maple; George Bluman for support and discussions, and providing systems of PDEs which were intractable enough
that the development of such a formal package was kept on practical lines; Patrick Doran-Wu and Alex Ma for helpful suggestions;
Rob Entwistle, Alistair Blachford, Tom Nicol and Bob Bajwa for expert computer support and Diana Kwan for her encouragement and support.

CONTENTS

1. INTRODUCTION
2. GETTING STARTED WITH THE PACKAGE
2.0 HOW TO USE THIS MANUAL
2.1 SIMPLE EXAMPLES
2.2 SYSTEMS OF ALGEBRAIC EQUATIONS
2.3 NOW YOU'RE READY FOR SYSTEMS OF DIFFERENTIAL EQUATIONS
2.4 IMPORTANT NOTE ABOUT OPTIONS FOR standard_form
3. MORE ADVANCED ASPECTS OF THE PACKAGE
3.1 HOW TO CHANGE THE ORDERING - the option janet_order
3.2 CALCULATION OF INITIAL DATA AND
TAYLOR SERIES SOLUTIONS OF SYSTEMS OF PDE
3.3 THE STRATEGY OPTION
3.3.1 DISPLAYING INTERMEDIATE INFO USING STRATEGY
3.3.2 CONTROLLING EQUATION FLOW USING STRATEGY'S KNOBS
3.3.3 THE STRATEGY OPTIONS INCONSOPT AND SIMPOPT
3.4 STORING INTERMEDIATE RESULTS AND THE USEREASY OPTION
3.4.1 STORING INTERMEDIATE RESULTS
3.4.2 SELECTING NICE EQUATION SUBSETS WITH USEREASY
3.5 CLASSIFICATION PROBLEMS - SYSTEMS OF PDE WITH UNSPECIFIED
FUNCTIONS OR PARAMETERS
3.5.1 LOW-TECH APPROACH
3.5.2 HI-TECH APPROACH - the class_eqns option
3.6 INTERFACE TO MAPLE'S DERIVATIVE NOTATION
3.7 INTERFACE TO MANSFIELD'S NOTATION FOR HER DIFFGROB PACKAGE
3.8 THE DIRECTED GRAPH OF A PDE SYSTEM
4. APPLICATION OF THE PACKAGE TO SYMMETRY ANALYSIS OF PDES
4.1 CALCULATION OF THE STRUCTURE OF LIE SYMMETRY ALGEBRAS
OF PDES FROM THEIR INFINITESIMAL DETERMINING SYSTEMS
4.1.1 SYMMETRIES OF THE HEAT EQUATION:  Ut = Uxx
4.1.2 A GROUP CLASSIFICATION PROBLEM: Uxx = Utt + H(U) Ut
4.2 INTERFACE TO THE MACSYMA PACKAGE SYMGRP.MAX
4.2.1 PRODUCING THE DETERMINING SYSTEM OF Ut = Uxx
4.2.2 PRODUCING THE DETERMINING SYSTEM OF Uxx = Utt + H(U) Ut
4.3 INTERFACE TO THE MAPLE PACKAGE SYMMETRY
4.4 INTERFACE TO THE MAPLE PACKAGE LIESYMM
5. INTEGRATION OF THE EQUATIONS
5.1 USE OF taylor_exp
5.2 INTERACTIVE INTEGRATION
6. CHANGING VARIABLES IN SYSTEMS OF PDE
7. TOOLBOX OF INTERACTIVE FUNCTIONS
7.1 Clearderivdep
7.2 polyd
7.3 sfdiff
8. STRATEGIES TO DEAL WITH MEMORY EXPLOSION
8.2 VARYING THE ORDERING
8.3 CLASSIFICATION VARIABLES AND DIFFERENTIAL EXTENSIONS
8.4 THE BIG GUNS - STORAGE AND USEREASY
9. NONTRIVIAL EXAMPLES
9.1 A BIG EXAMPLE: THE MHD EQUATIONS
9.2 A SMALL NASTY EXAMPLE: THE NONLINEAR TELEGRAPH SYSTEM
REFERENCES
INDEX OF FUNCTIONS AND OPTIONS

1.  INTRODUCTION

The Standard Form Package consists of procedures for obtaining
information about analytic systems of Partial and Ordinary
Differential Equations (PDEs and ODEs) using completely
algorithmic differentiation and algebraic processes.  These procedures do not
depend on the heuristic process of explicitly integrating
such systems.  The package can yield information about the number of
solutions possessed by a system (if any), canonical forms for such
systems and structural properties of symmetry groups
of differential equations.  Paradoxically such information sometimes leads to new ways of solving PDE systems [Reid91c], [Reid92b].
The algorithms of the package may also be used
as precursors for the application of more heurisitic
methods (e.g. canonically simplifying systems of PDEs so that
existing heuristic or numerical methods have a better chance of
solving such systems).  We have applied the algorithms of this
package for the most part to problems concerning symmetry analysis of differential
equations.  We stress however that the algorithms can be used to
obtain information about analytic systems of PDEs regardless of their origin.

The theoretical background of this package lies in
the formal theory of differential equations, the
algorithms of this theory and their application to
symmetries of PDEs developed in [Reid90, Reid91a, Reid91b,
Reid91c, Reid92a].
The general area of formal analysis of differential equations is
strewn with mathematics of every variety.  Amongst these:  the abstract
cohomological approaches of Spencer and his co-workers [Spenc69];
the algebraic work of Ritt and his descendents [Ritt50, Kol73];
Pommaret's geometric formal methods [Pomm78, Pomm90, Schu92];
the differential Buchberger algorithms of Carra-Ferro [Carra87],
Ollivier [Olliv91], and Mansfield [Mans91]; approaches through the
use of exterior differential systems [Cart45, BCGGG91, HaE71, HartTker91]; applications to Lie symmetry groups of differential equations [Reid90, Reid91b, Reid92a], [Sch92a, Sch92b]; Wolf's resultant algorithm [Wolf89].
This area lies at the intersection
of geometrical, algebraic, analytical and algorithmic techiques.
Formal theory concerns the existence of formal series solutions
of PDE systems. By a formal series solution of a system we mean a series that
satisfies the system assuming that all operations - differentiation and summation of series etc - can be carried out formally term by term.  This
theory has its roots in the works of Tresse [Tress94], Riquier [Riq10] and Janet [Jan20] and has been developed and applied by the workers above
and many others.

The heart of the Standard Form Package is
the function standard_form, which is a function for
canonically simplifying systems of PDEs and ODEs.
The simplified system has the same formal solutions as the
original system away from a certain degenerate set of points which
are described by a set of pivots returned by the algorithm.

One of the major ideas of the algorithm [Jan20], [Reid91a]
is to split the set of all derivatives of a
solution at a point in the space of independent variables of a given
PDE system into two disjoint classes:  parametric derivatives
(those derivatives whose values at the point can be freely ascribed),
and principal derivatives (those derivatives whose values at a point
are uniquely determined through the system once the values of the
parametric derivatives have been given).  The dimension of the space
of formal solutions at the point is equal the number of parametric
derivatives and can be determined by an algorithm which does not
depend on explicitly solving the system, once the system is in standard form.
The collection of parametric derivatives at a point is called the
parametric initial data for the system.
If the system is an analytic function of its independent, dependent variables
and derivatives at a point, the parametric initial data
is analytic and the ordering used in the calculation is of
certain type (e.g. dominated by total degree), then Riquier [Riq10, Thom29]
showed that the corresponding formal solution about the given point
is the unique analytic solution of the system.

One of the major obstacles in implementing formal methods is that
the size of the calculations grows explosively with the order of the
derivatives in the systems of PDEs to which it is applied.
While developing this package, we
have been primarily motivated by our efforts in solving real
problems, provided by our own interests, and those of other
applied researchers we have worked with.  Consequently, the
program has been developed, not in a vacuum, but in a practical
environment. We have attempted to make the existing code have enough flexibility to allow users to deal with the difficulties
that encountered in applied problems (see, in particular,
section 8, for ways of dealing with memory explosion).

The function standard_form is most useful for simplifying
systems of large numbers of linear PDEs, which are overdetermined
or have many redundant equations.  The strong point of the
package is that it is algorithmic for linear systems of PDEs
- no hit or miss heuristics are involved. The canonical
form is guaranteed after a finite number of steps. For nonlinear
systems the package is not fully algorithmic but may partially simplify
such systems.  Buchberger's algorithm has provided a method for canonically
simplifying polynomially nonlinear systems of algebraic equations.
Although valuable progress has been made for polynomially nonlinear systems of PDEs [Carra97, Mans91, Olliv91], as yet no fully effective
differential analogue of Buchberger's algorithm is known.

The procedure standard_form can be regarded as a generalization of
the Gaussian Reduction algorithm for systems of linear algebraic
equations to systems of linear PDEs (to which it reduces if only
linear algebraic equations are involved).

If you are dealing with a system of linear PDEs then the
procedure standard_form can in principal guarantee you:

(1) determination of whether the system is consistent or inconsistent
(i.e. has no formal or analytical solutions).

(2) reduction of the system to an ordered canonical form using the
procedure standard_form;

(3) determination of the the dimension of the solution space of the
system by using the procedure initial_data;

(4) determination of the Taylor series solution of the system to a
given finite degree using the procedure taylor_exp;

at points which do not belong to a certain degenerate set of points
described by a set of equations (the "pivots") which is also returned
by standard_form.

The second part of the package, involves powerful algorithms for the symmetry analysis of PDEs. In particular our packages enable the dimension and
structure of Lie symmetry algebras of differential equations to be
directly calculated from their determining systems, without needing to
integrate such systems [Reid91b, Reid92a].

The conventional approach to infinitesimal symmetries of differential
equations has been to
try to integrate the determining equations, explicitly determining the
basis generators of the algebra, and then substituting the results into the
commutation relations above to determine the structure constants
[BlK89, Olv86, Ovs82].

Since this process is based on the heuristic process of integration
it is not always guaranteed to succeed.

Our approach to finding the structure of Lie symmetry algebras of
PDEs is algorithmic:

(1) The determining system for the infinitesimals by one of the computer
packages designed for this purpose (see [Here93, Champ91]).

(2) Algorithm standard_form is used to reduce the determining system
to standard form.
(3) Application of initial_data to the standard form produced in Step (2)
algorithmically determines the dimension of the algebra
of infinitesimal symmetries, without having to determine the
explicit form of these generators.
(4) Algorithm fin_com_table uses the standard form from Step (2)
and its initial data from Step (3) to determine the structure of the
of the finite-dimensional subalgebra of infinitesimal symmetries
corresponding to the parameters in the finite initial data.
This algorithm depends only on differentiation and substitution.
(5) taylor series approximations of the infinitesimals or basis generators
can be found to any given finite degree using algorithms taylor_exp, or
finite_generators respectively.

All of Section 2 is essential reading.  To get started, those
interested in symmetry analysis of PDEs should carefully cover
all of sections 2 and 4.1.

2.  GETTING STARTED WITH THE PACKAGE

2.0 HOW TO USE THIS MANUAL

The main intention of this manual is to provide a low-level,
example-based guide to our package, which assumes a minimum
of mathematical and computer-algebra background, but tries to
sneak in examples which illustrate some of the ideas behind the
package.  Consequently the flavour is loose and sacrifices
precision in favour of analogy.  A computer-oriented Maple-styled
short guide is planned together with works that concentrate on the
mathematics involved.
(or cut and paste) the examples we present (or a similar example
of your choice) into a Maple session.  We encourage experimentation.
Many examples are simple enough to work by hand.  Nontrivial examples
are given in Section 9.  We welcome your suggestions,

The aim of the manual is to provide easy continuous reading.
It is written in askey so that the find option of your editor
can be used to
locate all occurrences of that topic
(you may find it helpful to first look at the table
of contents, or at the index of functions and options at the end
of the manual).

2.1 SIMPLE EXAMPLES

EXAMPLE A

Consider the system of nonlinear PDEs:

y Zxxx - x Zxyy  =  Zyy - y Zy

2     2    2
2 y x Zxxx Zxyy + x Zxxx + x y Zxyy  =  0

2    2
y Zxyy - x W + 2 x  y Z  =  0

2    2
Zyy - y Zy  + 2 x  y W  =  x W

where the dependent variables W and Z are functions of the
independent variables x and y, and Zxxx denotes the third partial
derivative of Z with respect to x etc.

We start a Maple session (Maple 5.0 or higher
version) and read in the Standard Form Package:

(Maple command lines begin with > and are terminated with ; or :).
To use the programs we input the system as a list of equations with
independent variables x1,x2,x3, etc. and dependent
variables V(1), V(2), ... etc.  So denoting
x by x1, y by x2, W(x,y) by V(1) and Z(x,y) by V(2)
the equations above are written as a list:

> sys:= [x2*V(2,1,1,1)-x1*V(2,1,2,2) = V(2,2,2)-x2*V(2,2),
2*x2*x1*V(2,1,1,1)*V(2,1,2,2) + x1*V(2,1,1,1)^2 + x1*x2^2*V(2,1,2,2)^2 = 0,
x2*V(2,1,2,2) - x1*V(1) + 2*x1^2*x2*V(2)^2 = 0,
V(2,2,2) - x2*V(2,2) + 2*x1^2*x2*V(2)^2 = x1*V(1)]:

where we have terminated the command line with : instead of ; to suppress
the list sys from being displayed.
Derivatives are denoted by using a repeated index notation specific
to the program.  For example, V(2,2) denotes the x2 partial of V(2), and,
V(2,1,2,2) denotes the partial derivative with respect to x1
of the x2 partial of the x2 partial derivative of V(2)
(unlike Maple's usual notation for differentiation).
Notice that Maple is by default, case sensitive, so that x1 is
distinct from X1 etc.

The program depends on a complete ordering  which is defined on derivatives.
By default the ordering is the following standard one, denoted >s:

(1) First ordered by total order of the derivative.

(2) Then lexiographically by derivative:
(d/dx1 >s d/dx2 >s ... >s d/dxm )

(3) Then by function number (V(1) >s V(2) >s ... >s V(n))

Other orderings are possible and assigning different orderings can be useful
for reasons of efficiency or for forcing the system or its equations to take
certain pleasant forms [Reid92b].

Hence, in the standard ordering, the derivatives in our example are ordered

V(2,1,1,1) >s V(2,1,2,2) >s V(2,2,2) >s V(2,2) >s V(1) >s V(2)

There is a unique derivative (or dependent
variable) which is largest in each equation of a system
with respect to the above ordering.

To reduce this system to standard form we use the procedure
standard_form:

> sf:= standard_form(sys);

and the resulting output is:

The system parameters are as follows:

There are 2 independent variables:

[x1, x2]

There are 2 dependent variables:

[V(1), V(2)]

There are no classification variables:

sf := [V(2, 1, 1, 1) = 0, V(2, 1, 2) = 0, V(2, 2, 2) = x2 V(2, 2),

2
V(1) = 2 x1 x2 V(2) ]

The standard form above will have the same formal solutions as the
original system sys at points where none of the pivots vanish.
The pivots for the calculation are stored in a global variable _pivs
and are:

> _pivs;

2
2                x2                     x1 (x2  + x1 + x2)
[x1 = 0, - x2  - x1 = 0, - -------- = 0, - x2 = 0, - ------------------ = 0]
2                               x2
x2  + x1

Consequently the calculation is valid at all points (x1,x2) with
x1#0, x2#0, x2^2+x1#0, x2^2+x1+x2#0.  The pivots e.g. x1, -x2^2 - x1 etc
above are quantities which appear in divisions during the computation
(for nonlinear systems they may involve dependent variables
and their derivatives).  It may also yield valid results
at points other than these, but this conclusion can not be deduced from
the output standard form, sf, above.  The sfprogs command shortenpiv will
remove repeated factors from the list _pivs (in the order that they occurred):

> spivs:= shortenpiv(_pivs);
2             2
spivs := [x1, x2  + x1, x2, x2  + x1 + x2]

In our original notation the (considerably) simplified system is:
2
Zxxx = 0, Zxy = 0, Zyy = y Zy, W = 2 x y Z

Once a system of PDEs sys is in standard form
then its initial data can be calculated using algorithm initial_data.
As mentioned in the Introduction the initial data consists of a set
of assignments of values to parametric derivatives at an initial point
which is denoted by X=(X1,X2,...,Xm).
In particular if V(r) is one of these parametric derivatives, with
r a multiindex, r=(i,j,...k),
and b is the value of this derivative at the initial point
x=X (i.e. x1=X1, x2=X2, ... , xm=Xm) then the corresponding initial
condition has the form:

V(r)|      =  b
x=X

For simplicitly of display we will write the above condition simply

V(r) =  b

(with evaluation at x=X being understood) and denote the different arbitrary
constants (b's) by a1, a2, a3, ... etc.
For example applying the procedure initial_data to the standard form of the
system above gives:

> id:= initial_data(sf);

yielding

id := [[V(2) = a1, V(2, 2) = a2, V(2, 1) = a3, V(2, 1, 1) = a4], []]

The initial data list id above is equivalent to the initial conditions:

V(2)|     =  a1, V(2,2)|    = a2, V(2,1)|    = a3, V(2,1,1)|    = a4
x=X                x=X              x=X                x=X

posed at x=X (i.e. x1=X1, x2=X2) where a1, a2, a3 and a4 are arbitrary.
Or in terms of our original notation (with (X1,X2)=(x0,y0))

Z(x0,y0) = a1, Zy(x0,y0) = a2, Zx(x0,y0) = a3, Zxx(x0,y0) = a4

Since the initial data contains 4 arbitrary constants the general solution
of the system depends on 4 parameters and the solution space is 4-dimensional.
If the solution space is infinite dimensional then there would be an infinite
number of initial conditions of the form above.  Writing these down in a
finite manner requires the clever packing of these infinities into a finite
number of analytic functions (see the wave equation example B later in this
section and section 3.2).

Once we have computed the standard form of the original system, and computed the initial data, the taylor_series option will compute the formal taylor series
solution of the system to any finite degree N, with the command

> taylor_exp(sf,id,N)

where sf and id are the standard form and initial data of the system,
respectively.  The output will be the Nth degree Taylor polynomial for the
dependent variables about an initial point X = (X1,...,Xm).
In these polynomials hi = xi - Xi, i=1...,m and
their coefficients are functions of the the initial point X and the initial data.

To form the taylor series to degree 2 about a general initial point
(X1,X2) for the above example we use the command:

> ts:= taylor_exp(sf,id, 2);

which yields
2           2                                2    2
ts := [V(1) = 2 X1 X2 a1  + (2 X1 X2  a1 a2 + 4 X1 a1 a2 + 2 X1 X2 a2 ) h2

2               2
+ (4 X1 X2 a1 a3 + 2 X2 a1 ) h1 + (2 X1 a1  + 4 X1 X2 a1 a2) h2

2    2
+ (2 X1 X2 a1 a4 + 4 X2 a1 a3 + 2 X1 X2 a3 ) h1

2
+ (4 X1 X2 a2 a3 + 2 a1  + 4 X2 a1 a2 + 4 X1 a1 a3) h1 h2,

2               2
V(2) = a1 + a3 h1 + a2 h2 + 1/2 a4 h1  + 1/2 X2 a2 h2 ]

where h1 = (x1 - X1) and h2 = (x2 - X2) (note that the initial data must be
in the format presented by the procedure initial_data, so that a1, a2, ... are
not given particular values).  So the taylor series solution
about x0=0, y0=0 to order 2 for our original system is

2
W(x,y) = 2 a1  x y + 0(3)

2
Z(x,y) = a1 + a3 x + a2 y + a4 x  /2 + 0(3)

Note that the coefficients are expressed entirely in terms of
initially specified information: the initial point, X and the
initial conditions, a1,...a4.  Hence, given the initial data, we
can always compute the taylor series solution to any finite
degree in terms of these values.  In general the initial data
specified by algorithm initial_data is sufficient to uniquely determine
formal solutions of systems of PDEs. For those familiar with the
underlying theory, this is equivalent to the theorem by
Riquier and Janet [Riq10, Jan20] that if a system is in passive form
(i.e, all integrability conditions are satisfied) then specifying values
for the parametric derivatives uniquely specifies values for the
remaining principal ones.  If the standard form sf is an analytic
function of its arguments and the dimension of the solution space
is finite then this series will converge to a unique analytic solution
of the system in some neighbourhood of the initial point.

EXAMPLE B
Consider the wave equation, Uxx - Utt = 0 written in Maple's usual derivative
notation:

> msysb:= [diff(U(t,x),t,t) - diff(U(t,x),x,x) = 0];

It can be converted to the repeated index notation used by the program
by hand or by using the program described in Section 6.
Applying the program of Section 6 (which works in Maple 5 release 2 or later):

> sysB:=map2rep(msysB,[u],[t,x]);

sysB := [V(1, 1, 1) - V(1, 2, 2) = 0]

Taking standard_form of this equation gives:

> sfB:= standard_form(sysB);

sfB := [V(1, 1, 1) = V(1, 2, 2)]

and initial data

> idB:= initial_data(sfB);

idB := [[], [V(1) = f1(x2), V(1, 1) = f2(x2)]]

This corresponds to the usual initial data for the wave equation

U(t0,x) = f1(x), Ut(t0,x) = f2(x)

Notice how this really corresponds to an infinite sequence of initial
conditions of the form V(r)|    =  b:
x=X

U(t0,x0) = f1(x0), Ux(t0,x0) = f1'(x0), Uxx(t0,x0) = f1''(x0), ...

Ut(t0,x0) = f2(x0), Utx(t0,x0) = f2'(x0), Utxx(t0,x0) = f2''(x0), ...

so the program has a clever of packing the infinite number of initial conditions
into the functions f1(x) amd f2(x).

In general the initial data id is presented as a list which has 2 sublists.

The first entry of id (i.e. the list id[1]) contains a list of a finite number
(J say) of initial conditions each of form

V(r) = aj, j=1,...,J,

The second list contains K (infinite) initial conditions each of form:

V(s) = fk(y) where y is a sublist of [x1,x2,...,xm]

If the second sublist is empty (i.e. id[2]=[]) then the solution only
depends on the finite number of parameters a1, a2, ... ,aJ.
A more detailed description (especially of the infinite data case) is given
in section 3.2.

2.2  SYSTEMS OF ALGEBRAIC EQUATIONS

In this section we apply the programs to some simple systems of
algebraic equations.  Maple's solve command
and Grobner package are far more suitable for treating such systems of algebraic
(non-differential) equations than standard_form.
We consider this case first to
offer the user some informal insight into the workings of our package,
and how it deals with nonlinear equations.
We also explain the options divpiv and the global list _pivs, which records
pivoting operations of the package.  This discussion should be helpful for understanding the differential equations case.

Example A

We enter a system of linear algebraic equations.  As usual the equations must
be written in terms of dependent variables V(1),V(2),... which are assumed to
be functions of the (independent) variables x1,x2, ... that appear in the
equations.  So the system below is regarded as a linear system for
V(1),V(2),V(3),V(4) which are functions of x1, with k being a parameter.

> sysA:= [ 5*k*V(1) - 2*V(2) = 3*V(3)-V(4),
6*k*V(1) + 5*V(2) - x1*V(4) = 0]:
We now find the standard form:

> sfA:= standard_form(sysA);

(- 5 + 2 x1) V(4)    15  V(3)
sfA := [V(1) = 1/37 ----------------- + ---- ----,
k            37    k

18
V(2) = (6/37 + 5/37 x1) V(4) - ---- V(3)]
37

The standard form algorithm solved for those equations that are highest in
the standard ordering, as described in 2.1.  Since V(1) >s V(2) >s V(3)
>s V(4), it solved the two equations for V(1) and V(2).

The initial_data algorithm yields:
> idA:= initial_data(sfA);

idA := [[], [V(3) = f1(x1), V(4) = f2(x1)]]

so that V(3) and V(4) are arbitrary functions of x1, which makes sense
since there are no restrictions on V(3) and V(4) in the above standard
form.

The input system contained a parameter k so we are dealing with
a class of equations and
the standard form obtained may not be valid for exceptional values
of k.  The way the standard form was achieved was through the
standard gaussian elimination algorithm.  Thus the program
first solved the first equation for V(1) (so that k arose as a pivot),
subsituted into the second equation, and then solved the second equation
for V(2).  However, if a pivot is equal to zero, then the algorithm has commited an error, since it has divided by it.  Such cases require separate treatment.
To detect such exceptional classification cases the pivots are stored in a
global list called _pivs.

> _pivs;

[5 k = 0, 37/5 = 0]

So we see that exceptional cases occur only if any of the pivot equations
are satisfied.  Standard Forming the only exceptional case k=0

> sfA1:= standard_form(subs(k=0, sysA));

yields

sfA1 := [V(2) = 1/5 x1 V(4), V(3) = (1/3 - 2/15 x1) V(4)]

with initial data

> idA1:= initial_data(sfA1);

idA1 := [[], [V(1) = f1(x1), V(4) = f2(x1)]]

and pivots

> _pivs;

_pivs :=  [-2 = 0, -15/2 = 0]

so that no further branching is possible.
In the PDE case the program has sophisicated classification capabilities
and these are discussed in Section 3.5.

Example B

We consider the system of nonlinear algebraic equations:

> sysB := [ x1*V(1) - log(x2)*V(2)^2 = 0,
V(3)^2 - 2*V(3)*V(4) + V(4)^2 = 0,
V(3)^4*V(2)^2 - 2*x4*log(V(4))*V(3)^2*V(2) + x4^2*log(V(4))^2 =0]:

> sfB:= standard_form(sysB);

yields:

******************************************************
********** The system is not fully reduced ***********
** The non-linear equations are in the second list: **
******************************************************
2
ln(x2) V(2)
sfB := [[V(1) = ------------, V(3) = V(4)],
x1

4     2                     2          2         2
[0 = - V(4)  V(2)  + 2 x4 ln(V(4)) V(4)  V(2) - x4  ln(V(4)) ]]

with pivots:

> _pivs;

[x1 = 0, 1 = 0]

The standard ordering >s on the dependent variables, as given in 2.1, is:

V(1) >s V(2) >s V(3) >s ...

which allows us to identify the largest variable (wrt >s)
in any equation.  If an equation is of form
N
p (V(i) - r)   = 0      (*)

with largest variable V(i)
where p is a function only of the independent variables
and r is not a function of V(i) (N a nonzero positive integer),
then procedure standard_form will solve the equation for V(i) to obtain

V(i) = r           (**)

Thus in Example B it solved the first equation for V(1)
and since the second equation
was (V(3) - V(4))^2 = 0, it solved this equation to obtain V(3) = V(4).
It did not solve the last equation, even though this equation has form (*)

4 /        x4 ln(V(4)) \ 2
V(3)  | V(2) - ----------- |    = 0
|              2     |
\          V(3)      /

since p (= V(3)^4) was a function of a dependent variable.
Consequently the system was not fully reduced.  ( Algorithm initial_data
can not be applied here since the system is not in standard form.)
If the option divpiv
is added to the argument list of standard_form then the procedure will
solve equations of form (*) to obtain (**) even if p depends on
some of the dependent variables and their derivatives, provided they do not
depend on V(i) or derivatives of V(i).

Thus for our example:

> sfB1:= standard_form(sysB, divpiv);

fully reduces the system

2         2
ln(x2) x4  ln(V(4))          x4 ln(V(4))
sfB1 := [V(1) = --------------------, V(2) = -----------, V(3) = V(4)]
4                         2
V(4)  x1                  V(4)

with pivots:

>  _pivs;

4
[x1 = 0, V(3)  = 0, 1 = 0]

alerting us that V(3) = 0 is a special case that needs separate treatment.
Since the system is now in standard form, algorithm initial_data can
be applied:

> id:= initial_data(sfB1);

which yields:

id :=  [ V(4) = f2(x1,x2,x3,x4) ]

It's easy to see that allowing division by such pivots is a dangerous
operation and so we've made divpiv an option.  By returning a
system of unsimplified leading nonlinear equations the opportunity is
open to then use other packages to simplify these nonlinear equations
(e.g. Maple's algebraic Grobner package or Mansfield's package [Mans91]).

Example C.

Consider

> sysC :=  [2*V(1) - a*x1 = 0, V(1) - x1 = 0];
> standard_form(sysC);

Error, (in maxder) Inconsistent equation in maxder:, 0 = 1/2*x1*(a-2)
which indicates that if a#2 then the system is inconsistent and has no solution
(no relations are allowed amongst the independent variables).  The case a=2
has to be considered separately.

Example D.

Consider

> sysD:=  [ V(1) - 4*x1 = 0, V(2)/(V(1)-4*x1) = 0]:

> sfD:= standard_form(sysD);

Error, (in clearderivdep) division by zero

which indicates the system has no solutions for obvious reasons.

2.3  NOW YOU'RE READY FOR SYSTEMS OF DIFFERENTIAL EQUATIONS

The algorithm standard_form can be applied to systems of PDEs
with independent variables x1,x2, ... , xm and dependent variables
V(1),V(2), ... , V(n).  We give an informal description through
examples of some of the major aspects of the program.

Example A.

Consider the system of PDEs
Z - Wxy + y Wyy - xW = 0             (1)
Zx - Wxxy + y Wxyy - W - x Wx = 0    (2)
for W(x,y) and Z(x,y).  The largest derivative in equation (1)
is Wxy and the largest derivative in equation (2) is Wxxy.

The first stage of the standard_form algorithm reduces a system as
much as possible by a gaussian-elimination type process to a solved-form
consisting of equations of the form

derivative = rhs

where the rhs derivatives are lower in the ordering than
the derivative on the lhs.  At the same time the program
records the pivots in the global list _pivs.
After this stage
no derivative can appear explicitly on both the LHS and RHS of the system,
and derivatives only appear once on the LHS.

For the above system this triangular form is accomplished by solving the
first equation for Wxy and the second equation for Wxxy
giving:
Wxy  = Z + y Wyy - x W                (1)'
Wxxy = Zx + y Wxyy - W - x Wx         (2)'
This is the maximum elimination possible if we regard Wxy and
Wxyy as completely different algebraic variables.
However they aren't independent since Wxxy = (Wxy)x
and the system above is equivalent to (has the same solutions as) the system:

(1),   (2) - d/dx1 (1).

The second equation is an identity, so the system reduces to
just one equation:
Wxy  = Z + y Wyy - x W
And this is the standard form of the system.

To carry out the calculation in Maple we read the package sfprogs

and let W=V(1), Z=V(2), x=x1, y=x2, so that the system for example A is:

> sysA:= [ V(2) - V(1,1,2) + x2*V(1,2,2) - x1*V(1) = 0,
V(2,1) - V(1,1,1,2) + x2*V(1,1,2,2) - V(1) - x1*V(1,1) = 0];

and then

> sfA := standard_form(sysA);

yields

sfA:= [ V(1,1,2) = V(2) + x2 V(1,2,2) - x1 V(1) ]

agreeing with the result obtained above.

Example B.

In general such differential substitutions, as in the previous
example will lead to some nontrivial relations
amongst the derivatives of the system.  These new equations can be added to
the original equations, and the resulting system can again
be triangularized.  After a finite number of steps
the system will be obtained in a form where no derivative on the
LHS of the system is a derivative of any other derivative on the
LHS of the system.  For example:
Px = P
Py = x P
where P=P(x,y).
However this system may still not be fully reduced.  We see that
for consistency we need the integrability condition:
(Px)y - (Py)x  = 0
(P)y  - (x P)x = 0 (modulo system)
Py   - P - x Px = 0
We substitute for Px and Py from the first two equations
and the new equation reduces to P = 0.  We add this to the system.
But with this equation, the whole system reduces to
P = 0
and this is the standard form of the system.
In maple we have

> sysB:= [V(1,1)=V(1), V(1,2)=x1*V(1)];

for the above system (with x=x1, P=V(1)) and applying standard_form

> sfB:= standard_form(sysB);

yields

sfB:= [ V(1) = 0 ]

as above.

We encourage the reader to experiment with their
own examples to get a feel for the algorithms.

Example C.

Consider the linear system of PDEs:

Uxt - r x t = 0,  Uxx - Ux = 0

Where U = U(x,t).

Putting U = V(1), x=x1, t=x2 then

and the system is entered as:

> sysC:= [ V(1,1,2) - r*x1*x2 = 0,
V(1,1,1) - V(1,1) = 0];

Applying standard_form:

> sfC1:= standard_form(sysC);

gives the message:

Error, (in bigstandard) System is inconsistent:, r*x2*(-1+x1) = 0

Written in the original variables, the condition  rt(-1+x) = 0
has resulted from taking the integrability condition (Uxt)x=(Uxx)t,
after the algorithm has reduced it to the solved form
Uxt = r x t ,  Uxx = Ux
So the system is inconsistent (i.e. has no solutions) unless r=0.

When r=0

> sfC2:= standard_form(subs(r=0, sysC));

sfC2 := [V(1, 1, 1) = V(1, 1), V(1, 1, 2) = 0]

and the initial data is

> idC2:= initial_data(sfC2);

idC2 := [[V(1, 1) = a1], [V(1) = f1(x2)]]

so that the general solution depends on one arbitrary function
of x2 and one arbitrary constant (i.e. the solution space
is infinite-dimensional).  In particular the initial data is

V(1,1)|           =    a1   and    V(1)|     =  f1(x2)
x1=X1,x2=X2                      x1=X1

where a1 is an arbitrary constant and f1 is an arbitrary function of x2.
Written in our original notation these initial conditions are:

Ux(x0,t0) = a1  and U(x0,t) = f1(t)

where (x0,t0) is a fixed point.  Also

> _pivs;
[-1 = 0, 1 = 0]

so that there are no nontrivial pivots.

The taylor series expansion of the solution to third order is

> ts:= taylor_exp(sfC2, idC2, 3);
/  d         \               2
ts := [V(1) = f1(X2) + a1 h1 + |----- f1(X2)| h2 + 1/2 a1 h1
\ dX2        /

/   2         \                        /   3         \
|  d          |   2            3       |  d          |   3
+ 1/2 |------ f1(X2)| h2  + 1/6 a1 h1  + 1/6 |------ f1(X2)| h2 ]
|    2        |                        |    3        |
\ dX2         /                        \ dX2         /

which is the taylor solution about a fixed initial point (X1,X2)
giving the values of V(1) at X1 + h1, X2 + h2.  There are some nonzero
third order terms in the series so the series may not terminate
at this order (it would have terminated if there had been been no
third order terms - see Section 5.1).

Example D.

Consider the nonlinear system of PDEs for U=U(x,t)

2
Uxt = 0,  Ux Uxx - Ux  - Ut Uxx  + Ut Ux = 0

that is putting U=V(1), x=x1, t=x2

> sysD:= [ V(1,1,2) = 0,
V(1,1)*V(1,1,1)-V(1,1)^2 - V(1,2)*V(1,1,1)+V(1,2)*V(1,1) = 0];

> sfD:= standard_form(sysD);

******************************************************
********** The system is not fully reduced ***********
** The non-linear equations are in the second list: **
******************************************************

sfD := [[V(1, 1, 2) = 0],

2
[0 = - V(1, 1)V(1, 1, 1) + V(1, 1)  + V(1, 2)V(1, 1, 1) - V(1, 2) V(1, 1)]]

does not reduce the system completely to standard form.  Adding the
option divpiv forces standard_form to solve for V(1,1,1)
in the second equation and the system is completely reduced:

> sfD1:= standard_form(sysD, divpiv);

sfD1 := [V(1, 1, 1) = V(1, 1), V(1, 1, 2) = 0]

Great care must be taken when dealing with nonlinear systems of PDE.
In particular the pivots should always be examined to identify
special cases.  Here

> _pivs;

[-1 = 0, V(1, 1) - V(1, 2) = 0]

So a special case occurs if V(1,1) - V(1,2) = 0
(reflecting the fact that the sfD1 was obtained by the dangerous
operation of dividing by the pivot V(1,1) - V(1,2) in the second equation
to solve for V(1,1,1)).
To treat this case we append this equation to sysD

> sysD2:= [op(sysD), _pivs[2]];

sysD2 := [V(1, 1, 2) = 0,

2
V(1, 1) V(1, 1, 1) - V(1, 1)  - V(1, 2) V(1, 1, 1) + V(1, 2) V(1, 1) = 0,

V(1, 1) - V(1, 2) = 0]

and reduce the resulting system to standard form:

> sfD2:= standard_form(sysD2);

sfD2 := [V(1, 2, 2) = 0, V(1, 1) = V(1, 2)]

> _pivs;

[-1 = 0, 1 = 0, 1 = 0]

so there are no further special cases.

2.4 IMPORTANT NOTE ABOUT OPTIONS FOR standard_form

The function standard_form can be used with a list of options:

> standard_form(sys, optionlist);

You've already met one of those options (divpiv) and many more
will be described throughout the manual.

Any number of the options can be used in any order.
Effective use of the program will often entail using quite
a number of them.  For example

> standard_form(sys, class_eqns=[[V(4),3], V(4,3,3) = V(4)], store,
strategy = [1.75], divpiv);
where the 3 options above can be given in any order.
See the index of functions and options for a list of these options.

3.  MORE ADVANCED ASPECTS OF THE PACKAGE

3.1 HOW TO CHANGE THE ORDERING - the option janet_order

The standard form of a system depends on the ordering chosen.

For purposes of efficiency, or to help expose simple equations
in the system it is useful to be able to change the ordering.

The ordering on derivatives must be left invariant under
differentiation.  A wide class of orderings having this
property was given by Riquier [Riq10] (also see Janet [Jan20]).

Such orderings can be characterised by a matrix R with m+n
columns (m=#indep vars and n=#dep vars), with nonnegative integer
entries.  There are orderings beyond those of Janet - see
[CarraSitt93], [Rust93], [Weisp93].

To define such janet orderings we represent derivatives of V(p):

a1      a2          am
(d/dx1) (d/dx2) ... (d/dxm)  V(p)

by the corresponding vectors

[a1, a2, ... , am, 0,...,0,1,0,...,0]
(m+p)-th

(e.g. V(2,1,1,2) where m=2, n=2 is represented by [2,1,0,1]).

The usual lexicographic ordering >lex on vectors is defined by:

p=[p1,p2, ... , pk]  >lex  q=[q1,q2, ..., qk]
iff
the first nonzero pj - qj is positive, that is, the first
nonzero element of the vector p-q is positive.

The ordering on derivatives >  defined by a matrix R is given by:
R

a1      a2          am                 b1      b2          bm
(d/dx1) (d/dx2) ... (d/dxm)  V(p)  >   (d/dx1) (d/dx2) ... (d/dxm) V(q)
R
iff
T            T
R v   >lex   R w

where T denotes transpose and

v=[a1,...,am,0,...,1,...,0],  w=[b1,...,bm,0,...,1,...,0]
m+p                           m+q

The standard ordering given in Section 2.3 has ordering matrix

x1 x2 x3 ... xm V(1) V(2) V(3) ....  V(n)
1   1  1      1  0    0    0          0
1   0  0      0  0                    0
R =   0   1  0      ......                  0
0   0  1      ......                  0
0   0  0      ......                  0
0   0  0      1    ......             0
0   0  0      0  n   (n-1) (n-2)      1

Example A.

New ordering matrices are entered as list of lists option
to standard_form (i.e. not in maple's matrix or array format).

For the heat equation:

> sf:= standard_form([V(1,2)-V(1,1,1) = 0]);

yields

sf:= [V(1,1,1)=V(1,2)]

in the standard ordering.
To execute standard_form with nonstandard orderings these are
entered as an option to the standard_form procedure, in the form:

janet_order = list of lists
where each list has length m+n (m=#indeps n=#deps)

For example to reduce the heat equation in in the ordering defined by
the ordering matrix

x1 x2 V(1)
0  1  0
R =    1  0  0
0  0  1

(i.e. ordering first by the independent variable x2) the command is:

> sf:= standard_form([V(1,2)-V(1,1,1)=0], janet_order=[[0,1,0],
[1,0,0],
[0,0,1]]);
yields

sf:= [V(1,2) = V(1,1,1)]

Entering at least one row is enough to define a complete ordering.
If your entered ordering fails to distinguish between 2 derivatives
the program will use the default ordering to break the tie.

For example, above, we could have obtained the result above by:

> sf:= standard_form([V(1,2) - V(1,1,1) = 0], janet_order=[[0,1,0]]);

that is

sf:= [V(1,2) = V(1,1,1)]

The effect on systems of PDEs by using new orderings can
be dramatic.  The most useful (but computationally most
expensive orderings) are the lexicographical orderings
which order first by function number rather than order of derivative:
The first row of R is then

[0,0,...,0,perm]

where perm is some permutation of the integers {1,2,...,n}.

The lexicographic orderings have the effect of sequentially decoupling
a system of PDEs.  An interesting extension of the work in [Reid91c] would
be to investigate how orderings affect the directed graph of a
PDE system.

3.2 CALCULATION OF INITIAL DATA AND
TAYLOR SERIES SOLUTIONS OF SYSTEMS OF PDE

Once a system of PDEs sys is in standard form

> sf:= standard_form(sys);

then its initial data can be calculated using algorithm initial_data:

> id:= initial_data(sf);

The first entry of id (i.e. the list id[1]) contains a list of J
(finite) initial conditions each of form

V(r) = aj

which can be interpreted as:

V(r)|      =  aj
x=X

i.e. the derivative V(r) evaluated at the fixed point X=(X1,X2,...,Xm)
has the value aj.

The second list contains K (infinite) initial conditions each of form:

V(s) = fk(y) where y is a sublist of [x1,x2,...,xm]

which can be interpreted as:

V(s)|   = fk(y)
H

i.e. the derivative V(s) evaluated on the hyperplane H takes
the values fk(y) where fk is an arbitrary function of y and

H:= {(x1,x2,...,xm)| xj = Xj, for xj not in y,
xj arbitrary for xj in y}.

Example A.

For the heat equation Uxx - Ut = 0:

> sf1:= standard_form([V(1,1,1) - V(1,2) = 0]):

yields

sf1:= [V(1,1,1) = V(1,2)]

and

> id1:= initial_data(sf1);

yields the initial data

id1:= [[],[V(1) = f1(x2), V(1,1) = f2(x2)]]

that is
U(x0,t) = f1(t), Ux(x0,t) = f2(t)
and there are no finite initial conditions since the first list is empty.
To obtain the usual initial data for the heat equation
we have to use a nonstandard ordering (see 3.1)

> sf2:= standard_form([V(1,1,1)-V(1,2)=0], janet_order=[[0,1,0]]);

yields

sf2:= [V(1,2) = V(1,1,1)]

and

> id2:= initial_data(sf2);

gives

id2:= [[], [V(1) = f1(x1)]]

Which is the usual initial data for the heat equation
(i.e. U(x,t0) = f1(x)).

See 2.2 for an example where only finite initial conditions occur
and Section 2.3 Example C, for a case where both finite and
infinite initial conditions occur.

Once the standard form of a system of PDEs sys has been
determined and the initial data calculated from sf
then its taylor series expansion about a point (X1,X2,...,Xm)
in powers of h1(=x1-X1),...,hm(=xm-Xm) to an order N
can be calculated by algorithm taylor_exp from sf and id:

> sf:= standard_form(sys);

> id:= initial_data(sf);

> ts:= taylor_exp(sf, id, N);

Example B.

The taylor series expansion for the heat equation
to degree 4 is given by:

> ts:=  taylor_exp(sf2, id2, 4);

where sf2 and id2 have been calculated previously (above).

Note: at present taylor_exp though simple is very inefficient and it
would be worthwhile to develop strategies to increase its efficiency.

3.3 THE STRATEGY OPTION

The strategy option allows the user to control many details of the operation
of standard_form.  Invoking the strategy option involves entering
a list, at the same time as a call to the standard form program, whose entries
are new values for certain parameters that control the operation
of standard form.  These parameters include one that control the
amount of intermediate information that standard form displays as
it runs (infolevel below), and also basic decisions regarding how many new,
unclassified equations it should attempt to process in each
iteration (controlled by the values of p1 and p2 below).

To change any of the parameters in strategy from their default settings,
you enter a new strategy list as an argument to standard_form, in
the following way

> sf:=standard_form(sys,strategy=list)

where list is one of the 3 possibilities:

[infolevel]
[infolevel, p1, p2]
[infolevel, p1, p2, inconsopt, simpopt]

i.e. list has 1, 3 or 5 entries. Also infolevel, p1 and p2 have
numerical values and inconsopt, simpopt are boolean-valued.

default setting:
atrategy = [1.0, 1.0, 2.0, true, true]

3.3.1 DISPLAYING INTERMEDIATE INFO USING STRATEGY

The first entry of the strategy list, regardless of whether it contains
1,3 or 5 entries, controls the displaying of
intermediate information about the current run.

At any time during the running of standard form, the equations in the
system will be in one of 4 classes:

one-term eqns  Easy eqns  Leading Nonlinear eqns  Unclassified eqns
OT             Ez         NL                      UnC

The leading nonlinear equations are those equations which are nonlinear
in their highest derivative.

The program tries to take advantage of one term equations
as much as possible (since they always lead to simplification),
and then uses a criteria (the Easy-Unclassified choice)
for choosing an easy subset of the full
set of equations (the easy equations).  Nasty (because of their
length or nonlinearity) equations are stored in an unclassified equations
list.  Usually the algorithm has terminated when there are no
unclassified equations remaining.  The program has sophisticated strategies
for avoiding doing and redoing integrability conditions.
Integrability conditions are executed once at each iteration.

strategy = [infolevel]
(or [infolevel,p1,p2] or [infolevel,p1,p2,inconsopt,simpopt])
where:

infolevel < 1
no intermediate information will be displayed,
except the usual Maple memory and time reports.
1    <= infolevel < 1.25 (default setting)
standard_form will report the
number of equations in each class at the end of each basic iteration
1.25 <= infolevel < 1.5
standard_form will also report
the maximum length of equations in the Ez and UnC categories, and
the number of equations in the UnC category that have been classified
and put into the OT, Ez and NL categories.
1.5  <= infolevel < 1.75
in addition to the above info standard_form
gives a glossary of the symbols in the output and gives a more detailed
reporting of the steps of the algorithm
1.75 <= infolevel
in addition to the above info standard_form gives
detailed info about all its major steps (this is GR's favourite level
and a good way to undestand the algorithm).

Of course you can also use Maple's printlevel command.

3.3.2 CONTROLLING EQUATION FLOW USING STRATEGY'S KNOBS

The behaviour of the program can be changed by twisting two knobs;
the p1 knob and the p2 knob (that is by giving the parameters p1 and p2
different values) where:

strategy=[infolevel, p1, p2]

These (equation flow') parameters control the the quantity and
relative size of the equations being chosen from UnC and placed into
the easy list Ez.

Special Settings

default setting:
strategy = [1, 1, 2]                  (i.e. p1=1, p2=2)

Ultra-Conservative setting:
strategy=[infolevel, 1, 0.9]          (i.e. p1=1, p2=0.9)
Good for large or complicated systems (transfers 1 equation
at a time from UNC to Ez).

Ultra-HOG setting:
strategy=[infolevel, 100000, 100000]  (i.e. p1=100000, p2=100000)
Only use Ultra-HOG on very simple systems, when you're greedy
for a quick run (transfers everything from UnC to Ez).

The p1 knob
-----------
Our first parameter p1 controls an upper bound for the number of
equations being placed into Ez relative to the number already present.
This bound is defined as a function of p1 as follows:

# equations currently       Max # equations taken from UnC & put in Ez
in Ez (j)                  bound(j)
-----------------       ------------------------------------------
0 <=j <=15             10 x p1
16 <=j <=25              8 x p1
26 <=j <=35              6 x p1
36 <=j <=40              5 x p1
41 <=j <=45              4 x p1
45 <=j                   3 x p1

In our experience the above relationship has worked fairly well
for large systems, as the memory and time
required for simplification of larger systems behaves roughly as
the product of the number of equations in Ez and the number of new
equations.
Setting p1=10000 will mean that the number of equations currently
in Ez will be ignored, and the maximum number of equations
allowed by the p2 knob will be transferred from
UnC to Ez at each iteration.
Setting p1=0 means that only 1 equation will be transferred from
UnC to Ez at each iteration.
As the p1 knob is turned up the number of equations transferred from
UnC to Ez becomes less sensitive to the currrent number of equations in Ez.

Changing the way Ez equations are chosen from UnC by altering the numbers
(15,25, ...etc) in the above table - can be  accomplished by changing the
these numbers in the following piece of sfprogs code (par[7] is p1):

$$start of section of sfprogs code$$
#
# Here we select which Hard equations are to be put in Easy
#
t2:=array[1..4];
if nops(Hard)>0 then
if nops(Easy)<=15 then numeqn:=10*par[7];
elif nops(Easy)>15 and nops(Easy)<=25 then numeqn:=8*par[7];
elif nops(Easy)>25 and nops(Easy)<=35 then numeqn:=6*par[7];
elif nops(Easy)>35 and nops(Easy)<=40 then numeqn:=5*par[7];
elif nops(Easy)>40 and nops(Easy)<=45 then numeqn:=4*par[7];
elif nops(Easy)>45 then numeqn:=3*par[7];
fi;
$$end of section of sfprogs code$$

More extensive changes to the way equations are moved between lists
require great care (note: par[8] in sfprogs corresponds to p2);
but is an activity which warrants investigation and since it is central
to efficiency of the program.

The p2 knob
-----------
Our second parameter controls the flow of equations in two ways, and is
sensitive to the length of equations in UnC.
The length of an equation is given by Maple's length function and is a
simple measure of the complexity of an equation by how much memory is
required to store it.
The equations in UnC are placed in order according to their increasing
length.  If the length of the first equation in UnC is more than p2 times
longer than the length of the longest equation in Ez, then that equation
is the only equation added to Ez in that iteration.
If not, then equations of UnC are considered in order of increasing
length until either we have the maximum number of equations (as dictated by
the p1 knob) or the length of the next equation is more than p2 times longer
than the length of the first equation added to Ez.

Setting p2=10000 will make standard form ignore the lengths of
equations in UnC, and transfer the maximum number allowed by the
p1 knob.

This option is very useful for small complicated systems as the
simplification of complicated equations is put off until after the
simplification of the easier equations. This gives the easier
equations an opportunity to simplify the complicated (UnC) equations
before they are put in Ez.

A note of clarification.  In the process of moving equations from the UnC
list to the Ez list, standard_form will automatically move equations
that are not polynomials in their top derivatives into the NL list,
(example, sin(V(1,1,1))+V(1,1,1)=0), and these will not be counted
towards the maximum number of equations that can be transferred from
the UnC list, as determined by the p1 knob.  However, equations which
are polynomial in their top derivative will be transferred to the Ez
list, even if their degree is higher than 1, and all such equations will
be counted in the total number of such equations to be allowed.  However,
after reducing these equations, with respect to the other equations in
Ez, all those which are still polynomially nonlinear will be moved
from the Ez to the NL list, and standard_form will not transfer any
more equations from UnC to Ez, even if their are still some which are
short enough given the settings of the p2 knob.

Example 1:
strategy=[1.75, 1, 2]
p1=1, p2=2 (Default knob settings and infolevel set to its
highest value)
Current Max length in Ez: 117
Lengths in UnC (min to max):
110, 140, 151, 152, 240, ......
Current Number of equations in Ez: 7

In this case there are 7 equations in the Ez list at the
beginning of the iteration, so, by the above table
describing the operation of the p1 knob,
the maximum number of equations that will
be added to the Ez list is 10 x p1 = 10 x = 10.
But the fifth equation in UnC has length 240 > 2 x 110
= p2 x 110, and 110 is the length of the first equation
in UnC, therefore the first equation added to Ez, so
due to the second part of the p2 test only the first 4
equations, with lengths  110, 140, 151, 152, which are
all less than 2 x 110, are moved.

Example 2:
strategy=[1.75, 0.5, 2]
p1=0.5, p2=2
Current max length in Ez: 164
Lengths in UnC(min to max):
1310, 1476, ....
Current number of equations in Ez: 20

In this case there are 20 equations in Ez initially, so
the maximum number of equations that will be added is
8 x p1 = 8 x 0.5 = 4.  But the first equation in UnC has
length 1310 > 2 x 164 > p2 x 164, so since 164 is the longest
equation initially in Ez, only the first equation is moved,
by the first part of the p2 test.

Example 3:
strategy=[1.75, 0.5, 2]
p1=1, p2=2
Current max length in Ez: 54
Lengths in UnC(min to max):
58, 61, 62, 66, 75, 79, 80, 81, ......
Current number of equations in Ez: 43

In this case the maximum number of equations is 4 x p1
= 4 x 1 = 4, since their are 43 equations in Ez initially.
Since all of the first four equations in Unc have
lengths less than 2 x 54 = p2 x 54, all of them are chosen.

Special Settings:

Default Knob Settings: p1=1, p2=2

Ultra-Conservative (one equation at a time): p1=1, p2=0.9
With this setting the program will only move one equation at a time
into the Ez list, since the second equation is always
greater than 0.9 times the length of the first.
For example including strategy=[1.75, 1, 0.9] in the option list for
standard_form will set the infolevel at 1.75, and the knob settings
to the ultra-conservative settings of p1=1, p2=0.9.

Ultra-HOG (everything at once): p1=100000, p2=100000
Note: only use Ultra-HOG on very simple systems.

3.3.3 THE STRATEGY OPTIONS INCONSOPT AND SIMPOPT

The 4th and 5th elements of the strategy list must be boolean
valued, and their use is outlined in the table below; they have
rarely if ever been used.

strategy list element    value              result
4th (inconsopt)          true       standard_form will terminate if any
inconsistent equations are found.
false      standard_form will save the
inconsistent equation in the global
variable _incons, and continue

5th (simpopt)            true       standard_form will collect coefficents
of dependent variables prior to any
call of maples simplify procedure.
false      standard_form will apply maples
simplify procedure directly on the
equations, without first collecting
coefficients of the dependent variables

The default value for both of these parameters is true.

INCONSOPT: The 4th element of strategy is the inconsopt parameter.
When its value is true, the standard_form
will terminate as soon as an equation relating only the independent
variables, as illustrated in example C of section 2.2.3.  When the
inconsopt parameter has the value false, then any such equations that
standard_form finds will be stored in the global variable _incons.
Standard_form will then ignore the equation, and continue to perform the
algorithm on the remaining equations.  This allows the standard form to
be calculated when relations amongst the independent variables are
found that actually reduce to zero, but this fact is not recognized
by maple.  For instance, the equation sin(x1)^2+cos(x1)^2-1=0 reduces
to zero, but maple will not recognize this.  The standard form obtained
with inconsopt set to false is only valid if all the equations in incons
reduce to zero.

SIMPOPT: The 5th element of strategy is the simpopt parameter.
If it has its default value of true, then whenever maple simplifies equations, it will first collect the coefficients of any expressions involving the
dependent variables and their derivatives, and then will apply maples
simplify procedure.  If the simpopt parameter is set to false, then
standard_form will apply maple's simplify procedure on the equations
directly, without collecting coefficients first.

3.4 STORING INTERMEDIATE RESULTS AND THE USEREASY OPTION

3.4.1 STORING INTERMEDIATE RESULTS
If the option store is included in the argument list for standard_form, e.g.

> sf:= standard_form(sys, store);

it will store the status of the calculations at the end of each
iteration to a file called StandardStorage.m in your current directory.
This file will be in stored in maple's internal format (as indicated by
the .m suffix), so the file cannot be viewed as stored, but must be
read into maple.  To access or view the information in StandardStorage.m,
go through the following steps:

1) Start a new maple session.  If you want to view the information in
StandardStorage.m as it is updated, that is, while standard_form is
running, this will mean starting a new maple session independent of
the one that standard_form is running on, by, for example, starting
a new system shell.  Once this new session is started, the Standardform
package will be needed to read the file, so type

2) Now call the procedure readstored().  This returns the state of the calculations at the end of the last iteration (integrability condition
loop) in the form of a list.  For example, if we call the list L:

then the list L will have 4 elements as follows:
L[1] is the list of one term equations
L[2] is the list of easy equations
L[3] is the list of equations nonlinear in their highest derivatives
L[4] is the list of unclassified equations
L[5] is the list of pivots

as they were at the end of the last iteration.

The file that the intermediate results are stored in can be given a name
different from StandardStorage, if so desired, if, for example, the
results from several runs are needed.  To store the results in a file
named filename, type, when standard_form is invoked

> sf:=standard_form(sys,store=filename);

The results will actually be stored in a file filename.m,
in maple's internal format, and must be be read
into a maple session to view the most updated results:

> L:=readstored(filename);

(note that the .m is not be included).

The command

> standard_form(sys,storeall);

will save the results from each iteration in
a different file, and thus, to have a permanent record of all of them
(if you're not concerned with flooding you file-space).
In particular the results at the end
of the nth iteration will be stored in the file
StandardStorage_n.m.  To view the results of the nth iteration type

> L:=readstored(StandardStorage_n);

It is also possible to store results
to files with user-assigned names:

> standard_form(sys,storeall=filename);

The results at the end of the nth iteration will then be stored in a file
filename_n.m.  To view these files, start a maple session, read in
the standard form package, and type

> L:=readstored(filename_n);

This is useful for:
recovering intermediate results if there is a memory explosion,
and restarting them, (perhaps with the original system and
without some of the unclassified equations), and
also peeking at the current status of the system.

3.4.2 SELECTING NICE EQUATION SUBSETS WITH USEREASY

Including:

usereasy = N     (N a positive integer)

in the option list for standard_form, will force standard_form
to take the first N equations in the input system to standard_form
into the easy list Ez at standard_form's first iteration.

This over rides standard_form's Ez/UnC choice strategy at its first
iteration.

It is most useful when the user has their own premonition of a set
of equations that would help the program out (and places these as the
first N equations in the input system list).

One situaiton in which this can be helpful is when memory explosion
has occurred, and the results of previous run before the explosion
have been saved using one of the storage options.
The user chooses N easy equations, throws out the ugly ones, and
forms a new system, including the original system:

> sf:= standard_form(sys, store);
.....
.....

!!!! memory explosion !!!!

pick N easy eqns from stored eqns, throwing out the ugly ones
Advisable to start a new session if the program went down in flames

> sysnew:= [N easy, op(sys)];

> sf:= standard_form(sysnew, usereasy = N);

3.5  CLASSIFICATION PROBLEMS - SYSTEMS OF PDE WITH UNSPECIFIED
FUNCTIONS OR PARAMETERS

Consider a system of PDEs with dependent variables
V(1),...,V(n) which are functions of independent variables
x1,...,xm.  In addition suppose that the system involves
certain parameters or functions of unspecified form ("classification
parameters or functions").  In general the properties and
number of solutions of the system will depend on the values and
forms of these parameters and functions.   We will call the
"classification problem for the system" the problem of
determining the standard form of the system for every functional
form of the classification functions.  In general a "classification
tree" results, each branch of which is characterised by a certain
number of differential equalities and inequalities amongst the
classification functions and parameters.

3.5.1  LOW-TECH APPROACH TO CLASSIFICATION

Problems where there appear only parameters or say only one
unspecified function can often be handled by the standard form
procedure without embellishment.

Already we have classified the system of algebraic equations
sysA, given in Example A of section 2.2.

The resulting classification tree was:

sysA
|
k=0            |        k#0
-----------------------------------------------
|                                               |
|                                               |
sfA                                             sfA1
idA                                             idA1

where the classification parameter was k and
separate branches were picked out by examining the pivots list _pivs.

We also classified the system of PDEs sysC given in Example C of
Section 2.3 with classification parameter r which yielded the
classification tree:

sysC
|
r=0            |        r#0
-----------------------------------------------
|                                               |
|                                               |
sfC2                                      inconsistent
idC2                                       no solutions

Example A.

We will classify the system of PDEs:

g(t) Ux - g'(t) U = 0 ,     Ut = 0

where U = U(x,t), and g(t) is the classification function.

Then in our notation:

> sysA:= [ g(x2)*V(1,1) - diff(g(x2),x2)*V(1) = 0, V(1,2) = 0];

Case 1:  No constraints on the classification function g(t)
- the generic case.

> sfA1:= standard_form(sysA);

yields

sfA1 := [ V(1) = 0]

and initial data

> idA1:= initial_data(sfA1);

idA1:= [[],[]]

The pivots are

> _pivs;

/   2        \
|  d         |         /  d        \2
|------ g(x2)| g(x2) - |----- g(x2)|
|    2       |         \ dx2       /
\ dx2        /
[-1 = 0, g(x2) = 0, ------------------------------------- = 0]
2
g(x2)

and these can be given a more transparent form by using the sfprogs
procedure shortenpiv which removes repeated factors from the pivots:

> spivsA1:= shortenpiv(_pivs);

/   2        \
|  d         |         /  d        \2
spivsA1 := [g(x2), |------ g(x2)| g(x2) - |----- g(x2)| ]
|    2       |         \ dx2       /
\ dx2        /

So g(x2) = 0 and

2
g''(x2) = g'(x2) /g(x2)

are special cases.

Case 2:  g(x2) = 0

> g:= proc(x2) 0 end;

> sfA2:= standard_form(sysA);

yields

sfA2 := [V(1,2) = 0]

with initial data

> idA2 := initial_data(sfA2);

idA2 := [[], [V(1) = f1(x1)]]

Case 3: g''(x2) = g'(x2)^2/g(x2)

To invoke the above simplification rule we can use Maple's inbuilt
capability for defining how diff acts on a function.
We employ Maple's diff/ capability in two steps:

> diff/g:= proc(x2,x2) h(x2) end;         # i.e. g'(x2) = h(x2)
> diff/h:= proc(x2,x2) h(x2)^2/g(x2) end; # i.e. h'(x2) = h(x2)^2 /g(x2)
> sfA3:= standard_form(sysA);

yields

sfA3 :=  [ V(1,1) = h(x2) V(1)/g(x2),  V(1,2) = 0 ]

with initial data

> idA3:= initial_data(sfA3);

idA3 := [[V(1)=a1], []]

and (shortened) pivots

> spivsA3:= shortenpiv(_pivs);

spivsA3:= [g(x2)]

This is just Case 2 again, so there are no further branches.

3.5.2  HI-TECH APPROACH TO CLASSIFICATION - the class_eqns option

General classification problems involve the classification of
systems of PDEs in dependent variables V(1), ... , V(n) and
independent variables x1, ... , xm and classification functions
g1, ... , gs.

In the hi-tech classification approach the classification functions
g1, ... , gs are regarded as additional dependent variables
("classification variables") V(n+1), ... , V(n+s)
which satisfy an consistent system of PDEs called the classifying equations.

It is important to note that:

THE CLASSIFYING EQUATIONS MUST BE IN STANDARD FORM.

THE DEPENDENCIES OF THE CLASSIFYING VARIABLES MUST BE DECLARED.

THE CLASSIFICATION VARIABLES MUST BE NUMBERED HIGHER THAN THE OTHER
DEPENDENT VARIABLES

In particular the classifying equations are entered as an option:

> standard_form(sys, class_eqns = class)

Here class is the list of the classifying equations in
standard form and the dependencies of each classifying variable V(r)
declared in the form of the list:

[ V(r), k,l,...]

meaning that V(r) depends on xk, xl, etc.

Example A.

Consider again the system of PDEs dealt with by the low-tech approach
in Example A of Section 3.5.1:

g(t) Ux - g'(t) U = 0 ,     Ut = 0

where U = U(x,t), and g(t) is the classification function.
In the hi-tech approach we set U=V(1) and define the classification
variable V(2)=g(t) (the classification variables must be numbered
higher than the nonclassification variables).  So in Maple

> sysA:= [ V(2)*V(1,1) - V(2,2)*V(1) = 0, V(1,2) = 0 ]:

Case 1:  No constraints on the classification variable V(2)
- the generic case

Then the class equations are entered as a list:

> classA1:= [[V(2),2]];

At the moment there are no differential constraints on V(2), so the
only information we enter is the dependency info that V(2) only
depends on x2.  A valid alternative here would have been
classA2:=[[V(2),2],V(2,1)=0], but not classA2:=[V(2,1)=0] since standard_form
uses the dependency declaration to determine which are the classification
variables.  Reduction of the system to standard form gives:

>  sfA1:= standard_form(sysA, class_eqns = classA1);

The system parameters are as follows:

There are 2 independent variables:

[x1, x2]

There are 1 dependent variables:

[V(1)]

There are 1 classification variables:

[V(2)]

sfA1 := [V(1) = 0]

which agrees with Case 1 of Example A in 3.5.1.  Notice that, standard_form
has reassured us by alerting us that V(2) has been declared
as a classification variable.

To calculate the initial data, we apply initial_data with the option
class_eqns = classA1, alerting initial_data that V(2) is a
classification variable:

> idA1:= initial_data(sfA1, class_eqns=classA1);

idA1 := [[], [V(2) = f1(x2)]]

Note that the idA1 includes the initial data for
the classification variable.

The pivots are

>  _pivs;

2
V(2, 2, 2)   V(2, 2)
[-1 = 0, V(2) = 0, ---------- - -------- = 0]
V(2)            2
V(2)

and in their shortened form:

>  spivsA1:= shortenpiv(_pivs);

2
spivsA1 := [V(2, 2, 2) V(2) - V(2, 2) , V(2)]

So V(2) = 0 and
2
V(2, 2, 2) V(2) - V(2, 2)  = 0

are special cases.

Case 2:  V(2) = 0

We need to give both the dependency info and the class equations in
standard form in the list class_eqns.  V(2) as a constant, does not
depend on anything, and to reflect this [V(2)] is included in class_eqns.
Obviously the standard form of the class equations is V(2) = 0.

Reduction of this case to standard form gives:

>  sfA2:= standard_form(sysA, class_eqns = [[V(2)],V(2)=0]);

sfA2 := [V(1, 2) = 0]

and

>  spivsA2:= shortenpiv(_pivs);

spivsA2 := []

so there are no nontrivial shortened pivots.

The initial data is:

>  idA2:= initial_data(sfA2, class_eqns= [[V(2)], V(2)=0]);

idA2 := [[], [V(1) = f1(x1)]]

To alert initial_data is a classification variable, we have inserted the
class_eqns option in the same way as for standard_form.

2
Case 3:   V(2, 2, 2) V(2) - V(2, 2)  = 0

This case corresponds to the vanishing of the first pivot in spivsA1:

>  class_sys:= [spivsA1[1]=0];

2
class_sys := [V(2, 2, 2) V(2) - V(2, 2)  = 0]

To reduce the classifying system to standard form we need the option
divpiv so that standard_form will divide by the pivot V(2):

> sfclass_sys:= standard_form(class_sys, divpiv);
2
V(2, 2)
sfclass_sys := [V(2, 2, 2) = --------]
V(2)

To give the classifying equations in proper format
we include the dependency info for V(2):

>  classA3:= [[V(2),2], op(sfclass_sys)];
2
V(2, 2)
classA3 := [[V(2), 2], V(2, 2, 2) = --------]
V(2)

Reduction of the original system to standard form for this case gives:

> sfA3:= standard_form(sysA, class_eqns = classA3);

V(2, 2) V(1)
sfA3 := [V(1, 1) = ------------, V(1, 2) = 0]
V(2)

> spivsA3:= shortenpiv(_pivs);

spivsA3 := [V(2)]

so that no further branching occurs.

> idA3:= initial_data(sfA3, class_eqns = classA3);

idA3 := [[V(1) = a1, V(2) = a2, V(2, 2) = a3], []]

Note that the initial data includes the initial data for the
classification equation for V(2).

In summary we have the classification tree:

sysA

g(t) Ux - g'(t) U  = 0

Ut = 0

|
|
g(t)=0       |       g(t)#0
----------------------------------------
|                                        |
|                                        |
|
Ut = 0                                     |
|
U(x,t0) = f1(x)                            |
|
Case 2                                     |
|
|
2       |               2
g''(t)=g'(t) /g(t)  |   g''(t)#g'(t) /g(t)

-------------------------------------------
|                                           |
|                                           |

Ut = 0                                      U = 0

Ux = g'(t)U/g(t)                            Case 1

U(x0,t0) = a1

Case 3

Figure:  Classification tree for the system of PDEs
g(t) Ux - g'(t) U = 0,  Ut = 0 for U=U(x,t) with
classification function g(t).
Each branch gives the pivot conditions, initial data
and standard form of the system.

3.6 INTERFACE TO MAPLE'S DERIVATIVE NOTATION

To convert a system of PDEs written in the repeated index notation
specific to our programs we use the convprogs routine rep2map:

Example A

To change the following system which is in our repeated index notation
to Maple's usual notation for derivative:

> sysA := [V(2,3) = 0, V(2,1) = 0, V(1,3) = 0, x1*x2*V(3,3,3) = 0,
V(2,2)-2*V(1,1) = 0, 2*V(3,1,3)+V(1,2)-V(1,1,1) = 0, V(3,2)-V(3,1,1) = 0];

sysA := [V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, x1 x2 V(3, 3, 3) = 0,

V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 1, 3) + V(1, 2) - V(1, 1, 1) = 0,

V(3, 2) - V(3, 1, 1) = 0]

we read the conversion package (sfprogs is not needed):

To use the conversion function rep2map (repeated index to maple)
which takes as its first argument the system of PDEs written as
a list of equations in our repeated-index notation, and as second
argument the number of independent variables (=3 here):

> msys:= rep2map(sysA, 3);

which yields:

msys := [D[3](g)(x1, x2, x3) = 0, D[1](g)(x1, x2, x3) = 0,

D[3](f)(x1, x2, x3) = 0, x1 x2 D[3, 3](h)(x1, x2, x3) = 0,

D[2](g)(x1, x2, x3) - 2 D[1](f)(x1, x2, x3) = 0,

2 D[1, 3](h)(x1, x2, x3) + D[2](f)(x1, x2, x3) - D[1, 1](f)(x1, x2, x3) = 0,

D[2](h)(x1, x2, x3) - D[1, 1](h)(x1, x2, x3) = 0]

This function is useful for interfacing our package sfprogs with
other packages and procedures e.g.  Maple's dsolve command.

A procedure, map2rep is also available which does the opposite
conversion.  If you have a system written in maple's standard
derivative notation, in some dependent and independent variables,
with arbitrary names, map2rep will convert it to a system in
repeated index notation.

The syntax for the map2rep procedure is

> map2rep(expr, deps, indeps)

where expr is the expression to be converted, and deps and indeps are
lists of the names used in expr for the dependent and independent
variables, in the order that they are to be labelled in the repeated
index notation.

In the situation where the independent variables are already
written as xi, i=1..n, just include a list of independents
consisting of the xi.  The independent variables appearing
explicitly in the expression will then be unaffected by the
conversion procedure.

As an example, we do the reverse of the above example.  To
show the most general situation, we enter the system with
different names for the independents: we make the substitutions
x1=x,x2=y,x3=z.  The system, in maple notation, is then

sysC := [D[3](g)(x, y, z) = 0, D[1](g)(x, y, z) = 0, D[3](f)(x, y, z) = 0,

x y D[3, 3](h)(x, y, z) = 0, D[2](g)(x, y, z) - 2 D[1](f)(x, y, z) = 0,

2 D[1, 3](h)(x, y, z) + D[2](f)(x, y, z) - D[1, 1](f)(x, y, z) = 0,

D[2](h)(x, y, z) - D[1, 1](h)(x, y, z) = 0]

To convert it to repeated index notation, we use the program
map2rep as follows:
>  map2rep(sysC,[f,g,h],[x,y,z]);

[V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, x1 x2 V(3, 3, 3) = 0,

V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 1, 3) + V(1, 2) - V(1, 1, 1) = 0,

V(3, 2) - V(3, 1, 1) = 0]

and this is the system we started with originally.

3.7 INTERFACE TO MANSFIELD'S NOTATION FOR HER DIFFGROB PACKAGE

As we have already mentioned, the program standard_form is basically a
linear solver, but with some capabilities for nonlinear systems of PDEs.
An approach and a Maple program for simplifying polynomially nonlinear
systems of PDEs has been given by
Mansfield [Mans91].  We illustrate by example how to convert back and forth
between Mansfield's notation (recently she changed her notation, but we
have not yet changed our conversion program to work for this changed
notation) and our own:

Example A

To change the following system which is in our repeated index notation
to the notation used by Mansfield's program:

> sysA := [V(2,3) = 0, V(2,1) = 0, V(1,3) = 0, x1*x2*V(3,3,3) = 0,
V(2,2)-2*V(1,1) = 0, 2*V(3,3,1)+V(1,2)-V(1,1,1) = 0, V(3,2)-V(3,1,1) = 0];

sysA := [V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, x1 x2 V(3, 3, 3) = 0,

V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 3, 1) + V(1, 2) - V(1, 1, 1) = 0,

V(3, 2) - V(3, 1, 1) = 0]

To use the conversion function rep2man (repeated index to mansfield)
we first have to declare the global variables vars: (the independent
variables) and ukns (the names of the dependent variables).

> vars:= [x1,x2,x3]:

> ukns:= [v1,v2,v3]:

The system in Mansfield's notation is then:

> msysA:= rep2man(sysA);

yielding

msysA := [p(v2, [x1, x2, x3], [0, 0, 1]), p(v2, [x1, x2, x3], [1, 0, 0]),

p(v1, [x1, x2, x3], [0, 0, 1]), x1 x2 p(v3, [x1, x2, x3], [0, 0, 2]),

p(v2, [x1, x2, x3], [0, 1, 0]) - 2 p(v1, [x1, x2, x3], [1, 0, 0]),

2 p(v3, [x1, x2, x3], [1, 0, 1]) + p(v1, [x1, x2, x3], [0, 1, 0])

- p(v1, [x1, x2, x3], [2, 0, 0]),

p(v3, [x1, x2, x3], [0, 1, 0]) - p(v3, [x1, x2, x3], [2, 0, 0])]

To convert a system written in Mansfield's notation to repeated index
notation we use the function man2rep (again both vars - the names of the
independent variables in the Mansfield system (these may be different to
x1,x2,...) and ukns (the names of the dependent variables):

> vars:= [x1,x2,x3]:

> ukns:= [v1,v2,v3]:

> rsysA := man2rep(msysA);

yielding

rsysA := [V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, x1 x2 V(3, 3, 3) = 0,

V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 1, 3) + V(1, 2) - V(1, 1, 1) = 0,

V(3, 2) - V(3, 1, 1) = 0]

3.8 THE DIRECTED GRAPH OF A PDE SYSTEM

This procedure creates the directed graph of a PDE system
as described in [Reid91c].  This procedure is still
preliminary and has not been thoroughly checked.  But have
fun anyway!

Example:  xmaple (Maple V release 2) is needed for this application

# from sfprogs

> with(networks);

> read heat.sys:  # the determining system for the Heat equation

> sf:= standard_form(sys, janet_order=[[0,0,1,0,0,0]]);
# the standard form and the graph will depend on
# the ordering which is used

> msf:= make_vset(sf);
#returns the vertices and edges of the graph

> G:= graph(msf):  # makes the data structure for the graph

> show(G);         # shows the graph as a data structure

> draw(G);         # attempts to draw the graph is directed -
# but xmaple (Maple V release 2) does
# not yet draw arrows on the edges

4. APPLICATION OF THE PACKAGE TO SYMMETRY ANALYSIS OF PDES

A symmetry of a system of PDEs is a transformation which leaves
invariant the family of solutions of the system.  Symmetries have
proved to be a useful tool for deriving analytical and geometrical
properties of systems of PDE:  invariant solutions, mappings between
equations, connections with conservation laws (Noether's Theorem), and
the relationship of higher order symmetries (Lie-Backlund
symmetries) with the inverse-scattering method.

Lie point symmetries are the simplest type of symmetries of systems
of PDEs.  Given a system of PDEs with independent variables
x=(x1,x2,...,xp) and dependent variables u=(u1, u2, ... , uq) a
1 parameter (a)-Lie symmetry group is a group of transformations
of the form:

x' = X(x,u; a)

u' = U(x,u; a)

i.e. the transformation parameterised by a, maps a point (x,u)
to a point (x', u') in the space of independent and dependent variables
of the system of PDEs.

Instead of dealing with the group transformations in their (nonlinear)
form above Lie dealt with the transformations in their (linearized)
infinitesimal form:

2
x' = X(x,u; a) = x + V(x,u) a + O(a )

2
u' = U(x,u; a) = u + W(x,u) a + O(a )

with associated infinitesimal generator

L = V(1)d/dx1 + V(2)d/dx2 + ... + V(p)d/dxp + W(1)d/du1 + ... + W(q)d/duq

in the infinitesimals V(1),...,V(p),W(1),...,W(q).

In order to input such infinitesimals into our programs we'll replace
uk, W(k) by x.(p+k), V(p+k) (x.1 = x1, and u1, W(1) is replaced by
x.(p+1), V(p+1) etc. ).  In this notation the infinitesimal operator is

L = V(1) d/dx1 + ... + V(r) d/dxr where r = p+q.

According to a fundamental theorem of Lie the infinitesimals
V(1),...,V(r) uniquely characterise the
group.  Lie provided a simple algorithm [BlK89], [Olv86], [Ovs82] for
deriving the equations which must be satisfied by the infinitesimals.
Today there are many effective symbolic programs for deriving these
so-called "defining" or "determining equations" [Here93].

The determining equations of a given system of PDEs form a linear
homogeneous system of PDEs (regardless of whether the original
system was linear or not).  This linearization of the group
problem was a key to the success of Lie's methods.
In addition such determining systems are usually very overdetermined.

Another fundamental result of Lie is that if M and N are two
infinitesimal symmetry generators of a PDE system then so too is the
commutator of M and N:

[ M, N ] = M N - N M

So these infinitesimal generators form a Lie algebra with respect to [,].
If this Lie algebra is finite-dimensional then according to Lie it is
characterised by its "structure constants" [Ovs82, Olv86, BlK89]

k
[L(i), L(j)] = C   L(k)     (sum on k)
ij

where L(i) is the infinitesimal basis generator corresponding to
the i-th basis solution of the determining equations.
(i.e. L = a1 L(1) + a2 L(2) + ... + aK L(K) where a1, a2, ... ,aK
are the parameters in the general solution of the determining system).

The structure of this algebra provides important coordinate-independent
(geometric) information about the original system of PDEs:  information
about reductions to ODEs, mappings to other more tractable PDEs etc.

4.1 CALCULATION OF THE STRUCTURE OF LIE SYMMETRY ALGEBRAS
OF PDES FROM THEIR INFINITESIMAL DETERMINING SYSTEMS

The conventional approach to infinitesimal symmetries has been to
try to integrate the determining equations, explicitly determining the
basis generators L(i), and then substituting the results into the
commutation relations above to determine the structure constants.

Since this process is based on the heuristic process of integration
it is not always guaranteed to succeed.

Our approach to finding the structure of Lie symmetry algebras of
PDEs is completely algorithmic:

(1) The determining system for the infinitesimals is either
produced by hand or by one of the many computer packages
automating this process [Here93, Champ91].
This system is converted into the xi,V(j) notation used
by our sfprogs package (see 4.2,4.3,4.4).
(2) Algorithm standard_form is used to reduce the determining system
to standard form.
(3) Application of initial_data to the standard form produced in Step (2)
algorithmically determines the dimension of the algebra
of infinitesimal symmetries.
(4) Algorithm fin_com_table uses the standard form from Step (2)
and its initial data from Step (3) to determine the structure of the
of the finite-dimensional subalgebra of infinitesimal symmetries
corresponding to the parameters a1,a2,... in the finite initial data.
This algorithm depends only on differentiation and substitution.
(5) taylor series approximations of the infinitesimals or basis generators
are found to any chosen finite degree using algorithms taylor_exp,
or finite_generators respectively.

Notes
1.  Commutators of generators belonging to infinite-dimensional
Lie symmetry algebras can also be computed using the procedure
lie_symm (instead of fin_com_table).  However the output of lie_symm is
complicated and currently does not have a good structural interpretation.
2.  Classes of PDEs with variable coefficients can be treated
using the class_eqns option (see Section 4.4.2).
3.  In many applications the explicit forms of the infinitesimals are
required (e.g. in the calculation of invariant solutions) so our algorithms
do not make redundant the packages which attempt to explicitly solve
determining equations.  The role of our algorithms is instead a
complementary one and we are interested in developing interfaces
with such packages.

4.1.1 SYMMETRIES OF THE HEAT EQUATION:  Ut = Uxx

We will not violate a longstanding tradition so we first
consider the heat equation: Ut - Uxx = 0 where U=U(x,t).

Step (1):
In Section 4.2.1 it is shown how to produce the determining
system for the heat equation using the macsyma package symgrp.max.
The converted system has been stored in the file called heat.sys.

We next read the file heat.sys containing the infinitesimal determining
system for the Heat equation

> read heat.sys;

sys := [V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, V(3, 3, 3) = 0,

V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 3, 1) + V(1, 2) - V(1, 1, 1) = 0,

V(3, 2) - V(3, 1, 1) = 0]

The infinitesimal generator of the algebra is given
by L = V(1) d/dx1 + V(2) d/dx2 + V(3) d/dx3, where
x1 = x (space), x2 = t (time) and x3 = U (temperature),
with corresponding infinitesimals V(1), V(2) and V(3).

Step (2):
We reduce the system to standard form:

> sf:= standard_form(sys);

sf := [V(3, 2, 2, 3) = 0, V(3, 1, 1) = V(3, 2), V(3, 1, 3) = - 1/2 V(1, 2),

V(1, 2, 2) = 0, V(2, 2, 2) = - 4 V(3, 2, 3), V(3, 3, 3) = 0,

V(1, 1) = 1/2 V(2, 2), V(2, 1) = 0, V(1, 3) = 0, V(2, 3) = 0]

(For fun also try:
> sffun:= standard_form(sys, janet_order=[[0,1,0,0,0,0]]);
> idfun:= initial_data(sffun);      )

Step (3):
We calculate the initial data of the system:

> id:= initial_data(sf);

id := [
[
V(1)=a1, V(1, 2)=a2, V(2)=a3, V(2, 2)=a4, V(3, 3)=a5, V(3, 2, 3)=a6
],
[V(3)=f1(x2), V(3, 1)=f2(x2)]]

See Section 3.2 for interpretation of the initial data.
So the dimension of the group of Lie symmetries of the system
is infinite-dimensional.  It is determined by 6 parameters
and 2 arbitrary functions of 1 variable.  Note that we did not
need to integrate the determining system to uncover this information.

Step (4):
To determine the commutation relations amongst the generators of
the finite (6-dimensional) component of the algebra we use the
procedure fin_com_table:

>  ct:= fin_com_table(sf, id);

[ 0  - 1/2 L(5)     0     1/2 L(1)   0      - 2 L(2)    ]
[                                                       ]
[ *       0      - L(1)  - 1/2 L(2)  0         0        ]
[                                                       ]
[ *       *         0       L(3)     0  - 4 L(4) + L(5) ]
ct := [                                                       ]
[ *       *         *         0      0        L(6)      ]
[                                                       ]
[ *       *         *         *      0         0        ]
[                                                       ]
[ *       *         *         *      *         0        ]

The i-j entry of this table corresponds to [L(i), L(j)]
(e.g. [L(1), L(2)] = -1/2 L(5)).  To look at ct later use op(ct) or print(ct).

Step (5):
To obtain the taylor series approximations of the infinitesimals
to degree 3 about (X1,X2,X3) in powers of h1,h2,h3:

>  ts:= taylor_exp(sf, id, 3);

ts := [V(1) = a1 + 1/2 a4 h1 + a2 h2 - 2 a6 h1 h2,

2
V(2) = a3 + a4 h2 - 2 a6 h2 ,

V(3) =
/  d         \                  /  d         \   2
f1(X2) + f2(X2) h1 + |----- f1(X2)| h2 + a5 h3 + 1/2 |----- f1(X2)| h1
\ dX2        /                  \ dX2        /

/   2         \
/  d         \                            |  d          |   2
+ |----- f2(X2)| h1 h2 - 1/2 a2 h1 h3 + 1/2 |------ f1(X2)| h2  + a6 h2 h3
\ dX2        /                            |    2        |
\ dX2         /

/   2         \
/  d         \   3       |  d          |   2               2
+ 1/6 |----- f2(X2)| h1  + 1/2 |------ f1(X2)| h1  h2 + 1/2 a6 h1  h3
\ dX2        /           |    2        |
\ dX2         /

/   2         \              /   3         \
|  d          |      2       |  d          |   3
+ 1/2 |------ f2(X2)| h1 h2  + 1/6 |------ f1(X2)| h2
|    2        |              |    3        |
\ dX2         /              \ dX2         /

]

By linearity we can express the general infinitesimal generator in terms
of basis generators L(1),...,L(6) corresponding to the initial data
a1,a2,...,a6 respectively and two basis generators M1(f1(x2)), M2(f2(x2))
corresponding to the arbitrary functions f1,f2 in the initial data.
That is: L = a1 L(1) + a1 L(2) + ... + a6 L(6) + M1(f1(x2)) + M2(f2(x2))

In particular we can get taylor approximations of L(1),...,L(6)
by making substitutions in the taylor series above:
For example: L(1) = T[1]*Dx(1) + T[2]*Dx(2) + T[3]*Dx(3)
where T[i] = coefficient a.i in ts[i] and Dx(1) = d/dx1 etc.
We've provided a procedure finite_generators, which does this
automatically.  For example to obtain the taylor
approximations to degree 4 of the finite (subgroup) basis generators:

>  fgs:= finite_generators(sf, id, 4);

fgs := [L(1) = Dx(1), L(2) = h2 Dx(1) - 1/2 h3 h1 Dx(3), L(3) = Dx(2),

L(4) = 1/2 h1 Dx(1) + h2 Dx(2), L(5) = h3 Dx(3),

2                           2
L(6) = - 2 h1 h2 Dx(1) - 2 h2  Dx(2) + (h3 h2 + 1/2 h3 h1 ) Dx(3)]

where Dx(1) represents d/dx1 etc.  In our original notation
expanding about (X1,X2,X3)=(0,0,0) these generators are:

L(1) = d/dx + O(5),   L(2) = t d/dx - 1/2 x U d/dU + O(5),
L(3) = d/dt + O(5),   L(4) = 1/2 x d/dx + t d/dt + O(5)
L(5) = U d/dU + O(5), L(6) = -2xt d/dx - 2 t^2 d/dt +(tU + x^2 U/2)d/dU + O(5)

Similarly, although we have not automated it yet, we can obtain
taylor approximations of the generators M(f1(x2)) and M2(f2(x2))
corresponding to the infinite initial data.
For the moment you can do this yourself for M1(f1(x2)) by
assigning f2:= proc(z) 0 end: and substituting a1=0,a2=0,...,a6=0 into ts.
And similarly for M2(f2(x2)).

As is well known the taylor approximations of the 6-dimensional finite
symmetry subalgebra of the Heat equation are at most cubic in x,t,U, so the
remainder terms O(5) in the above generators are exact.
There are no 4-th order terms in the above expansions and this motivates
the general question:  If there are no N-th order terms in the taylor
series approximations of the generators does this mean that they are exact
to order N-1 ? (Answer - Yes:  see Section 5.1).

We emphasise, however, that the termination of the Taylor series
had nothing to do with the effectiveness of the algorithm fin_com_table
for calculating the structure of the finite dimensional subalgebra above.
In particular, fin_com_table does not even use Taylor series in
the construction of structure constants (instead it uses the approach
of [Reid92a] which makes the Taylor series approach of [Reid91b] obsolete).

Finally the command

> lie_symm(sf, id);

will display the entire 8 X 8 commutator table, between the operators
L(1),L(2),L(3),L(4),L(5),L(6),M1(f1(x2)),M2(f2(x2)), in a somewhat
obscure vector-format (how about that for great documentation!).
In particular the algorithms can calculate the commutators for infinite-dimensional
Lie symmetry groups in terms of their initial data [Reid90, Reid91b, Reid92b]
(again no knowledge of the solutions is required).  However we have not
as yet developed a good structural (coordinate-independent) interpretation of this table.

Despite its name lie_symm
is not a part of or an application of the differential forms package liesymm.

4.1.2 A GROUP CLASSIFICATION PROBLEM:  Uxx = Utt + H(U) Ut

The group classification problem for the family of nonlinear
telegraph equations

Uxx = Utt + H(U) Ut where U=U(x,t)    (*)

is to determine the size and structure of the group of Lie
symmetries of (*) for every functional form of H(U) (see [Reid91b] for
complete results).  To illustrate the
method We will follow a couple of branches of the classification.

Step (1):
In Section 4.2.2 it is shown how to produce the determining
system for (*) by using the macsyma package symgrp.max.
The converted system has been stored in the file called telegraph.sys.

We next read the file telegraph.sys containing the infinitesimal determining
system of (*):

> read telegraph.sys;

sys := [V(1, 3) = 0, V(2, 3) = 0, V(3, 3, 3) = 0, V(2, 1) - V(1, 2) = 0,

V(2, 2) - V(1, 1) = 0, V(3, 2) V(4) - V(3, 2, 2) + V(3, 1, 1) = 0,

V(1, 2) V(4) - 2 V(3, 3, 1) - V(1, 2, 2) + V(1, 1, 1) = 0,

V(3) V(4, 3) - V(2, 2) V(4) + 2 V(1, 1) V(4) + V(2, 2, 2) - V(2, 1, 1)

- 2 V(3, 3, 2) = 0                                              ]

The infinitesimal generator of the algebra is given
by L = V(1) d/dx1 + V(2) d/dx2 + V(3) d/dx3, where
x1 = x (space), x2 = t (time) and x3 = U,
with corresponding infinitesimals V(1), V(2) and V(3).
The classification variable V(4) depends only on x3.

Case 1:  No constraints on V(4).

Step (2):
Since we are treating V(4) as a classification variable (see 3.5.2)
we have to enter its dependency info in the list class_eqns.
Reducing the system to standard form

> sf1:= standard_form(sys, class_eqns=[[V(4),3]]);

sf1 := [V(1, 1) = 0, V(2, 1) = 0, V(1, 2) = 0, V(2, 2) = 0, V(1, 3) = 0,

V(2, 3) = 0, V(3) = 0]

Step (3):
Calculation of the initial data of the system gives

> id1:= initial_data(sf1, class_eqns=[[V(4),3]]);

id1 := [[V(1) = a1, V(2) = a2], [V(4) = f1(x3)]]

where we have entered the class_eqns option to initial_data
since there is a classification variable (see 3.5.2).
The group of Lie symmetries of the system is 2-dimensional.
Since we are dealing with a classification problem we look
at the (shortened) pivots to identify special cases:

> spivs1:= shortenpiv(_pivs);

2
spivs1 := [- 5/2 V(4, 3)  + V(4) V(4, 3, 3),

2                                                  2
- 2 V(4) V(4, 3, 3)  + V(4) V(4, 3, 3, 3) V(4, 3) + V(4, 3, 3) V(4, 3) ,

V(4), V(4, 3)]

So the above standard form and initial data are valid if these pivots are
not identically zero.

Step (4):
The commutator table is given by:

>  ct1:= fin_com_table(sf1, id1, class_eqns=[[V(4),3]]);

[ 0  0 ]
ct1 := [      ]
[ *  0 ]

Note that again we had to enter the class_eqns option to a
procedure since we are dealing with a classification problem.

Step (5):
Neither taylor_exp nor finite_generators have a class_eqns option
yet.  So for the time being you have to manipulate your problems
into a form with a non-classification form to use these functions
and the right initial data.

From the pivots spivs1 we see that a special case is:

2              2
Case 2:  H(U) H'(U) H'''(U) - 2 H(U) H''(U)  + H''(U) H(U)  = 0

that is spivs1[2] = 0.  The class equations must be in solved form

> sfclass2:= standard_form([spivs1[2]=0], divpiv);
2
V(4, 3) V(4, 3, 3)     V(4, 3, 3)
sfclass2 := [V(4, 3, 3, 3) = - ------------------ + 2 -----------]
V(4)              V(4, 3)

where the option divpiv had to be used to force standard_form to
divide by the nonlinear pivot V(4)V(4,3) and solve for V(4,3,3,3).

The class equations for Case 2 are:

> class2:= [[V(4),3], op(sfclass2)];
2
V(4, 3) V(4, 3, 3)     V(4, 3, 3)
class2 := [[V(4), 3], V(4, 3, 3, 3) = - ------------------ + 2 -----------]
V(4)              V(4, 3)

Step (2):
Standard forming the system for Case 2 gives

> sf2 :=  standard_form(sys, class_eqns=class2);
V(3) V(4, 3)
sf2 := [V(1, 1) = - ------------, V(2, 1) = 0, V(3, 1) = 0, V(1, 2) = 0,
V(4)

V(3) V(4, 3)
V(2, 2) = - ------------, V(3, 2) = 0, V(1, 3) = 0, V(2, 3) = 0,
V(4)

V(3) V(4, 3)   V(3) V(4, 3, 3)
V(3, 3) = ------------ - ---------------]
V(4)           V(4, 3)

with shortened pivots

> spivs2:= shortenpiv(_pivs);

spivs2 := [V(4), V(4, 3)]

Step (3):  The initial data is:

> id2:= initial_data(sf2, class_eqns=class2);

id2 := [

[V(1)=a1, V(2)=a2, V(3)=a3, V(4)=a4, V(4, 3)=a5, V(4, 3, 3)=a6], []

]

Note that the initial data a4,a5,a6 for second order ODE satisfied by the
classification variable V(4) was included in id2.
The initial data corresponding to the infinitesimals V(1), V(2), V(3)
is a1, a2, a3 so the symmetry algebra for Case 2 is 3-dimensional.

Step (4):  Calculation of the commutator table gives:

> ct2 :=  fin_com_table(sf2, id2, class_eqns=class2);

[         a5 L(1) ]
[ 0  0  - ------- ]
[            a4   ]
[                 ]
ct2 := [         a5 L(2) ]
[ *  0  - ------- ]
[            a4   ]
[                 ]
[ *  *      0     ]

Notice how initial data for the classification variable V(4) has
entered into the commutator table.

Step (5):
As we've already mentioned taylor_exp and finite_generators don't have
a class_eqns option so they can't be directly applied to this case.
Cunning such as that exhibited below is required to apply taylor_exp:

> ts2:= taylor_exp([op(sf2), op(sfclass2), V(4,1)=0, V(4,2)=0], id2, 2);

a3 a5 h1              a3 a5 h2
ts2 := [V(1) = a1 - --------, V(2) = a2 - --------,
a4                    a4

/a3 a5   a3 a6\                                  2
V(3) = a3 + |----- - -----| h3, V(4) = a4 + a5 h3 + 1/2 a6 h3 ]
\  a4      a5 /

Since there are no second order terms in the taylor expansion of the
infinitesimals V(1), V(2), V(3) we can assume that these expansions have terminated at
first order yielding by inspection of ts2 the exact (see 5.1) symmetries:

L(1) =  Dx(1),  L(2) = Dx(2),

L(3) = - (a5/a4) h1 Dx(1) - (a5/a4) h2 Dx(2) + (a5/a4 - a6/a5) h3 Dx(3),

where a4 = H(U0), a5 = H'(U0) and a6 = H''(U0).

4.2 INTERFACE TO THE MACSYMA PACKAGE SYMGRP.MAX

of Champagne-Hereman-Winternitz (available by ftp)

In [Cha90, Here93] a macsyma program symgrp.max is described which
produces the determining
equations for the infinitesimal symmetries of systems of PDEs.
A feedback mechanism is also described for the attempted solution
of such systems.  We illustrate through examples how the
determining systems produced by symgrp.max can be converted into
a form suitable for input into our maple program standard_form.

4.2.1 PRODUCING THE DETERMINING SYSTEM OF Ut = Uxx

Production of the infinitesimal determining system for the
Heat Equation Ut - Uxx = 0 where U=U(x,t).

To produce the determining system we batch the macsyma-syntax files
heat.dat and maxconv in macsyma, both of which reside
in the same directory as sfprogs:

(c1) batch("heat.dat");
(c2) batch("maxconv");

heat.dat is typical of the macsyma batch files needed to produce
infinitesimal determining systems of PDEs - in this case the
heat equation (have a peek at it).
You'll notice that heat.dat requires the macsyma program
symgrp.max which is available by ftp.
We have provided a macsyma program, "maxconv",  which takes the
determining system produced by the typical macsyma dat files
and writes them to a file called lode in your current directory
which is in a form more suitable for conversion into maple syntax.

After executing the above commands in macsyma, a file
called lode should appear in your directory which has contents:

Listing of lode:

[2,
['DIFF(ETA(2),U(1),1),'DIFF(ETA(2),X(1),1),
'DIFF(ETA(1),U(1),1),'DIFF(PHI(1),U(1),2),
'DIFF(ETA(2),X(2),1)-2*'DIFF(ETA(1),X(1),1),
2*'DIFF(PHI(1),U(1),1,X(1),1)+'DIFF(ETA(1),X(2),1)
-'DIFF(ETA(1),X(1),2),
'DIFF(PHI(1),X(2),1)-'DIFF(PHI(1),X(1),2)]];

Before lode can be read into maple
it must be stripped of ' using your editor
(maybe this could be avoided by extracting lode in a
smarter way from the macsyma program, but no one has figured out
how to do this yet).

After stripping, lode is:

[2,
[DIFF(ETA(2),U(1),1),DIFF(ETA(2),X(1),1),
DIFF(ETA(1),U(1),1),DIFF(PHI(1),U(1),2),
DIFF(ETA(2),X(2),1)-2*DIFF(ETA(1),X(1),1),
2*DIFF(PHI(1),U(1),1,X(1),1)+DIFF(ETA(1),X(2),1)
-DIFF(ETA(1),X(1),2),
DIFF(PHI(1),X(2),1)-DIFF(PHI(1),X(1),2)]];

ETA(1),ETA(2) correspond to the infinitesimals for the
independent variables X(1),X(2) (i.e. space x and time t),
and PHI(1) corresponds to the infinitesimal for the
dependent variable U(1) (the temperature U(x,t)).
To produce the determining equations in a form
suitable for simplification by standard_form
we start a maple session and read our maple package convprogs:

Then we read the determining system and give it a name:

> msys:= ":

The conversion of msys to repeated index notation is
accomplished by applying the function max2rep:

> sys:= max2rep(msys);

obtaining

sys := [V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, V(3, 3, 3) = 0,

V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 3, 1) + V(1, 2) - V(1, 1, 1) = 0,

V(3, 2) - V(3, 1, 1) = 0]

Notice the correspondence:  ETA(1)=V(1),ETA(2)=V(2),PHI(1)=V(3)
and X(1)=x1, X(2)=x2, U(1)=x3.

You can either save sys to a file e.g.

> save sys, heat.sys;

or read in the sfprogs package, and apply the functions in that package
(see Section 4.1).

Note:  The message Warning - DIFF had an argument
which was neither ETA nor PHI'
indicates that your system contains unspecified functions
other than the infinitesimals ETA and PHI.
These will require individual attention for conversion
(i.e. you're probably dealing with a classification problem
see the following example).

4.2.2 PRODUCING THE DETERMINING SYSTEM OF Uxx = Utt + H(U) Ut

Production of the infinitesimal determining system for the
class of nonlinear telegraph equations Uxx = Utt + H(U) Ut where U=U(x,t)
and H is a classification function.

We make up a dat file which is very similar to
the heat.dat file of Section 4.2.1 (we have called this file
telegraph.dat and its available for inspection in the same
directory as heat.dat).

When the PDE involves a classification function (e.g. H here)
then that information is entered in the macsyma dat file
by using the macsyma depends command
(e.g. here the line
depends([H], [U[1]]);

Following the same procedure as in Example A, we execute the
commands:

(c1) batch("telegraph.dat");
(c2) batch("maxconv");

The resulting file, lode which has been written to your current
directory and which contains the determining system
in macsyma format (after removing  using your editor), is:

listing of lode:

[2,
[DIFF(ETA(1),U(1),1),DIFF(ETA(2),U(1),1),DIFF(PHI(1),U(1),2),
DIFF(ETA(2),X(1),1)-DIFF(ETA(1),X(2),1),
DIFF(ETA(2),X(2),1)-DIFF(ETA(1),X(1),1),
DIFF(PHI(1),X(2),1)*h-DIFF(PHI(1),X(2),2)+DIFF(PHI(1),X(1),2),
DIFF(ETA(1),X(2),1)*h-2*DIFF(PHI(1),U(1),1,X(1),1)-DIFF(ETA(1),X(2),2)
+DIFF(ETA(1),X(1),2),
PHI(1)*DIFF(h,U(1),1)-DIFF(ETA(2),X(2),1)*h+2*DIFF(ETA(1),X(1),1)*h
+DIFF(ETA(2),X(2),2)-DIFF(ETA(2),X(1),2)
-2*DIFF(PHI(1),U(1),1,X(2),1)]];

To convert this into a format that maple and standard_form can handle,
we start a maple session:

We read in the file that was created by our macsyma session

and give the determining system a name:

> msys:= ":

To use the hi-tech classification where h is regarded as a
classification variable we substitute h with PHI(2):

> msys:= subs(h=PHI(2), msys);

msys := [2,

[DIFF(ETA(1), U(1), 1), DIFF(ETA(2), U(1), 1), DIFF(PHI(1), U(1), 2),

DIFF(ETA(2), X(1), 1) - DIFF(ETA(1), X(2), 1),

DIFF(ETA(2), X(2), 1) - DIFF(ETA(1), X(1), 1),

DIFF(PHI(1), X(2), 1) PHI(2) - DIFF(PHI(1), X(2), 2) + DIFF(PHI(1), X(1), 2),

DIFF(ETA(1), X(2), 1) PHI(2) - 2 DIFF(PHI(1), U(1), 1, X(1), 1)

- DIFF(ETA(1), X(2), 2) + DIFF(ETA(1), X(1), 2),

PHI(1) DIFF(PHI(2), U(1), 1) - DIFF(ETA(2), X(2), 1) PHI(2)

+ 2 DIFF(ETA(1), X(1), 1) PHI(2) + DIFF(ETA(2), X(2), 2)

- DIFF(ETA(2), X(1), 2) - 2 DIFF(PHI(1), U(1), 1, X(2), 1)]

]

This ensures that when the conversion program is applied to msys:

> sys:= max2rep(msys);

sys := [V(1, 3) = 0, V(2, 3) = 0, V(3, 3, 3) = 0, V(2, 1) - V(1, 2) = 0,

V(2, 2) - V(1, 1) = 0, V(3, 2) V(4) - V(3, 2, 2) + V(3, 1, 1) = 0,

V(1, 2) V(4) - 2 V(3, 3, 1) - V(1, 2, 2) + V(1, 1, 1) = 0,

V(3) V(4, 3) - V(2, 2) V(4) + 2 V(1, 1) V(4) + V(2, 2, 2) - V(2, 1, 1)

- 2 V(3, 3, 2) = 0                                              ]

h and its derivatives will be replaced by V(4) and the appropriate
derivatives of V(4).  At this stage you can either save sys to a file e.g.

> save sys, telegraph.sys;

or read in sfprogs and continue (see Section 4.1).

4.3 INTERFACE TO THE MAPLE PACKAGE SYMMETRY

Mark Hickman [Hick93] has written a package for producing infinitesimal determining
equations and solving them interactively.  Using the map2rep procedure,
explained in section 3.6, it is possible to convert the equations produced
by this package into repeated index notation, interactively, to use
the standard form programs on them.  Currently, Hickman's package only
is only available for release 2 of Maple 5
His package exploits the postscript interface so that greek letters etc.
can be displayed and calculations can be done in much the way one would do them by hand.

As an example we use Hickman's package symmetry'' to generate the determining
system for the class of nonlinear telegraph systems:

Vt = Ux ,

Vx = C(U) Ut  +  B(U)

with dependent variables U = U(x,t), V = V(x,t), and classification
functions C(U) and B(U).  First we read in Hickman's symmetry package.

To use Hickman's package, the dependencies of the arbitrary functions
B and C on u must be entered by using his domain command

> domain([u],B):

> domain([u],C):

Hickman's package has the procedure setup for setting up the
maple environment to use his symmetry tools:

> setup();

Welcome to Symmetry

Please remember to end each input with a ;

The independent variables are x,t;
The dependent variables are u,v;

The coefficients of the independent variables are
> X,T;

The coefficients of the dependent variables are
> Y,Z;

The DEs are (use the Diff operator) Diff(v,t)=Diff(u,x),Diff(v,x)=C*Diff(u,t)+B

The independent variables are

[x, t]

The dependent variables are

[u, v]

The symmetry generator is

/  d   \     /  d   \     /  d   \     /  d   \
P -> X |---- P| + T |---- P| + Y |---- P| + Z |---- P|
\ dx   /     \ dt   /     \ du   /     \ dv   /

The differential equations are

d        d       d        /  d   \
---- v = ---- u, ---- v = C |---- u| + B
dt       dx      dx        \ dt   /

Are these correct? y or n (do not use a ;)

y
Note that the jet form of the equations are given by the variable eqn

> eqn;

[vt = ux, vx = C ut + B]

To find the determining system, we use the procedure geneqn:

> eq:=op(geneqn(eqn));

/  d   \   /  d   \      /  d   \
eq := - |---- Y| + |---- X| ux + |---- T| ut
\ dx   /   \ dx   /      \ dx   /

//  d   \   /  d   \      /  d   \   \
- ux ||---- Y| - |---- X| ux - |---- T| ut|
\\ du   /   \ du   /      \ du   /   /

//  d   \   /  d   \      /  d   \   \   /  d   \
- (C ut + B) ||---- Y| - |---- X| ux - |---- T| ut| + |---- Z|
\\ dv   /   \ dv   /      \ dv   /   /   \ dt   /

/  d   \              /  d   \
- |---- X| (C ut + B) - |---- T| ux
\ dt   /              \ dt   /

//  d   \   /  d   \              /  d   \   \
+ ut ||---- Z| - |---- X| (C ut + B) - |---- T| ux|
\\ du   /   \ du   /              \ du   /   /

//  d   \   /  d   \              /  d   \   \
+ ux ||---- Z| - |---- X| (C ut + B) - |---- T| ux|,
\\ dv   /   \ dv   /              \ dv   /   /

/  /  d   \      /  d   \\   /  d   \   /  d   \              /  d   \
Y |- |---- C| ut - |---- B|| + |---- Z| - |---- X| (C ut + B) - |---- T| ux
\  \ du   /      \ du   //   \ dx   /   \ dx   /              \ dx   /

//  d   \   /  d   \              /  d   \   \
+ ux ||---- Z| - |---- X| (C ut + B) - |---- T| ux|
\\ du   /   \ du   /              \ du   /   /

//  d   \   /  d   \              /  d   \   \    /  d   \
+ (C ut + B) ||---- Z| - |---- X| (C ut + B) - |---- T| ux| - (|---- Y|
\\ dv   /   \ dv   /              \ dv   /   /    \ dt   /

/  d   \      /  d   \         //  d   \   /  d   \      /  d   \   \
- |---- X| ux - |---- T| ut + ut ||---- Y| - |---- X| ux - |---- T| ut|
\ dt   /      \ dt   /         \\ du   /   \ du   /      \ du   /   /

//  d   \   /  d   \      /  d   \   \
+ ux ||---- Y| - |---- X| ux - |---- T| ut|) C
\\ dv   /   \ dv   /      \ dv   /   /

In Hickman's package, the names for the infinitesimals are indexed
quantities in the maple environment.  Such a quantity (M say)
can be converted into Maple's function notation by making the substitution
M=op(0,M)(op(M)).  Thus to carry out this conversion on the system above
we make the substitutions:

> neweq:=subs(X=op(0,X)(op(X)),T=op(0,T)(op(T)),Y=op(0,Y)(op(Y)),Z=op(0,Z)(op(
Z)),B=op(0,B)(op(B)),C=op(0,C)(op(C)),[eq]);

neweq := [

/  d               \   /  d               \      /  d               \
- |---- Y(x, t, u, v)| + |---- X(x, t, u, v)| ux + |---- T(x, t, u, v)| ut -
\ dx               /   \ dx               /      \ dx               /

//  d               \   /  d               \      /  d               \   \
ux ||---- Y(x, t, u, v)| - |---- X(x, t, u, v)| ux - |---- T(x, t, u, v)| ut|
\\ du               /   \ du               /      \ du               /   /

- (C(u) ut + B(u))

//  d               \   /  d               \      /  d               \   \
||---- Y(x, t, u, v)| - |---- X(x, t, u, v)| ux - |---- T(x, t, u, v)| ut|
\\ dv               /   \ dv               /      \ dv               /   /

/  d               \   /  d               \
+ |---- Z(x, t, u, v)| - |---- X(x, t, u, v)| (C(u) ut + B(u))
\ dt               /   \ dt               /

/  d               \          /  d               \
- |---- T(x, t, u, v)| ux + ut (|---- Z(x, t, u, v)|
\ dt               /          \ du               /

/  d               \                    /  d               \
- |---- X(x, t, u, v)| (C(u) ut + B(u)) - |---- T(x, t, u, v)| ux) + ux (
\ du               /                    \ du               /

/  d               \   /  d               \
|---- Z(x, t, u, v)| - |---- X(x, t, u, v)| (C(u) ut + B(u))
\ dv               /   \ dv               /

/  d               \
- |---- T(x, t, u, v)| ux),
\ dv               /

/  /  d      \      /  d      \\   /  d               \
Y(x, t, u, v) |- |---- C(u)| ut - |---- B(u)|| + |---- Z(x, t, u, v)|
\  \ du      /      \ du      //   \ dx               /

/  d               \                    /  d               \
- |---- X(x, t, u, v)| (C(u) ut + B(u)) - |---- T(x, t, u, v)| ux + ux (
\ dx               /                    \ dx               /

/  d               \   /  d               \
|---- Z(x, t, u, v)| - |---- X(x, t, u, v)| (C(u) ut + B(u))
\ du               /   \ du               /

/  d               \                         /  d               \
- |---- T(x, t, u, v)| ux) + (C(u) ut + B(u)) (|---- Z(x, t, u, v)|
\ du               /                         \ dv               /

/  d               \                    /  d               \
- |---- X(x, t, u, v)| (C(u) ut + B(u)) - |---- T(x, t, u, v)| ux) - (
\ dv               /                    \ dv               /

/  d               \   /  d               \      /  d               \
|---- Y(x, t, u, v)| - |---- X(x, t, u, v)| ux - |---- T(x, t, u, v)| ut +
\ dt               /   \ dt               /      \ dt               /

//  d               \   /  d               \      /  d               \   \
ut ||---- Y(x, t, u, v)| - |---- X(x, t, u, v)| ux - |---- T(x, t, u, v)| ut|
\\ du               /   \ du               /      \ du               /   /

+

//  d               \   /  d               \      /  d               \   \
ux ||---- Y(x, t, u, v)| - |---- X(x, t, u, v)| ux - |---- T(x, t, u, v)| ut|
\\ dv               /   \ dv               /      \ dv               /   /

) C(u)

]

The infinitesimals are now expressed in maple function notation.
The above command could have also been carried out
using > neweq:= subs(map(proc(z) z=op(0,z)(op(z)) end, [X,T,Y,Z,B,C]), [eq]);

neweq should be saved to a file now, and a new Maple
session started since Hickman's array enviroment otherwise causes problems
for later calculations.
Once in the new session, we read in the conversion
programs:

to apply the map2rep procedure (see 3.6) it is necessary to have B and
C as functions of all the variables before using the procedure, so

> sys1:= subs(C(u)=C(x,t,u,v), B(u)=B(x,t,u,v), neweq):

> sys1:=map2rep(sys1,[X,T,Y,Z,B,C],[x,t,u,v,ux,ut]);

sys1:=

[- V(3, 1) + V(1, 1) x5 + V(2, 1) x6 - x5 (V(3, 3) - V(1, 3) x5 - V(2, 3) x6)

- %1 (V(3, 4) - V(1, 4) x5 - V(2, 4) x6) + V(4, 2) - V(1, 2) %1

- V(2, 2) x5 + x6 (V(4, 3) - V(1, 3) %1 - V(2, 3) x5)

+ x5 (V(4, 4) - V(1, 4) %1 - V(2, 4) x5),

V(3) (- V(6, 3) x6 - V(5, 3)) + V(4, 1) - V(1, 1) %1 - V(2, 1) x5

+ x5 (V(4, 3) - V(1, 3) %1 - V(2, 3) x5)

+ %1 (V(4, 4) - V(1, 4) %1 - V(2, 4) x5) - (V(3, 2) - V(1, 2) x5

- V(2, 2) x6 + x6 (V(3, 3) - V(1, 3) x5 - V(2, 3) x6)

+ x5 (V(3, 4) - V(1, 4) x5 - V(2, 4) x6)) V(6)              ]

%1 :=                          V(6) x6 + V(5)

To apply the standard form package these need to be equations:

> sys1:= map(proc(z) z=0 end, sys1):

Finally, to obtain the determining equations, we polynomially
decompose the above system on x5 and x6 (i.e, ux and ut) using
the sfprogs procedure polyd (see section 7.2):

> polyd(sys1,[x5,x6]);

[- V(3, 1) - V(5) V(3, 4) - V(1, 2) V(5) + V(4, 2) = 0,

V(1, 1) - V(3, 3) + V(4, 4) - V(2, 2) = 0,

- V(1, 2) V(6) + V(2, 1) + V(5) V(2, 4) - V(6) V(3, 4) + V(4, 3) - V(1, 3) V(5)

= 0,

V(1, 3) - V(2, 4) = 0, V(6) V(2, 4) - V(1, 3) V(6) = 0,

- V(3) V(5, 3) - V(3, 2) V(6) + V(4, 1) - V(1, 1) V(5)

+ V(5) (V(4, 4) - V(1, 4) V(5)) = 0,

- V(2, 1) - (- V(1, 2) + V(3, 4)) V(6) - V(1, 3) V(5) - V(5) V(2, 4) + V(4, 3)

= 0,

- V(3) V(6, 3) - (- V(2, 2) + V(3, 3)) V(6) - V(1, 1) V(6) - V(5) V(1, 4) V(6)

+ V(6) (V(4, 4) - V(1, 4) V(5)) = 0,

2
- V(2, 3) + V(1, 4) V(6) = 0, V(2, 3) V(6) - V(6)  V(1, 4) = 0,

- V(6) V(2, 4) - (- V(1, 3) - V(2, 4)) V(6) - V(1, 3) V(6) = 0]

These then are the determining equations for the symmetry generators
of the nonlinear telegraph system.  They are, in fact, quite difficult
to reduce to standard form, and we go through this computation in
section 9.2.  We also mention that the 2 original determining equations
(before polynomial decompostion) could have been entered directly into the
standard_form provided that the expressions:
[V(1),1,2,3], [V(2),1,2,3], [V(3),1,2,3], [V(4),1,2,3],
expressing the dependency info was appended to the above equations
(i.e. that each V(i) is independent of x4 and x5) are added to the system list.
In addition the dependency declarations [V(5),3], [V(6),3]
and are inserted in the class_eqns list since this is a classification problem
(see section 9.2).

Finally the Hickman package can also be used in batch mode:

Batch template for the Hickman Program Symmetry:

domain([u],B);          # dependency of classification variable
domain([u],C);          # dependency of classification variable
varlist   := [x,t];     # independent variables
INDEPCOEFF:= [X,T];     # their corresponding infinitesimals
dependvar := [u,v];     # dependent variables
DEPCOEFF  := [Y,Z];     # their corresponding infinitesimals
eqn:= [Diff(v,t)=Diff(u,x),Diff(v,x)=C*Diff(u,t)+B];
# differential equations (solved-form)

vars:=op(varlist),op(dependvar);
domain([vars],op(INDEPCOEFF),op(DEPCOEFF));
INDEPCOEFF:=[seq(op(i,INDEPCOEFF)[vars],
i=1..nops(INDEPCOEFF))];
DEPCOEFF:=[seq(op(i,DEPCOEFF)[vars],
i=1..nops(DEPCOEFF))];
eqn:=jetform(eqn);
eqn:= geneqn(eqn);     # unpolynomially decomposed determining system

end of template

4.4 INTERFACE TO THE MAPLE PACKAGE LIESYMM

The map2rep procedure explained in section 3.6 allows for an easy
interface between standard form and the maple symmetry package
liesymm.  It can be used to convert the determining equations
produced by the liesymm package into repeated index notation
which can then be analyzed by the standard form programs.
Documentation for the liesymm package can be found in [HaE71,
CamDevFee92],  as well as online documentation.  We illustrate the
procedure with the example below.  Note that this is in release 2 of
Maple V.

As an example, we find the determining equations for the
linear heat equation, using the liesymm package, and then
use the map2rep procedure to put them in repeated index notation.
To begin with, we start a maple session, and read in the
liesymm package using with:

> with(liesymm);

Our equation is the linear heat equation:

> sys:=[Diff(u(x,t),t)-Diff(u(x,t),x,x)=0];

/  d         \
sys := [|---- u(x, t)| - Diff(u(x, t), x, x) = 0]
\ dt         /

Note that when using liesymm, it is necessary to write the DEs
using maple's inert differential operator Diff.  To find the
determining system we use the liesymm procedure determine.

> dets:=determine(sys,V,[u(x,t)],k);
2                      2
d                      d
dets := {----- V2(x, t, u) = 0, ----- V1(x, t, u) = 0,
2                      2
du                     du

/    2              \
d                  |   d               |
---- V1(x, t, u) = - |------- V2(x, t, u)|,
du                  \ du dx             /

2
d                   d
----- V3(x, t, u) = ---- V3(x, t, u),
2                 dt
dx

2                  /    2              \
d                   |   d               |   /  d             \
----- V1(x, t, u) = 2 |------- V3(x, t, u)| + |---- V1(x, t, u)|,
2                  \ du dx             /   \ dt             /
dx

2
d                 /  d             \     /  d             \
----- V2(x, t, u) = |---- V2(x, t, u)| - 2 |---- V1(x, t, u)|,
2                \ dt             /     \ dx             /
dx

d                     d
---- V2(x, t, u) = 0, ---- V2(x, t, u) = 0,
du                    dx

2                  /    2              \
d                   |   d               |
----- V3(x, t, u) = 2 |------- V1(x, t, u)|}
2                  \ du dx             /
du

The determine procedure takes 4 arguments.  The first is the set or
list of differential equations.  The second is a name to give to
the infinitesimal generators (i.e. the dependent variables in the
determining equations).  We have used V.  The third is a list of
the dependent and independent variables entered in the form of
dependent variables as functions of the independent variables
(as above).  The fourth is a name for the extended variable names.
The infinitesimals are numbered such that,
if there are n independents, and m dependents, then the first n
infinitesimals are for the independents, numbered in the order
they appear in the argument lists of the dependents, and the
remaining m are for the dependents in the order they were enetered in
the list.

Now that we have the determining equations, we can convert them
to repeated index notation using the map2rep procedure.  However,
we cannot enter dets directly, because the liesymm procedure
uses an inert differentiation operator, that appears in the
determining equations dets.  To force these operators to become
active, we must use the liesymm procedure value().  So, we enter

> newdets:=map2rep(value(dets),[V1,V2,V3],[x,t,u]);

newdets := {V(1, 3, 3) = 0, V(2, 3, 3) = 0, V(2, 1) = 0, V(2, 3) = 0,

V(2, 1, 1) = - 2 V(1, 1) + V(2, 2), V(1, 3) = - V(2, 1, 3),

V(3, 2) = V(3, 1, 1), V(3, 3, 3) = 2 V(1, 1, 3),

V(1, 1, 1) = V(1, 2) + 2 V(3, 1, 3)}

and these are the determining equations in repeated index notation.

5. INTEGRATION OF THE EQUATIONS

The focus of this package has been on canonically simplifying
systems of PDEs.  Other packages have gone a long way in
implementing methods for integrating PDEs.  It is hoped to
develop interfaces to such packages in the future to take advantage
of their integrators (e.g. the powerful package of Wolf [Wolf89, Wolf91, Wolf92]).
In [Reid91c] we present an algorithm for integrating systems on the
basis of directed graphs associated with their standard forms,
and in [Reid92b] we give an algorithm for solving systems of PDEs in
their coefficient fields by using standard_form to recursively
decouple and solve ODEs.

5.1 USE OF taylor_exp

Occasionally systems of PDEs possess polynomial solutions.
In particular often the determining systems for infinitesimals
possess solutions which are second order (or sometimes third
order) polynomial functions.  By using taylor_exp to third order
(or sometimes 4-th order) one can identify whether the series
has terminated at order 2 (or order 3) in view of the observation that:

If all the order K terms vanish in the taylor series representation
to order K of a function f(x) about an arbitrary point x=x0,
(i.e. they vanish for all values of x0)
then the series terminates at order K-1.

e.g. in the one variable case:  no K-th order terms for any x0 in
(K)
f(x) = f(x0) + f'(x0) (x-x0) + ... + f  (x0) (x-x0)^K / K! + O(K+1)

(K)                          (J)
implies f   (x0) = 0 for all x0.  So f  (x0) = 0 for J > K
and f(x) is a polynomial of degree at most K.

e.g. y''(x) = 0 has initial data y(x0) = a1, y'(x0) = a2
and taylor series to order 2: y(x) = a1 + a2 (x-x0)  + O(3)

This obviously isn't true if the expansion is about
particular point e.g. the series for sin(x) = x + x^3/6 + O(5)
doesn't terminate even though it has no 4-th order terms.

A completely algorithmic method for finding all polynomial/rational
solutions of linear PDE systems with polynomial/rational coefficients
or more generally all solutions belonging to some extension of the
rational function field is given in [Reid92b].

5.2 INTERACTIVE INTEGRATION
At any stage you can integrate simple PDEs by hand or by
a package (e.g. maple's dsolve).

e.g. V(3,1,1) = 0 integrates to
V(3) = V(5) + V(6)*x1,  V(5,1) = 0, V(6,1) = 0
Appending the above 3 equations to the system you wish to
standard_form or solve may help in the process.

6. CHANGING VARIABLES IN SYSTEMS OF PDEs

The procedure transform can be used to perform a general (invertible)
transformation of dependent and independent coordinates on a system
of PDEs.  The basic formula for calculating the prolonged transformations
on partial derivatives of the dependent coordinates with respect to the
independent coordinates is given in [BlK89].  This
highly recursive formula is the core of the transform procedure, which
can express a given PDE in any coordinates given an invertible
transformation from the original coordinates.  It is also possible
to display the little and big jacobians of the transformation.  The
little jacobian is the total differential (i.e. regarding the dependent
variables as functions of the independents) of the transformation for
the independent variables with respect to the new independent variables.
The big jacobian is the jacobian of the entire transformation
with respect to the independent and dependent coordinates.

The transform procedure works in regular maple notation, not
repeated index notation.  Hence coordinate transformations
can be done in a more natural way.  The transform procedure can
be used together with the conversion procedures described in
section 5.6 to use the standard form programs on the new
equations, or to change variables on the standard form of a
system.

The basic syntax of the transform procedure is

> transform(eqns,newindeps,indeptrans,deptrans);

where:

eqns is an equation or system of equations to be transformed;
newindeps is a list of the new independent variables;
indeptrans is a list of the transformations on the independent variables;
deptrans is a list of the transformations for the dependent variables

The transformations in indeptrans
and deptrans must be expressed as equalities, with the transformation,
in terms of the new variables, on the right hand side, and the old
variable to be replaced, on the left hand side.  For instance,
if the existing independent variable and dependent variable
was x and u(x) respectively,
and we wished to make a general linear change to new independent and
dependent variables y and v(y), then indeptrans:=[x=a11*y+a12*v(y)],
and deptrans:=[u(x)=a21*y+a22*v(y)].  The dependencies must
be included for both the new and old dependent variables.
Procedure transform will return the system of equations eqns, expressed
in the new coordinates y, v(y).

To have transform return either the little or big jacobian,
a fifth argument is included whose value is either little, big,
or biglittle to have either the little jacobian, big jacobian, or
both, respectively, returned.  If either big or biglittle is the
fifth argument, then a sixth argument is included, which must
be a list of the new dependent variables, written without dependencies
in this case.  If any of these options are invoked, transform will return
a list of 2 or 3 elements whose first element is the transformed equations
and whose remaining elements are the requested jacobians.  These options
and the general use of transform are demonstrated in the examples below.

Example A:

As a first example, we wish to transform the equation ut - 2 t uxx=0,
under the coordinate transformation x=ay, t=K(s), and u(x,t)=v(y,s),
where K(s) is an arbitrary function and a is a parameter.

> sys:=[diff(u(x,t),t)-2*t*diff(u(x,t),x,x)=0]; # define the system

sys := [D[2](u)(x, t) - 2 t D[1, 1](u)(x, t) = 0].

Then, we invoke transform, to create the transformed system in
the new coordinates

>  tsys:=transform(sys,[y,s],[x=a*y,t=K(s)],[u(x,t)=w(y,s)]);

D[2](w)(y, s)     K(s) D[1, 1](w)(y, s)
tsys := [------------- - 2 --------------------- = 0]
d                        2
---- K(s)                 a
ds

Example B:

Here, we will transform the equation H(Yx) Yxx = Yt where H is
an arbitrary function, by the hodograph x=Y', t=t', Y=x'.
Also, we will display the little and big jacobians.

We enter the system

> sys:=[H(diff(Y(x,t),x))*diff(Y(x,t),x,x)=diff(Y(x,t),t)];

sys := [H(D[1](Y)(x, t)) D[1, 1](Y)(x, t) = D[2](Y)(x, t)]

and then invoke transform

> transform(sys,[y,s],[x=Z(y,s),t=s],[Y(x,t)=y],biglittle,[Z]);

table([
1
H(-------------) D[1, 1](Z)(y, s)
D[1](Z)(y, s)                       D[2](Z)(y, s)
1 = [- --------------------------------- = - -------------]
3              D[1](Z)(y, s)
D[1](Z)(y, s)
[        1           ]
[  -------------   0 ]
[  D[1](Z)(y, s)     ]
2 = [                    ]
[   D[2](Z)(y, s)    ]
[ - -------------  1 ]
[   D[1](Z)(y, s)    ]
[ 0  0  1 ]
[         ]
3 = [ 0  1  0 ]
[         ]
[ 1  0  0 ]
])

The first element is the transformed equation, and the second and
third elements are the little and big jacobians respectively. Note
that the system is linearized when H(Zy)=1/(Zy)^2.

7. TOOLBOX OF INTERACTIVE FUNCTIONS
In future we'll be making some of the internal routines of
the sfprogs package available for interactive work.
Clearderivdep, below, is a useful example of such a routine.

7.1 Clearderivdep

Given a system in standard form sf, then the command

> Clearderivdep(sf, sys);

clears the list of equations sys of all dependencies on sf.
(sf and sys must be entered as lists of equations in repeated
index notation).

Example:

> Clearderivdep([V(3)=V(5)+V(6)*x1,V(5,1)=0,V(6,1)=0],
[V(3,1,1)=0, V(3,1) = V(3,1)^2 + V(4)]);

yields:

2
[0 = V(6)  + V(4) - V(6)]

(note that the first equation simplified to 0=0 and was omitted
from the output).  Also see Section 5.2.

7.2 polyd

Procedure polyd

Input: list of equations which is polynomial in pvars, and
the list of variables pvars.

Output: list of polynomially decomposed eqns with respect to pvars.

Example:

> polyd([t1^2*y+t1*t2*k=h ,(t1+t2)*u=t1*y],[t1,t2]);

yields

[- h = 0, y = 0, k = 0, u = 0, u - y = 0]

7.3 sfdiff

sfdiff(exp, string)
where exp is an expression or equation in repeated index notation
and string is a string of independent variables from the set of
independent variables {x1,x2, ... ,xm}.

Example:

> sfdiff(V(1,1,1)^2, x1);

2 V(1, 1, 1) V(1, 1, 1, 1)

> sfdiff(V(2)^2 + V(1,3)=0, x1,x3);

2 V(2, 3) V(2, 1) + 2 V(2) V(2, 1, 3) + V(1, 1, 3, 3) = 0

8. STRATEGIES TO DEAL WITH MEMORY EXPLOSION

In general, memory explosion seems to occur when the system to which standard_form is applied is
very complex (algebraically), or has a large number of equations.
We discuss strategies for dealing with these problems
(also see the the interactive programs of Kersten [Ker87] and the
automatic programs of Wolf [Wolf89, Wolf91, Wolf92] for interesting
approaches to the problem of controlling memory swell).

Changing the equation-flow settings to their ideal can be a
system dependent process. Often adjusting the knobs to
their Ultra-Conservative setting (i.e. one equation at a time)
will lead to termination using a reasonable quantity of memory,
However, some which cases blow up in the conservative setting
terminate in a reasonable time using little memory in the default
setting. This occurs in systems having one or more large equations which
dramatically reduce the complexity of the system once solved and
substituted for. These situations are not easy to detect, and
currently the best choice is an experimental approach.
Another (desparate strategy!) would be to change that part of the
code which controls the way in which easy equations are chosen from
unclassified equations (see 3.3.2).

8.2 VARYING THE ORDERING

Changing the ordering using the janet_order option can lead
to drastic reductions in memory and run time.  We are
not aware of any general strategy for choosing the ordering to obtain such
reductions. The main reason for this is the difficulty in predicting
the path on which the solution will be arrived at. Again the best
method currently seems to be experimental.

Sometimes ordering with respect to the dependent variables
last, and differentiation variable first, minimizes memory expansion
and produces the best run times.  Sometimes total ordering after the
differentiation variable ordering can be effective.

e.g. 3 independent variables , 3 dependent variables
set
> ordertab1:= [[1,0,0, 0,0,0],
[0,1,0, 0,0,0],
[0,0,1, 0,0,0]]:
> standard_form(  , janet_order = ordertab1 )

and if this does not work, one could permute order of the three
lists inside ordertab1.

8.3 CLASSIFICATION VARIABLES AND DIFFERENTIAL EXTENSIONS

(or polynomials are nicer)

Maple and other computer algebra systems can always reduce
polynomial or rational functions to normal form. However, when
a system possesses a number of non-polynomial functions such as
sin(x1), exp(x3),and x4^(1/2), there is a great deal of memory
consumption involved in simplification and manipulation of such
expressions (in general Maple has no canonical
way of simplifying expressions involving such functions).

One way to get eliminate a nonpolynomial function is to introduce a new
classification variable which represents the function, and
specify its behavior with a differential constraint ( i.e.
define the function by a differential extension).
For example consider a system of PDEs with 5 dependent variables, 5
independent variables, and a painful recurrence of the factor
x4^(1/2). Substituting x4=V(6)^2 into the system and introducing
the class equation

class_eqns=[[V(6),4],V(6,4)=1/(2*V(6))]

which defines the derivative properties of V(6) eliminates the
square root and replaces it with pleasant rational occurences of
V(6).

Alternatively, a non-standard method (courtesy of ADW),
giving a linear form for
the differential constraint could be used. For example:
First substitute x4^(1/2)=V(6) then set:

class_eqns=[[V(6),4],V(6,4)=V(6)/(2*x4)]

(i.e. V(6)=c2*x^(1/2) where c2 is an arbitrary constant).
However, one must ensure that none of the pivots
vanish.  For example, if one of the pivots was: V(6)^2-x4=0 then an invalid
operation has occurred and the calculation must be done again with a
suitably different class equation.

8.4 THE BIG GUNS - STORAGE AND USEREASY

Okay... So we're desperate now. Don't sweat it yet, there are still
a couple of ways to obtain the answer. However, do not read further if
you don't expect to have to work for it.

Still Here? Okay, here it is you glutton for punishment:
One can use the storage options (see 3.4) to store the system at the
most recent iteration (store) or if you don't mind clogging up
your filespace at every iteration (storeall).  You can also set
the equation flow parameters at ultra-consevative. And hack
at the system a piece at a time. In order to do this, one requires
a working knowledge of Maple's list handling features such as
op(), susbop(), map() to manipulate the lists (also see section 7).
Oh, and lest we forget that many of the options mentioned previously
in this manual will come to bear on this strategy, so go read the first
7 chapters again.

One strategy which GR has used with some success is to use the store option
to recover the system in at some chosen stage or at the last iteration before
the WOD (our favorite phrase for memory
explosion).  Then choose the simplest equations from the system, throw the
ugly ones out and start the original system together with these simplest
equations and perhaps use the usereasy option (see 3.4.2).

You're Back Already?!? OK I trust you.
The storeall option (3.4) allows one to retrieve the stored system
before the WOD was reached. So you've now run the code in ultra-conservative setting (see (3.3.2)) with the infolevel at its maximum
(i.e. strategy = [1.75, 0.9, 1.0]).
You should have as many StandardStorage_i.m files in your directory as iterations the program completed (Usually 50~80 of them for large systems) DON'T PANIC! This is normal. So now you begin another Maple session (Note: it is important to start with a fresh slate if the program went down in flames) load the code, and use readstored to load the most recent version (highest numbered store-file) with:

> L:=readstored(StandardStorage_'i'):

where 'i' is the highest number available.
So now you have the system, i.e. OT,Ez,UnC,_pivs and possibly NL all stored
in a list with name L (you can refer to 4.3 to see which is which).
A beneficial
move at this point would be to issue the command:

> assign(L);

which will now assign the list of one term equations to the variable OT,
the list of Easy equations to the variable Ez, etc.
Now comes the fun part! Try:

> map(length,Ez);
> map(length,OT);

You should now have a couple of lists displayed, and they possess the
length (a crude measure of complexity) of the corresponding elements
in Ez and OT respectively. You can look at each individual element
(say the i'th element of Ez) with:

> Ez[i];
or
> op(i,Ez);

Good now we're ready for the real work. Since there was an explosion,
there is likely a number of elements of Ez with relatively large length.
Start a new variable:

> tough:=NULL;

and pull undesirable or long elements out of Ez from the last to first
using:

> tough:=tough,Ez[i]:
> Ez:=subsop(i=NULL,Ez):

In this way you are moving long or complicated equations out of Ez
and appending them to tough ( it is imperative not to lose equations on
the way ) Okay, all the real nasties out? You can check this with
another map command:

> map(length,Ez)

Good! Now turn Ez into an appendable variable by typing:

> Ez:=op(Ez):
now
> map(length,UnC);

Look for the simple equations in the UnClassified (it's good to
manually check these also; sometimes even long equations can be simple)
and a good rule is the fewer the classification variables, and the smaller
the coefficients (as far as the xi polynomials go) the better.
To transfer them to the Ez list do as follows:

> Ez:=Ez,UnC[i]:
> UnC:=subsop(i=NULL,UnC):

Put it all together and give it a test drive! Okay,
first step : We know that we like the Ez equations and the OT equations
so let's put them together:

> Ez:=op(OT),Ez:

Now we want to know how many good equations there are, we use:

> nops([Ez]);

Remember the number you now see, we call it numeasy for reference.
Now we put everything together into a variable called sys:

> sys:=[Ez,op(UnC),tough,op(NL)]:
> firstpivs:=_pivs:
and finally:
> save sys,firstpivs,temporaryfile;

We got it right where we want it now! Remember that number? numeasy?
You need it in the following statement:

> sf:=standard_form(sys,usereasy=numeasy,storeall=Attempt2):

This will now fire off the system and store the system at the i'th
iteration in Attempt2_i.m, assuming that the first numeasy number of
equations are fairly easy, and will solve those first. Hopefully
the nasty equations which caused the blowup will now simplify to a
reasonable, workable level.

SPECIAL NOTES:
If it is important to know the pivots of the system, some of those are
now stored in temporaryfile in the variable firstpivs, and the others
since the re-run are being built up in _pivs. Both lists are important
and result from the current problem being solved.

In the new code we have not yet found a linear problem of practical
origin for which we've had to resort to this drastic measure, so be sure
to exhaust the other possibilities first. Also if you find a problem of
this magnitude, we would be interested in knowing of it.

9. NONTRIVIAL EXAMPLES

9.1 A BIG EXAMPLE: THE MHD EQUATIONS

A description on how to do this type of problem
was already introduced in a simple example in sections 4.2.1 and 4.1.1
of the manual.

As an example of our programs working on a nontrivial (large)
example we generate the determining system for symmetries of the
the magneto-hydro-dynamic (mhd) equations (below).  The computation
of the standard form, although a large one, terminates in about
50 minutes, on a Silicon Graphics machine, and remarkably never
uses more than about 2.2 Megabytes of RAM (Sparc Stations give
comparable performance).

MHD EQUATIONS:

f   +  (v.Grad) f  + f Div v  =  0
t

f ( v   +  (v.Grad) v )  +  Grad ( p  +  (H.H)/2 ) - (H.Grad ) H = 0
t

H    +  (v.Grad) H  +  H (Div v)   -  (H.Grad) v  =  0
t

Div H  =  0

K                        K
( p / f  )    +  (v.Grad) (p / f  )  =  0
t

Above pressure=p=p(x,y,z,t), mass density = f = f(x,y,z,t),
coefficient of viscosity = K, vector fluid velocity v=v(x,y,z,t)=(Vx,Vy,Vz)
and magnetic field H=H(x,y,z,t)=(Hx,Hy,Hz).
Also Grad = (d/dx,d/dy,d/dz), and . denotes dot product.

For the macsyma program we use the correspondence:
x=x[1], y=x[2], z=x[3], t=x[4], f = u[1], p = u[2],
Vx=u[3], Vy=u[4], Vz=u[5],
Hx=u[6], Hy=u[7], Hz=u[8].
with corresponding symmetry generator:
L =   ETA[1] d/dx  + ETA[2] d/dy + ETA[3] d/dz + ETA[4] d/dt
+ PHI[1] d/df  + PHI[2] d/dp
+ PHI[3] d/dVx + PHI[4] d/dVy + PHI[5] d/dVz
+ PHI[6] d/dHx + PHI[7] d/dHy + PHI[8] d/dHz

The corresponding data file for the macsyma program is stored under
mhd.dat in the same directory as sfprogs.  Batching a file containing
the (macsyma) commands:

(c1) batch("mhd.dat");
(c2) batch("maxconv");
(c3) quit();

produces a file lode containing the mhd determining equations in
macsyma syntax.
Following the same path as in section 4.2.1 and 4.1.1 of the manual:
Once the mhd determining equations have been generated the file lode is
stripped of ' using an editor.
We start a maple session by loading the conversion programs and the
stripped lode:

Then give the determining system a name:

> msys:= ":

The conversion of msys to repeated index notation is
accomplished by applying the function max2rep:

> sys:= max2rep(msys);

Notice the correspondence:
ETA(i)=V(i), i=1,2,3,4
PHI(1)=V(5),PHI(2)=V(6),...,PHI(8)=V(12),
and
x1=x[1], x2=x[2], x3=x[3], x4=x[4],
x5=u[1], x6=u[2], x7=u[3],   ...  , x12=u[8]

and save the converted system to a file (called mhd.sys here),
which you will find in the same directory as sfprogs.

> save sys, mhd.sys;

We start a new maple session, read the standard form package

and the converted system

> read mhd.sys;

Next the determining system is reduced to standard form:

> sf:= standard_form(sys);

The resulting standard form obtained in about 1 hour of cpu
on a SGI machine, which we have not displayed
has 133 equations - 102 of which are one-term equations
The original system had 42 one-term equations.

The shortened pivots are

> spivs:= shortenpiv(_pivs);
2      2      2
[K, x11, x10, x12, x5, x8, x7, x6, x10  + x12  + x11 ,

- x8 x12 + x9 x11,- 2 + K, - x7 x12 + x9 x10,

4        2    2      2    2        2    2        4
x11  + 4 x11  x12  - x10  x11  - 2 x12  x10  - 2 x10 ,

- 1 + K, - 3/2 + K,

2      2          2     2          2
x11  + x12  + 1/4 x10 , x11  + 1/2 x12 ,

4          2    2      2    2          4      2    2
x11  - 1/2 x11  x12  - x10  x11  + 5/2 x10  + x12  x10 ,

2        2
x11  - 2 x10 ]

So this standard form is valid if these shortened pivots
do not vanish identically that is the standard form is
valid for K # 0,1,3/2,2.

We calculate its initial data:

> id:= initial_data(sf);

yielding

[[V(1) = a1, V(2) = a2, V(3) = a3, V(4) = a4, V(4, 4) = a5, V(5) = a6,
V(7) = a7, V(8) = a8, V(9) = a9, V(10) = a10, V(11) = a11,
V(11, 12) = a12, V(12) = a13],                                []]

So if K # 0,1,3/2,2 the Lie group is 13-dimensional (i.e. the Lie symmetry
group is finite-dimensional).

The commutator table for the symmetries L(1),L(2),...,L(13) corresponding
to this initial data a1,a2,...,a13 respectively is given by:

> ct:= fin_com_table(sf,id);

The table is a 13 X 13 array, which we have not displayed because of
its complexity.  The quantities X1,X2,X3,...,X11,X12 appearing in the
table are the coordinates
of the arbitrary point through at which the initial data is posed.
This table was calculated before the exact form of the generators was
known and we will give it in simplified form at the end of the session.

To calculate the generators we use the command finite_generators,
which uses the taylor_exp routine (see 4.1.1).

In particular

> fg1:= finite_generators(sf,id,1);

gives the taylor expansion of the generators to order 1:

fg1 := [L(1) = Dx(1), L(2) = Dx(2), L(3) = Dx(3), L(4) = Dx(4),

L(5) = h1 Dx(1) + h2 Dx(2) + h3 Dx(3) + h4 Dx(4),

/       h1        X7 h4\         /    X8 h4        h2 \
L(6) = |- 1/2 ---- + 1/2 -----| Dx(1) + |1/2 ----- - 1/2 ----| Dx(2)
\       X5          X5 /         \      X5         X5 /

/       h3        X9 h4\         /     h5 \             h7 Dx(7)
+ |- 1/2 ---- + 1/2 -----| Dx(3) + |1 + ----| Dx(5) - 1/2 --------
\       X5          X5 /         \     X5 /                X5

h8 Dx(8)       h9 Dx(9)
- 1/2 -------- - 1/2 --------,
X5             X5

L(7) = h4 Dx(1) + Dx(7), L(8) = h4 Dx(2) + Dx(8),

L(9) = h4 Dx(3) + Dx(9), L(10)= ..., L(11)=..., L(12)=..., L(13)=...]

where h1=(x1-X1), h2=(x2-X2), etc and L(10), L(11), ... , L(13) are
long expressions.

These generators contain terms of order 1 in the hi's so by section 5.1
of the manual we can't be sure that they are exact (i.e. that there are
no higher order terms in the expansion).

Expanding the generators to order 2

> fg2:= finite_generators(sf,id,2);

we find that there are no order 2 terms in the generators
(by inspection or by checking that fg1=fg2 using the Maple
boolean command evalb(fg1=fg2)).
So the generators fg1 are exact (see section 5.1 of the manual).

In terms of the independent variables x1,...,x12 and the
point (X1,X2,...,X12) at which the initial data is posed
the generators are given by

> subslist1:= [h1=x1-X1,h2=x2-X2,h3=x3-X3,h4=x4-X4,
h5=x5-X5,h6=x6-X6,h7=x7-X7,h8=x8-X8,
h9=x9-X9,h10=x10-X10,h11=x11-X11,h12=x12-X12]:
> fg1:= subs(subslist1, fg1):

and look rather horrible.

To simplify their appearance we use our freedom to choose the
point (X1,X2,...,X12) at which the initial data is posed.
Usually X1=0,X2=0,...,X12=0 is the best choice but the generators
and standard form are singular when X5=0 and X10=0, so we put
X5=1, X10=1 and the remaining Xi's to zero:

> subslist2:= [X1=0,X2=0,X3=0,X4=0,
X5=1,X6=0,X7=0,X8=0,
X9=0,X10=1,X11=0,X12=0]:

> fg1:= subs(subslist2, fg1);

obtaining the generators

[L(1) = Dx(1), L(2) = Dx(2), L(3) = Dx(3), L(4) = Dx(4),

L(5) = x1 Dx(1) + x2 Dx(2) + x3 Dx(3) + x4 Dx(4),

L(6) = - 1/2 x1 Dx(1) - 1/2 x2 Dx(2) - 1/2 x3 Dx(3)

+ x5 Dx(5) - 1/2 x7 Dx(7)

- 1/2 x8 Dx(8) - 1/2 x9 Dx(9),

L(7) = x4 Dx(1) + Dx(7), L(8) = x4 Dx(2) + Dx(8), L(9) = x4 Dx(3) + Dx(9),

L(10) = x1 Dx(1) + x2 Dx(2) + x3 Dx(3) + 2 x6 Dx(6) + x7 Dx(7) + x8 Dx(8)

+ x9 Dx(9) + x10 Dx(10) + x11 Dx(11) + x12 Dx(12),

L(11) = - x2 Dx(1) + x1 Dx(2) - x8 Dx(7) + x7 Dx(8) - x11 Dx(10) + x10 Dx(11),

L(12) = x3 Dx(2) - x2 Dx(3) + x9 Dx(8) - x8 Dx(9) + x12 Dx(11) - x11 Dx(12),

L(13) = - x3 Dx(1) + x1 Dx(3) - x9 Dx(7) + x7 Dx(9) - x12 Dx(10) + x10 Dx(12)]

These results agree with those given by other investigators (hereman [Here93]).
[in more physical notation x1=x,x2=y,x3=z,x4=t,x5=f,x6=p,x7=Vx,x8=Vy,x9=Vz,
x10=Hx,x11=Hy,x12=Hz so that L(5)=x d/dx + y d/dy + z d/dz + t d/dt etc]

Recall that the commutator table ct was calculated already.  We simplify
it by wisely choosing our initial point using the same substitution as
above:

> L:= proc(k) cat(L,k) end:   # this procedure shortens L(1) to L1, etc.

> ct:= map(simplify, subs(subslist2, op(ct)));

L1 L2 L3 L4 L5      L6       L7      L8      L9   L10  L11    L12    L13

L1 [ 0  0  0  0  L1  - 1/2 L1     0       0       0    L1   L2     0      L3  ]
[                                                                          ]
L2 [ *  0  0  0  L2  - 1/2 L2     0       0       0    L2  - L1   - L3    0   ]
[                                                                          ]
L3 [ *  *  0  0  L3  - 1/2 L3     0       0       0    L3    0     L2    - L1 ]
[                                                                          ]
L4 [ *  *  *  0  L4      0       L1      L2      L3     0    0     0      0   ]
[                                                                          ]
L5 [ *  *  *  *   0      0        0       0       0     0    0     0      0   ]
[                                                                          ]
L6 [ *  *  *  *   *      0     1/2 L7  1/2 L8  1/2 L9   0    0     0      0   ]
[                                                                          ]
L7 [ *  *  *  *   *      *        0       0       0    L7   L8     0      L9  ]
[                                                                          ]
L8 [ *  *  *  *   *      *        *       0       0    L8  - L7   - L9    0   ]
[                                                                          ]
L9 [ *  *  *  *   *      *        *       *       0    L9    0     L8    - L7 ]
[                                                                          ]
L10[ *  *  *  *   *      *        *       *       *     0    0     0      0   ]
[                                                                          ]
L11[ *  *  *  *   *      *        *       *       *     *    0   - L13   L12  ]
[                                                                          ]
L12[ *  *  *  *   *      *        *       *       *     *    *     0    - L11 ]
[                                                                          ]
L13[ *  *  *  *   *      *        *       *       *     *    *     *      0   ]

The i-j th entry in the above table corresponds to [L(i), L(j)].

In summary the exact symmetry generators for the mhd equations
were obtained by applying fully automatic procedures.
The above results are valid for k # 0,1,3/2, 2.  We have separately
run the cases K=1, 3/2, 2 and found that the standard forms
were just the same as substituting the given values of K in the
standard form above - and that we get the same size of group and
so these cases have the same generators as above.

When we ran the program on the case K=0 we found a 14-parameter
invariance group, with the same generators above and an extra generator translation invariance in the pressure:  d/dp.

9.2 A SMALL NASTY EXAMPLE: THE NONLINEAR TELEGRAPH SYSTEMS

The nonlinear telegraph system produces a determining system
which appears quite simple and easy to deal with, but leads to
memory explosion in the default settings for standard_form.
We believe that this system is a tough challenge for any
automatic PDE simplification or integration program (and the
solutions are interesting mathematically and physically! [Reid91b]).

The nonlinear telegraph system is

Vt = Ux ,

Vx = C(U) Ut  +  B(U) .

where C and B are arbitary functions of U.  The determining
system for it is written in repeated index notation, as it
would be entered into standard form, was derived in Section 4.3
and is given below:

sys:= [V(6)*V(3,2)+V(5,3)*V(3)+V(5)^2*V(1,4)+V(5)*(V(1,1)-V(4,4))-V(4,1)=0,
V(6)*(V(3,4)-V(1,2))+V(5)*(V(2,4)+V(1,3))-V(4,3)+V(2,1)=0,
V(6,3)*V(3)+2*V(5)*V(6)*V(1,4)+V(6)*(V(3,3)+V(1,1)-V(2,2)-V(4,4))=0,
V(2,3)-V(6)*V(1,4)=0,
V(5)*(V(3,4)+V(1,2))-V(4,2)+V(3,1)=0,
V(2,2)+V(3,3)-V(1,1)-V(4,4)=0,
V(1,3)-V(2,4)=0,
V(6)*(V(3,4)+V(1,2))+V(5)*(V(1,3)-V(2,4))-V(4,3)-V(2,1)=0];

Note that it is a classification problem.  In addition to the 4
infinitesimals corresponding to the 2 independent and 2 dependent variables
of the telegraph system, there are the additional dependent
variables V(5) and V(6),  corresponding to the arbitrary functions C and B.
They are classification variables, and so are entered in the
list class.

class:= [[V(5),3],[V(6),3]];

Now we will find the standard form.  Note that we are setting the
equation control knobs (see 3.3.2) in the strategy list
at their ultraconservative setting.

> sf:= standard_form(sys, class_eqns=class,strategy=[1.0,1.0,0.9]);

The system parameters are as follows:

There are 4 independent variables:

[x1, x2, x3, x4]

There are 4 dependent variables:

[V(1), V(2), V(3), V(4)]

There are 2 classification variables:

[V(5), V(6)]

bytes used=589668, alloc=393144, time=6.600
Initial State of System:
#OT: 0   #Ez: 0   #NL: 0   #UnC: 8

bytes used=681444, alloc=393144, time=7.666
******End iteration #1  Time:.2.683
#OT: 0   #Ez: 1   #NL: 0   #UnC: 7

bytes used=1012948, alloc=589716, time=10.333
******End iteration #2  Time:.2.050
#OT: 0   #Ez: 2   #NL: 0   #UnC: 6

bytes used=1174184, alloc=589716, time=12.466
******End iteration #3  Time:.3.067
#OT: 0   #Ez: 3   #NL: 0   #UnC: 6

bytes used=1448196, alloc=589716, time=15.633
******End iteration #4  Time:.3.400
#OT: 0   #Ez: 4   #NL: 0   #UnC: 6

bytes used=1737072, alloc=655240, time=19.083
******End iteration #5  Time:.3.300
#OT: 0   #Ez: 5   #NL: 0   #UnC: 5

bytes used=1997928, alloc=655240, time=22.416
******End iteration #6  Time:.5.683
#OT: 0   #Ez: 6   #NL: 0   #UnC: 6

bytes used=2512716, alloc=786288, time=28.266
******End iteration #7  Time:.5.633
#OT: 0   #Ez: 7   #NL: 0   #UnC: 5

bytes used=5970332, alloc=1376004, time=59.416
******End iteration #8  Time:.31.317
#OT: 0   #Ez: 8   #NL: 0   #UnC: 7

bytes used=8773096, alloc=1376004, time=83.216
******End iteration #9  Time:.24.550
#OT: 0   #Ez: 9   #NL: 0   #UnC: 9

bytes used=11655660, alloc=1441528, time=109.000
******End iteration #10  Time:.21.484
#OT: 0   #Ez: 10   #NL: 0   #UnC: 9

bytes used=15004008, alloc=1572576, time=138.116
******End iteration #11  Time:.29.966
#OT: 0   #Ez: 11   #NL: 0   #UnC: 10

bytes used=18418052, alloc=1572576, time=170.883
******End iteration #12  Time:.34.650
#OT: 0   #Ez: 12   #NL: 0   #UnC: 10

bytes used=23125076, alloc=1572576, time=215.566
******End iteration #13  Time:.45.584
#OT: 0   #Ez: 13   #NL: 0   #UnC: 11

bytes used=39105772, alloc=2555436, time=352.183
******End iteration #14  Time:.139.534
#OT: 0   #Ez: 14   #NL: 0   #UnC: 14

bytes used=61606716, alloc=2817532, time=538.583
******End iteration #15  Time:.181.583
#OT: 0   #Ez: 15   #NL: 0   #UnC: 17

bytes used=69920608, alloc=2817532, time=613.500
******End iteration #16  Time:.74.183
#OT: 1   #Ez: 15   #NL: 0   #UnC: 8

bytes used=71817492, alloc=2817532, time=630.733
******End iteration #17  Time:.22.134
#OT: 3   #Ez: 14   #NL: 0   #UnC: 2

bytes used=86684956, alloc=2817532, time=747.966
******End iteration #18  Time:.113.016
#OT: 10   #Ez: 7   #NL: 0   #UnC: 5

bytes used=87367372, alloc=2817532, time=755.000
******End iteration #19  Time:.5.284
#OT: 8   #Ez: 6   #NL: 0   #UnC: 0

bytes used=87688192, alloc=2817532, time=758.750
bytes used=87689740, alloc=2817532, time=759.050
sf := [V(1, 1) = 0, V(2, 1) = 0, V(4, 1) = 0, V(1, 2) = 0, V(2, 2) = 0,

V(4, 2) = 0, V(1, 3) = 0, V(2, 3) = 0, V(4, 3) = 0, V(1, 4) = 0,

V(2, 4) = 0, V(4, 4) = 0, V(3) = 0]

The output illuminates the process standard_form went through in completing
this awesome task.  The first 7 iterations are fairly routine short iterations,
with standard form adding equations to the Ez list 1 equation at a
time (recall that this is because we have set the knobs to their
ultra-conservative setting).  The next 6 iterations up to the 13th
are of moderate length in terms of CPU time, and the UnC list is
growing a bit faster than before, showing that as standard_form adds
the longer equations from the UnC list, more difficult integrability
conditions are arising.  But the 14th iteration takes 139.534
seconds of CPU time, and the UnC list grows by an additional 4, as it
does at the 15 iteration which takes 181.583 seconds of CPU.  The
equations arising at these iterations are causing far more work
form standard_form than in previous iterations.  At the 16th
iteration however, a "miracle" happens and we get a collapse in the
UnC list from 17 to 8 equations, and at the 17th, it collapses
to just 2.  Their is one more difficult iteration at #18,
which takes 113.016 CPU seconds, but the final system has only 13
one term equations.  The final standard form is very simple,
but it was a great deal of work to find it.  Note that while the
output system is quite simple, the pivots, if we were to look at them,
even in simplified form, would take up several pages.  But with Ian Lisle's
ultra sophisticated approach to classification problems given in [Lis92]
these pivots can be simplified to a few lines.  The generic result above
reveals only inspectional symmetries.  However by examining the pivots
and following special cases where some of the pivots vanish truly nontrivial
results are obtained [Reid91b].

For example a nontrivial calculation and result occurs
if B=1/(1+U^2), C=1/(1+U^2)^2

> sf1:= standard_form(sys, class_eqns=[[V(5),3],[V(6),3], V(5)=1/(1+x3^2),
V(6)=1/(1+x3^2)^2],strategy=[1.75,1.0,0.9]);
V(3)
sf1 := [V(1, 1) = ----, V(2, 1) = 0, V(3, 1) = 0, V(4, 1) = 0, V(1, 2) = 0,
x3

V(3)
V(2, 2) = - ----, V(3, 2) = 0, V(4, 2) = 0, V(1, 3) = 0,
x3

V(3)                V(3)
V(2, 3) = - 2 -------------, V(3, 3) = ----, V(4, 3) = 0,
2 2                x3
(1 + x3 )  x3

V(3)                                        V(3)
V(1, 4) = - 2 ----, V(2, 4) = 0, V(3, 4) = 0, V(4, 4) = - ----]
x3                                          x3

> id1:= initial_data(sf1, class_eqns=[[V(5),3],[V(6),3], V(5)=1/(1+x3^2),
V(6)=1/(1+x3^2)^2]);

id1 := [[V(1) = a1, V(2) = a2, V(3) = a3, V(4) = a4], []]

Thus there is a nontrivial 4-dimensional Lie symmetry group
in this case.

The commutator table is

> ct1:= fin_com_table(sf1, id1, class_eqns=[[V(5),3],[V(6),3], V(5)=1/(1+x3^2),
V(6)=1/(1+x3^2)^2]);

[        L(1)                 ]
[ 0  0   ----         0       ]
[         X3                  ]
[                             ]
[         L(2)                ]
[ *  0  - ----        0       ]
ct1 := [          X3                 ]
[                             ]
[                 L(1)   L(4) ]
[ *  *     0    2 ---- + ---- ]
[                  X3     X3  ]
[                             ]
[ *  *     *          0       ]

where X3 can be taken to be any fixed nonzero value (e.g. X3 = 1).

REFERENCES

\documentstyle[12pt]{article}
\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                 %
%                           BIBLIOGRAPHY                          %
%                                                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{thebibliography}{Abc99}

\bibitem{BCGGG89}
Bryant, R. L., Chern, S. S., Gardner R. B., Goldschmidt H. L.
and Griffiths, P. A. (1991), Exterior Differential Systems, Mathematical
Sciences Research Institute Publications {\bf 18} (Springer-Verlag, New York).

\bibitem{BlK89}
Bluman, G.W. and S. Kumei. 1989. {\em Symmetries and Differential
Equations.} Springer, New York.

\bibitem{Boch89}
Bocharov, A. V. and Bronstein, M. L. (1989), Efficiently implementing two
methods of the geometrical theory of differential equations'',
{\it Acta. Appl. Math.} {\bf 16} 143-16.

\bibitem{CamDevFee92}
{\sc J. Carminati, J.S. Devitt} and {\sc G.J. Fee},\,
{\em Isogroups of differential equations using algebraic computing},
J. Sym. Comp. {\bf 14} (1992) 103-120.

\bibitem{Carra87}
Carr\{a}-Ferro, G. 1987, Gr\"{o}bner bases and differential ideals," {\it Proc. AAECC5}, Menorca, Spain, Springer-Verlag, 129-140.

\bibitem{CarraSitt93}
Carr\a Ferro, G. and Sit, W. Y. (1993), On Term-Orderings and Rankings", to appear.

\bibitem{Cart45}
E. Cartan,
{\it
Les syst\{e}mes diff\'{e}rentiels ext\'{e}rieurs et leurs applications
g\'{e}om\'{e}triques}, (Paris, Hermann, 1945).

\bibitem{Champ91}
Champagne, B., Hereman, W. and Winternitz, P. (1991), The computer calculation of Lie point symmetries of large systems of differential equations," {\it Computer Physics Communications} {\bf 66}, 319-340.

\bibitem{Cha90} Char, B.W., K.O. Geddes, G.H. Gonnet, M.B. Monagan, and S.M.
Watt. 1990. (3rd ed.) {\em {\sc Maple} Reference Manual.}
Watcom, Waterloo.

\bibitem{Fed86}
Fedorova, R.N. and V.V. Kornyak. 1986. Determination of Lie B\"acklund
symmetries of differential equations using {\sc FORMAC}.
{\em Comput. Phys. Comm.} {\bf 39}: 93--103.

\bibitem{HaE71} Harrison, B.K. and F.B. Estabrook. 1971. Geometric approach to
invariance groups and solution of partial differential equations.
{\em J.~Math. Phys.} {\bf 12}: 653--666.

\bibitem{HartTker91}
Hartley, D. and Tucker, R. W. 1991. {\it A Constructive Implementation
of the Cartan-Kahler Theory of Exterior Differential Systems}.
J. Symb. Comp. {\bf 12}, 655-657.

\bibitem{Here93}
Hereman, W. 1993.  Review of Symbolic Software for the Computation of
Lie Symmetries of Differential Equations.
Preprint.

\bibitem{Hick93}
Hickman, M. 1993.  The use of MAPLE in the Search for Symmetries. Preprint.

\bibitem{Jan20}
Janet, M. 1920. Sur les syst\emes d'\'equations aux d\'eriv\'ees partielles.
{\em J. de math\/} {\bf 3}: 65--151.

\bibitem{Ker87}
Kersten, P.H.M. 1987. {\em Infinitesimal symmetries: a computational approach.}
Centrum voor Wiskunde en Informatica, Amsterdam.

\bibitem{Kol73}
Kolchin, E. R. (1973), {\it Differential Algebra and Algebraic Groups}, Academic Press, New York.

\bibitem{Lis92}
Lisle, I.G. 1992. {\em Equivalence Transformations for Classes of Differential
Equations\/.} Ph.D Thesis, University of British Columbia.

\bibitem{Mans91}
Mansfield, E. 1991.
{\em Differential Gr\"{o}bner Bases}, Ph.D Thesis, University of
Sydney.

\bibitem{Olliv91}
Ollivier, F. (1991), Standard bases of differential ideals'',
{\it Lecture Notes in Comp. Sci.}, {\bf 508}, 304-321.

\bibitem{Olv86}
Olver, P.J. 1986. {\em Applications of Lie Groups to Differential
Equations.}
Springer, New York.

\bibitem{Ovs82}
Ovsiannikov, L.V. 1982. {\em Group Analysis of Differential

\bibitem{Pank89}
E. V. Pankrat'ev, Acta Applicandae Mathematicae {\bf 16}, 167 (1989).

\bibitem{Pomm78}
Pommaret, J. F. 1978.
{\em Systems of  Partial Differential Equations and Lie Pseudogroups}
(Gordon and Breach science publishers, Inc.)

\bibitem{Pomm90}
Pommaret, J. F. 1990.
{\em Equations aux d\'{e}riv\'{e}es partielles et
th\'{e}orie des groupes} (Institute National de Recherche
en Informatique et en Automatique, Roquencourt).

\bibitem{Reid90}
Reid, G.J. 1990. A triangularization algorithm which determines the Lie
symmetry algebra of any system of PDEs. {\em J. Phys.} {\bf A23}:
L853--859.

\bibitem{Reid91a}
Reid, G.J. 1991. Algorithms for reducing a system of PDEs to standard form,
determining the dimension of its solution space and calculating its
Taylor series solution.
{\em Euro. J. Appl. Maths.\/} {\bf 2}: 293--318.

\bibitem{Reid91b}
Reid, G.J. 1991. Finding abstract Lie symmetry algebras of differential
equations without integrating determining equations.
{\em Euro. J. Appl. Maths.\/} {\bf 2}: 319--340.

\bibitem{Reid91c}
Reid, G. J. and Boulton, A. (1991c), Reduction of systems of differential equations to standard form and their integration using directed graphs," {\it Proceedings of ISSAC '91}, ACM Press, New York, 308-312.

\bibitem{Reid92a}
Reid, G.J., I.G. Lisle, A. Boulton and A. D. Wittkopf. 1992. Algorithmic
determination of commutation relations for Lie symmetry algebras of PDEs.
In {\em Proc. ISSAC~'92 }. ACM Press, Berkeley.

\bibitem{Reid92b}
Reid, G. J. and McKinnon, D. K. (1992b),
{\it Solving systems of linear PDEs in their coefficient field by recursively
decoupling and solving ODEs}, preprint.

\bibitem{Riq10}
Riquier, C. 1910. {\em Les Syst\{e}mes d'\'{E}quations aux D\'{e}riv\'{e}es
Partielles.\/} Gauthier-Villars, Paris.

\bibitem{Ritt50}

Ritt, J. F. (1950), {\it Differential Algebra}, (Dover).

\bibitem{Rust93}

Rust, C. 1993. On The Classification of Rankings of Partial Derivatives.
Preprint.

\bibitem{Schu92}
Sch\"{u}, J., W.M. Seiler, and J. Calmet. 1992. Algorithmic Methods for Lie Pseudogroups.
In {\em Proc. Modern Group Analysis.\/} Acireale.
\bibitem{Sch84}
Schwarz, F. 1984. The Riquier-Janet theory and its application to nonlinear
evolution equations.  {\em Physica\/} {\bf D11}: 243--251.

\bibitem{Sch85}
Schwarz, F. 1985. Automatically determining symmetries of differential equations.
{\em Computing.} {\bf 34}: 91--106.

\bibitem{Sch92a}
Schwarz, F. 1992. An Algorithm for Determining
the Size of Symmetry Groups.{\em Computing.} {\bf 49}: 95--115.

\bibitem{Sch92b}
Schwarz, F. 1992. Reduction and completion algorithms for Partial Differential Equations.
In {\em Proc. ISSAC '92}. 49--56. ACM Press, Berkeley.

\bibitem{Spenc69}
Spencer, D. C. 1969.  {\it Overdetermined systems of linear partial differential equations}. Bull. of the Amer. Math. Soc. {\bf 75}, 179-239.

\bibitem{Thom29}
Thomas, J.M. 1929. Riquier's existence theorems. {\em Annals of Maths.\/}
{\bf 30}: 285--310.

\bibitem{Top}
Topunov, L. 1989. Reducing systems of linear differential equations to a passive form.
{\em Acta Appl. Math.\/} {\bf 16}: 191--206.

\bibitem{Tress94}
Tresse, A. 1894. Sur les invariants diff\'{e}rentiels des groupes continus de
transformations. {\em Acta Mathematica\/} {\bf 18}: 1--88.
(English translation I. Lisle 1989, available from the author)

\bibitem{Wolf89}
Wolf, T. (1989), {\it Proc. EUROCAL 87, Lect. Notes in Comp. Sci.} {\bf 378}, 479.

\bibitem{Wolf91}
Wolf, T. (1991), The symbolic integration of exact PDEs," preprint.

\bibitem{Wolf92}
Wolf, T. (1992), CRACK Reference Manual", REDUCE Network Library.

\bibitem{Weisp93}
Weispfenning, V. (1993), Differential-Term Orders'',
to appear in {\it Proc. of ISAAC 1993} (Kiev, acm press).

\end{thebibliography}
\end{document}

INDEX OF FUNCTIONS AND OPTIONS

This is a preliminary effort at making an index
(it needs to be completed).
To use this index, use your editor's search capability to
find occurrences of the given expression.

FUNCTIONS

standard_form
initial_data
taylor_exp
finite_generators
lie_symm
polyd
Clearderivdep
sfdiff
make_vset
shortenpiv
map2rep
rep2map
man2rep
rep2man
transform
max2rep

OPTIONS

divpiv
class_eqns
usereasy
strategy
storeall
store
little
big
biglittle

PACKAGES

sfprogs
convprogs
maxconv
graphDE
symmetry
liesymm
symgrp.max

OTHERS

_pivs
_incons

`
Special Purpose Systems

webmaster@can.nl