reduce-historical/r34/lib/fide.doc

                             F  I  D  E
                             ==========




                  A REDUCE package for automation of

                 FInite difference method for partial

                     Differential Equation solving






                              Version 1.1




                             User's Manual
                             -------------




                              Richard Liska




          Faculty of Nuclear Science and Physical Engineering
                    Technical University of Prague
              Brehova 7, 115 19 Prague 1, Czechoslovakia
                     E-mail: tjerl@cspuni12.bitnet (EARN)
                     Fax: (42 - 2) 84 73 54
                     Tel: (42 - 2) 84 77 86



                                May 1991




                                         1









                               Abstract
                               --------


     The FIDE  package performs  automation of  the process of numerical
solving  partial  differential  equations  systems  (PDES)  by  means of
computer algebra.  For PDES solving finite difference method is applied.
The  computer  algebra  system  REDUCE  and  the  numerical  programming
language FORTRAN  are used in the presented methodology. The main aim of
this methodology is to  speed  up  the  process  of  preparing numerical
programs for  solving PDES.  This process is quite often, especially for
complicated systems, a tedious and time consuming task.
     In the process  one  can  find  several  stages  in  which computer
algebra  can  be  used  for  performing routine analytical calculations,
namely: transforming differential  equations  into  different coordinate
systems,   discretization   of   differential   equations,  analysis  of
difference  schemes  and  generation  of  numerical  programs.  The FIDE
package consists of the following modules:

  EXPRES  for transforming PDES into any orthogonal coordinate system.
  IIMET   for discretization of PDES by integro-interpolation method.
  APPROX  for determining the order of approximation of difference
          scheme.
  CHARPOL for calculation of amplification matrix and characteristic
          polynomial of difference scheme, which are needed in Fourier
          stability analysis.
  HURWP   for polynomial roots locating necessary in verifying the von
          Neumann stability condition.
  LINBAND for generating the block of FORTRAN code, which solves a
          system of linear algebraic equations with band matrix
          appearing quite often in difference schemes.

     Version  1.1  of  the  FIDE  package  is the result of porting FIDE
package to REDUCE 3.4. In comparison with Version 1.0 some features  has
been changed  in the  LINBAND module  (possibility to  interface several
numerical libraries).





References
----------

[1] R. Liska, L. Drska: FIDE: A REDUCE package for automation of  FInite
      difference method for  solving pDE. In  ISSAC '90, Proceedings  of
      the International Symposium on Symbolic and Algebraic Computation,
      Ed. S. Watanabe, M. Nagata. p. 169-176, ACM Press, Addison Wesley,
      New York 1990.


                                         2



                    Table of contents
                    =================

1   E X P R E S                                                      4

1.1 The specification of the coordinate system.......................4
1.2 The declaration of tensor quantities.............................5
1.3 New infix operators..............................................5
1.4 New prefix operators.............................................5
1.5 Tensor expressions...............................................6
1.6 Assigning statement..............................................7

2   I I M E T                                                        8

2.1 Specification of the coordinates and the indices
    corresponding to them............................................8
2.2 Difference grids.................................................9
2.3 Declaring the dependence of functions on coordinates............10
2.4 Functions and difference grids..................................11
2.5 Equations and difference grids..................................12
2.6 Discretization of basic terms...................................13
2.7 Discretization of a system of equations.........................18
2.8 Error messages..................................................20

3   A P P R O X                                                     21

3.1 Specification of the coordinates and the indices
    corresponding to them...........................................21
3.2 Specification of the Taylor expansion...........................21
3.3 Function declaration............................................22
3.4 Order of accuracy determination.................................23

4   C H A R P O L                                                   24

4.1 Commands common with the IIMET module...........................24
4.2 Function declaration............................................24
4.3 Amplification matrix............................................25
4.4 Characteristic polynomial.......................................26
4.5 Automatic denotation............................................26

5   H U R W P                                                       28

5.1 Conformal mapping...............................................28
5.2 Investigation of polynomial roots...............................28

6   L I N B A N D                                                   30

6.1 Program generation..............................................30
6.2 Choosing the numerical library..................................32
6.3 Completion of the generated code................................32






                                         3















                            1   E X P R E S
                                ===========



                A Module for Transforming Differential
         Operators and Equations into an Arbitrary Orthogonal
                           Coordinate System



     This module  makes it  possible to  express various scalar, vector,
and tensor differential equations  in any  orthogonal coordinate system.
All transformations  needed are  executed automatically according to the
coordinate  system  given  by  the  user.  The  module  was  implemented
according to the similar MACSYMA module from [1].


1.1 The specification of the coordinate system
    ------------------------------------------

     The coordinate system is specified using the following statement:

  SCALEFACTORS <d>,<tr 1>,...,<tr d>,<cor 1>,...,<cor d>;
  <d> ::= 2 | 3                       coordinate system dimension
  <tr i> ::= "algebraic expression"   the expression of the i-th
                                      Cartesian coordinate in new
                                      coordinates
  <cor i> ::= "identifier"            the i-th new coordinate

All evaluated  quantities are transformed into the coordinate system set
by the last SCALEFACTORS statement. By default, if this statement is not
applied, the  three-dimensional Cartesian coordinate system is employed.
During the evaluation of SCALEFACTORS statement the metric coefficients,
i.e. scale  factors SF(i),  of a  defined coordinate system are computed
and printed. If the WRCHRI switch  is ON,  then the  nonzero Christoffel
symbols of the coordinate system are printed too. By  default the WRCHRI
switch is OFF.






                                         4






1.2 The declaration of tensor quantities
    ------------------------------------

     Tensor  quantities  are  represented  by  identifiers.  The VECTORS
declaration declares  the identifiers  as vectors, the DYADS declaration
declares the identifiers as dyads. i.e. two-dimensional tensors, and the
TENSOR  declaration  declares  the  identifiers as tensor variables. The
declarations have the following syntax:

  <declaration> <id 1>,<id 2>,...,<id n>;
  <declaration> ::= VECTORS | DYADS | TENSOR
  <id i> ::= "identifier"

The value of the identifier V declared as vector in  the two-dimensional
coordinate  system  is  (V(1),  V(2)),  where V(i) are the components of
vector V. The value of the identifier T declared as a dyad  is ((T(1,1),
T(1,2)), (T(2,1),  T(2,2))). The value of the tensor variable can be any
tensor (see below). Tensor variables  can  be  used  only  for  a single
coordinate  system,  after  the  coordinate  system  redefining by a new
SCALEFACTORS statement, the tensor variables have to be re-defined using
the assigning statement.


1.3 New infix operators
    -------------------

     For  four  different  products  between  the tensor quantities, new
infix operators have been  introduced  (in  the  explaining  examples, a
two-dimensional  coordinate  system,  vectors  U,  V, and dyads T, W are
considered):

  .  - scalar product                U.V = U(1)*V(1)+U(2)*V(2)
  ?  - vector product                U?V = U(1)*V(2)-U(2)*V(1)
  &  - outer product                 U&V = ((U(1)*V(1),U(1)*V(2)),
                                            (U(2)*V(1),U(2)*V(2)))
  #  - double scalar product         T#W = T(1,1)*W(1,1)+T(1,2)*W(1,2)+
                                           T(2,1)*W(2,1)+T(2,2)*W(2,2)

The other usual arithmetic infix operators +, -,  *, **  can be  used in
all situations  that have  sense (e.g. vector addition, a multiplication
of a tensor by a scalar, etc.).


1.4 New prefix operators
    --------------------

     New  prefix  operators  have  been  introduced  to  express  tensor
quantities  in  its  components  and the differential operators over the
tensor quantities:

  VECT  - the explicit expression of a vector in its components
  DYAD  - the explicit expression of a dyad in its components

                                         5






  GRAD  - differential operator of gradient
  DIV   - differential operator of divergence
  LAPL  - Laplace's differential operator
  CURL  - differential operator of curl
  DIRDF - differential operator of the derivative in direction
          (1st argument is the directional vector)

The results  of the  differential operators  are written  using the DIFF
operator.  DIFF(<scalar>,<cor  i>)  expresses the derivative of <scalar>
with respect to the  coordinate <cor  i>. This  operator is  not further
simplified. If  the user wants to make it simpler as common derivatives,
he performs the following declaration:

  FOR ALL X,Y LET DIFF(X,Y) = DF(X,Y);  .

Then, however, we must realize that if the scalars or  tensor quantities
do not directly explicitly depend on the coordinates, their dependencies
have  to  be  declared  using  the  DEPEND  statements,   otherwise  the
derivative will  be evaluated  to zero.  The dependence of all vector or
dyadic components (as dependence of the name of  vector or  dyad) has to
appear  before  VECTORS  or  DYADS  declarations,  otherwise after these
declarations one has to declare the dependencies of  all components. For
formulating  the   explicit  derivatives   of  tensor  expressions,  the
differentiation operator DF can be used  (e.g. the  differentiation of a
vector in its components).


1.5 Tensor expressions
    ------------------

     Tensor expressions  are the  input into  the EXPRES  module and can
have a variety of forms. The output is then the formulation of the given
tensor expression  in the  specified coordinate system. The most general
form of a  tensor  expression  <tensor>  is  described  as  follows (the
conditions  (d=i)  represent  the  limitation  on  the  dimension of the
coordinate system equalling i):

  <tensor> ::= <scalar> | <vector> | <dyad>
  <scalar> ::= "algebraic expression, can contain <scalars>" |
               "tensor variable with scalar value" |
               <vector 1>.<vector 2> | <dyad 1>#<dyad 2> |
               (d=2)<vector 1>?<vector 2> | DIV <vector> |
               LAPL <scalar> | (d=2) ROT <vector> |
               DIRDF(<vector>,<scalar>)
  <vector> ::= "identifier declared by VECTORS statement" |
               "tensor variable with vector value" |
               VECT(<scalar 1>,...,<scalar d>) | -<vector> |
               <vector 1>+<vector 2> | <vector 1>-<vector 2> |
               <scalar>*<vector> | <vector>/<scalar> |
               <dyad>.<vector> | <vector>.<dyad> | (d=3)
               <vector 1>?<vector 2> | (d=2) <vector>?<dyad> |
               (d=2) <dyad>?<vector> | GRAD <scalar> |

                                         6






               DIV <dyad> | LAPL <vector> | (d=3) ROT <vector> |
               DIRDF(<vector 1>,<vector 2>) | DF(<vector>,"usual
               further arguments")
  <dyad>    ::= "identifier declared by DYADS statement" |
               "tensor variable with dyadic value" |
               DYAD((<scalar 11>,...,<scalar 1d>),...,(<scalar d1>,
               ...,<scalar dd>)) | -<dyad> | <dyad 1>+<dyad 2> |
               <dyad 1>-<dyad 2> | <scalar>*<dyad> | <dyad>/<scalar>
               | <dyad 1>.<dyad 2> | <vector 1>&<vector 2> |
               (d=3) <vector>?<dyad> | (d=3) <dyad>?<vector> |
               GRAD <vector> | DF(<dyad>,"usual further arguments")


1.6 Assigning statement
    -------------------

     The assigning statement for  tensor variables  has a  usual syntax,
namely:

  <tensor variable> := <tensor>
  <tensor variable> ::= "identifier declared TENSOR"  .

The assigning  statement assigns  the tensor  variable the  value of the
given tensor expression,  formulated  in  the  given  coordinate system.
After a change of the coordinate system, the tensor variables have to be
redefined.



References
----------

[1] M. C. Wirth, On the Automation of Computational Physics. PhDr
          Thesis. Report UCRL-52996, Lawrence Livermore National
          Laboratory, Livermore, 1980.


















                                         7















                             2   I I M E T
                                 =========



                 A Module for Discretizing the Systems
                   of Partial Differential Equations



     This program module makes it possible  to discretize  the specified
system of partial differential equations using the integro-interpolation
method,  minimizing  the  number  of  the  used  interpolations  in each
independent variable.  It can  be used for non-linear systems and vector
or tensor variables as well.  The user specifies the way of discretizing
individual terms  of differential equations, controls the discretization
and obtains various difference schemes according to his own wish.


2.1 Specification of the coordinates and the indices corresponding
    to them
    --------------------------------------------------------------

     The independent variables of differential equations  will be called
coordinates.  The  names  of  the  coordinates and the indices that will
correspond to the particular  coordinates in  the difference  scheme are
defined using the COORDINATES statement:

  COORDINATES  <coordinate 1>{,<coordinate i>} [ INTO
             <index 1>{,<index i>}];
  <coordinate i> ::= "identifier"  - the name of the coordinate
  <index i> ::= "identifier"       - the name of the index

This statement  specifies that the <coordinate i> will correspond to the
<index i>. A new COORDINATES statement cancels the  definitions given by
the preceding  COORDINATES statement.  If the  part [  INTO ... ] is not
included in the statement, the  statement  assigns  the  coordinates the
indices I, J, K, L, M, N, respectively. If it is included, the number of
coordinates and the number of indices should be the same.





                                         8






2.2 Difference grids
    ----------------

     In the discretization, orthogonal  difference  grids  are employed.
In addition to the basic grid, called the integer one, there is another,
the half-integer grid in each coordinate, whose cellular boundary points
lie in  the centers of the cells of the integer grid. The designation of
the  cellular  separating  points  and  centers  is  determined  by  the
CENTERGRID switch:  if it is ON and the index in the given coordinate is
I, the centers of the grid cells are designated by indices I, I + 1,...,
and the  boundary points of the cells by indices I + 1/2,..., if, on the
contrary, the switch is OFF,  the  cellular  centers  are  designated by
indices I  + 1/2,...,  and the  boundary points  by indices I, I + 1,...
(see Fig. 2.1).


                                         ON CENTERGRID
  I-1/2     I       I+1/2          I+1            I+3/2
---|--------|--------|--------------|--------------|----
   I       I+1/2    I+1            I+3/2          I+2
                                         OFF CENTERGRID

              Figure 2.1 Types of grid


In the case of  ON CENTERGRID,  the indices  i,i+1,i-1... thus designate
the centers  of the cells of the integer grid and the boundary points of
the cells of the half-integer grid, and, similarly,  in the  case of OFF
CENTERGRID,  the  boundaries  of  the  cells of the integer grid and the
central points of the half-integer grid. The meaning of the  integer and
half-integer  grids  depends  on  the CENTERGRID switch in the described
way. After the package is loaded, the CENTERGRID is ON.  Obviously, this
switch is significant only for non-uniform grids with a variable size of
each cell.
     The grids can be uniform, i.e. with a constant cell size - the step
of the grid. The following statement:

  GRID UNIFORM,<coordinate>{,<coordinate>};

defines  uniform  grids  in  all  coordinates  occurring  in  it.  Those
coordinates that do not occur in the GRID UNIFORM statement are supposed
to have non-uniform grids.
     In the  outputs, the grid step is designated by the identifier that
is made by putting the character  H before  the name  of the coordinate.
For a  uniform grid, this identifier (e.g. for the coordinate X the grid
step HX) has the meaning of a step of an  integer or  half-integer grids
that  are  identical.  For  a  non-uniform  grid,  this identifier is an
operator and has the  meaning of  a step  of an  integer grid,  i.e. the
length  of  a  cell  whose  center  (in  the  case  of ON CENTERGRID) or
beginning (in the case of OFF  CENTERGRID)  is  designated  by  a single
argument  of  this  operator.  For  each  coordinate s designated by the


                                         9






identifier i, this step of the  integer non-uniform  grid is  defined as
follows:

  Hs(i+j) = s(i+j+1/2) - s(i+j-1/2)      at ON CENTERGRID
  Hs(i+j) = s(i+j+1) - s(i+j)            at OFF CENTERGRID

for all integers j (s(k) designates the value of the coordinate s in the
cellular boundary point subscripted with the index k). The steps  of the
half-integer non-uniform grid are not applied in outputs.


2.3 Declaring the dependence of functions on coordinates
    ----------------------------------------------------

     In  the  system  of  partial  differential  equations, two types of
functions, in other words dependent  variables  can  occur:  namely, the
given  functions,  whose  values  are  known  before the given system is
solved, and the sought functions, whose  values are  not available until
the system  of equations is solved. The functions can be scalar, vector,
or tensor, for vector or tensor  functions the  EXPRES module  has to be
applied at  the same  time. The  names of  the functions employed in the
given system and their dependence on the coordinates are specified using
the DEPENDENCE statement.

  DEPENDENCE <dependence>{,<dependence>};
  <dependence> ::= <function>([<order>],<coordinate>{,
                   <coordinate>})
  <function> ::= "identifier"  - the name of the function
  <order> ::= 1|2 tensor order of the function (the value of
              the function is 1 - vector, 2 - dyad (two-
              dimensional tensor))

Every <dependence>  in the  statement determines  on which <coordinates>
the <function> depends. If the tensor <order> of the function  occurs in
the <dependence>,  the <function> is declared as a vector or a dyad. If,
however, the <function> has  been  declared  by  the  VECTORS  and DYADS
statements of  the EXPRES  module, the  user need not present the tensor
<order>. By default, a function without  any declaration  is regarded as
scalar. In the discretization, all scalar components of tensor functions
are replaced by  identifiers  that  arise  by  putting  successively the
function name  and the  individual indices  of the given component (e.g.
the tensor component T(1,2), written in the EXPRES module as  T(1,2), is
represented by  the identifier  T12). Before the DEPENDENCE statement is
executed, the coordinates  have  to  be  defined  using  the COORDINATES
statement. There  may be  several DEPENDENCE  statements. The DEPENDENCE
statement cancels all preceding determinations of which grids  are to be
used for differentiating the function or the equation for this function.
These determinations can be  either  defined  by  the  ISGRID  or GRIDEQ
statements, or computed in the evaluation of the IIM statement.
     The GIVEN statement:

  GIVEN <function>{,<function>};

                                        10







declares all  functions included  in it  as given functions whose values
are known to the user or can be computed. The CLEARGIVEN statement:

  CLEARGIVEN;

cancels all preceding GIVEN declarations. If  the TWOGRID  switch is ON,
the given  functions can  be differentiated  both on the integer and the
half-integer grids. If the TWOGRID switch is OFF, any given function can
be differentiated  only on  one grid.  After the  package is loaded, the
TWOGRID is ON.


2.4 Functions and difference grids
    ------------------------------

     Every scalar function or scalar component of a  vector or  a dyadic
function occurring  in the  discretized system can be discretized in any
of the coordinates either on the  integer or  half-integer grid.  One of
the tasks  of the  IIMET module  is to  find the optimum distribution of
each of these dependent variables  of  the  system  on  the  integer and
half-integer grids  in all variables so that the number of the performed
interpolations in the integro-interpolation method will be minimal.
     Using the statement

  SAME <function>{,<function>};

all functions given in one of these declarations will be  discretized on
the same  grids in all coordinates. In each SAME statement, at least one
of these functions in one SAME statement must be the sought one.  If the
given function occurs in the SAME statement, it will be discretized only
on one grid, regardless of the state of the TWOGRID switch. If  a vector
or a  dyadic function  occurs in  the SAME statement, what has been said
above relates to all  its  scalar  components.  There  are  several SAME
statements that can be presented. All SAME statements can be canceled by
the following statement:

  CLEARSAME;

The SAME statement can be successfully used, for example, when the given
function depends  on the  function sought  in a  complicated manner that
cannot be  included  either  in  the  differential  equation  or  in the
difference scheme explicitly, and when both the functions are desired to
be discretized in the same points so that the user will not be forced to
execute the interpolation during the evaluation of the given function.
     In  some  cases,  it  is  convenient  too to specify directly which
variable on which grid is to be discretized,  for which  case the ISGRID
statement is applied:

  ISGRID <s-function>{,<s-function>};
  <s-function> ::= <function>([<component>,]<s-grid>{,<s-grid>})
  <s-grid> ::= <coordinate> .. <grid>,

                                        11






  <grid> ::= ONE | HALF            designation of the integer
                                   (ONE) and half-integer (HALF)
                                   grids
  <component> ::= <i-dim> |        for the vector <function>
                  <i-dim>,<i-dim>  for the dyadic <function>
                                   it is not presented for the
                                   scalar <function>
  <i-dim> ::= *| "natural number from 1 to the space dimension
                the space dimension is specified in the EXPRES
                module by the SCALEFACTORS statement, * means all
                components

The statement  defines that the given functions or their components will
be discretized in the specified coordinates  on the  specified grids, so
that, for example, the statement ISGRID U (X..ONE,Y..HALF), V(1,Z..ONE),
T(*,1,X..HALF); defines that scalar U will be discretized on the integer
grid in  the coordinate X, and on the half-integer one in the coordinate
Y, the first component of vector V will  be on  the integer  grid in the
coordinate  Z,  and  the  first  column  of  tensor  T  will  be  on the
half-integer grid in the  coordinate  X.  The  ISGRID  statement  can be
applied  more  times.  The  functions  used in this statement have to be
declared before by the DEPENDENCE statement.


2.5 Equations and difference grids
    ------------------------------

     Every equation of the system of  partial differential  equations is
an equation  for some  sought function (specified in the IIM statement).
The correspondence between the sought  functions  and  the  equations is
mutually  unambiguous.   The  GRIDEQ  statement  makes  it  possible  to
determine on which grid  an individual  equation will  be discretized in
some or all coordinates

  GRIDEQ <g-function>{,<g-function>};
  <g-function> ::= <function>(<s-grid>{,<s-grid>})

Every  equation  can  be  discretized  in  any  coordinate either on the
integer   or   half-integer   grid.   This   statement   determines  the
discretization of the equations given by the functions included in it in
given coordinates, on given grids.  The  meaning  of  the  fact  that an
equation is discretized on a certain grid is as follows: index I used in
the DIFMATCH statements (discussed in the following section), specifying
the discretization  of the basic terms, will be located in the center of
the cell of this  grid,  and  indices  I+1/2,  I-1/2  from  the DIFMATCH
statement on the boundaries of the cell of this grid. The actual name of
the index in the  given coordinate  is determined  using the COORDINATES
statement, and its location on the grid is set by the CENTERGRID switch.





                                        12






2.6 Discretization of basic terms
    -----------------------------

     The discretization of a system of partial differential equations is
executed successively in individual  coordinates. In  the discretization
of an  equation in  one coordinate,  the equation is linearized into its
basic terms first that will be discretized independently then.   If D is
the  designation  for  the  discretization operator in the coordinate x,
this linearization obeys the following rules:

  1. D(a+b) = D(a)+D(b)
  2. D(-a)  = -D(a)
  3. D(p.a) = p.D(a)        (p does not depend on the coordinate x)
  4. D(a/p) = D(a)/p

     The linearization lasts as long  as  some  of  these  rules  can be
applied.   The   basic   terms   that  must  be  discretized  after  the
linearization have then the forms of the following quantities:

  1. The actual coordinate in which the discretization is performed.
  2. The sought function.
  3. The given function.
  4. The product of the quantities 1 - 7.
  5. The quotient of the quantities 1 - 7.
  6. The natural power of the quantities 1 - 7.
  7. The derivative of the quantities 1 - 7 with respect to the
     actual coordinate.

The way of discretizing these basic  terms, while  the functions  are on
integer  and  half-integer  grids,  is  determined  using  the  DIFMATCH
statement:

  DIFMATCH <coordinate>,<pattern term>,{{<grid specification>,}
           <number of interpolations>, <discretized term>};
  <coordinate> ::= ALL | "identifier" - the coordinate name from
                                        the COORDINATES statement
  <pattern term> ::= <pattern coordinate>|
       <pattern sought function>|
       <pattern given function>|<pattern term> *
       <pattern term>|<pattern term> / <pattern term>|
       <pattern term> ** <exponent>|
       DIFF(<pattern term>,<pattern coordinate>[,<order
       of derivative>])|
       <declared operator>(<pattern term>{,<pattern term>})
  <pattern coordinate> ::= X
  <pattern sought function> ::= U | V | W
  <pattern given function> ::= F | G
  <exponent> ::= N | "integer greater than 1"
  <order of derivative> ::= "integer greater than 2"
  <grid specification> ::= <pattern function>=<grid>
  <pattern function> ::= <pattern sought function>|
                         <pattern given function>

                                        13






  <number of interpolations> ::= "non-negative integer"
  <discretized term> ::= <pattern operator>(<index expression>)|
                "natural number"|DI|DIM1|DIP1|DIM2|DIP2|
                <declared term> | - <discretized term> |
                <discretized term> + <discretized term> |
                <discretized term> * <discretized term> |
                <discretized term> / <discretized term> |
                (<discretized term>) |
                <discretized term> **<exponent>
  <pattern operator> ::= X | U | V | W | F | G
  <index expression> ::= <pattern index> |
                         <pattern index> + <increment> |
                         <pattern index> - <increment>
  <pattern index> ::= I
  <increment> = "rational number"
  DIFCONST <declared term>{,<declared term>};
  <declared term> ::= "identifier" - the constant parameter of
                                     the difference scheme.
  DIFFUNC <declared operator>{,<declared operator>};
  <declared operator> ::= "identifier" - prefix operator, that can
                appear in discretized equations (e.g. SIN).

The first parameter of the DIFMATCH statement determines  the coordinate
for which the discretization defined in it is valid. If ALL is used, the
discretization  will   be   valid   for   all   coordinates,   and  this
discretization is  accepted when  it has  been checked whether there has
been no other discretization  defined for  the given  coordinate and the
given  pattern  term.  Each  pattern  sought  function, occurring in the
pattern term, must be included in  the specification  of the  grids. The
pattern  given  functions  from  the  pattern term can occur in the grid
specification, but in some cases  (see  below)  need  not.  In  the grid
specification the  maximum number  of 3 pattern functions may occur. The
discretization  of  each  pattern  term  has  to  be  specified  in  all
combinations   of   the   pattern   functions   occurring  in  the  grid
specification, on the  integer  and  half-integer  grids,  that  is 2**n
variants   for   the   grid   specification  with  n  pattern  functions
(n=0,1,2,3). The discretized term is the  discretization of  the pattern
term in  the pattern  coordinate X in the point X(I) on the pattern grid
(see  Fig.  2.2),  and  the  pattern  functions  occurring  in  the grid
specification are  in the  discretized term on the respective grids from
this  specification  (to  the  discretized  term  corresponds  the  grid
specification preceding it).

                                            integer grid
        X(I-1)                X(I)               X(I+1)
          |        DIM1        |       DIP1        |
---|------|------|-------------|-------------|-----|-----|---
   |    DIM2     |            DI             |    DIP2   |
 X(I-3/2)      X(I-1/2)                    X(I+1/2)     X(I+3/2)
                                            half-integer grid

                 Figure 2.2 Pattern grid

                                        14







The pattern grid steps defined as

  DIM2 = X(I - 1/2) - X(I - 3/2)
  DIM1 = X(I) - X(I - 1)
  DI   = X(I + 1/2) - X(I - 1/2)
  DIP1 = X(I + 1) - X(I)
  DIP2 = X(I + 3/2) - X(I + 1/2)

can occur in the discretized term.
     In  the  integro-interpolation  method,  the  discretized  term  is
specified by the integral

  <discretized term>=1/(X(I+1/2)-X(I-1/2))*DINT(X(I-1/2),X(I+1/2),
                                                <pattern term>,X),

where DINT is operator of definite integration DINT(from,  to, function,
variable).   The   number   of   interpolations   determines   how  many
interpolations were needed for  calculating this  integral in  the given
discrete form (the function on the integer or half-integer grid). If the
integro-interpolation method is not used,  the  more  convenient  is the
distribution of the functions on the half-integer and integer grids, the
smaller number is chosen by the user.  The parameters  of the difference
scheme defined  by the  DIFCONST statement  can occur in the discretized
expression  too  (for  example,  the  implicit-explicit  scheme  on  the
implicit layer  multiplied by  the constant C and on the explicit one by
(1-C)).  As a matter of fact, all DIFMATCH statements  create a  base of
pattern  terms  with  the  rules  of  how  to  discretize these terms in
individual coordinates under the assumption that the functions occurring
in  the   pattern  terms  are  on  the  grids  determined  in  the  grid
specification (all combinations must be included).
     The DIFMATCH statement does not check whether the  discretized term
is actually  the discretization  of the  pattern term  or whether in the
discretized term occur the functions from the grid  specification on the
grids  given  by  this  specification.  An  example can be the following
definition of the discretization of the first and  second derivatives of
the sought function in the coordinate R on a uniform grid:

  DIFMATCH R,DIFF(U,X),U=ONE,2,(U(I+1)-U(I-1))/(2*DI);
                         U=HALF,0,(U(I+1/2)-U(I-1/2))/DI;
  DIFMATCH R,DIFF(U,X,2),U=ONE,0,(U(I+1)-2*U(I)+U(I-1))/DI**2,
    U=HALF,2,(U(I+3/2)-U(I+1/2)-U(I-1/2)+U(I-3/2))/(2*DI**2);

     All DIFMATCH statements can be cleared by the statement

  CLEARDIFMATCH;

After this statement user has to supply its own DIFMATCH statements.
     But now back to the discretizing of the basic terms obtained by the
linearization of the partial differential equation, as mentioned  at the
beginning of  this section.  Using the  method of  pattern matching, for
each basic term a term representing its pattern is found in the  base of

                                        15






pattern  terms  (specified  by  the  DIFMATCH  statements).  The pattern
matching obeys the following rules:

  1. The  pattern for  the coordinate  in which  the discretization is
  executed is the pattern coordinate X.

  2.  The  pattern  for  the  sought  function  is some pattern sought
  function, and this correspondence is mutually unambiguous.

  3.  The  pattern  for  the  given  function  is  some  pattern given
  function, or,  in case  the EQFU  switch is  ON, some pattern sought
  function, and, again, the correspondence  of  the  pattern  with the
  given  function  is  mutually  unambiguous  (after  loading the EQFU
  switch is ON).

  4. The pattern for the products of quantities is the  product of the
  patterns of these quantities, irrespective of their sequence.

  5. The pattern for the quotient of quantities is the quotient of the
  patterns of these quantities.

  6. The pattern for the natural power of a quantity is the same power
  of the  pattern of  this quantity or the power of this quantity with
  the pattern exponent N.

  7. The pattern for the derivative of a quantity with  respect to the
  coordinate in which the discretization is executed is the derivative
  of  the  pattern  of  this  quantity  with  respect  to  the pattern
  coordinate X of the same order of differentiation.

  8. The  pattern for  the sum  of the  quantities that  have the same
  pattern with the identical  correspondence of  functions and pattern
  functions is  this common  pattern (so that it will not be necessary
  to multiply the parentheses during discretizing the products  in the
  second and further coordinates).

When  matching  the  pattern  of  one  basic term, the program finds the
pattern term and the functions corresponding  to the  pattern functions,
maybe also  the exponent  corresponding to the pattern exponent N. After
determining on which grids the individual  functions and  the individual
equations  will  be  discretized,  which  will  be discussed in the next
section, the program finds in the pattern term base the discretized term
either with  pattern functions  on the  same grids  as are the functions
from the basic  term  corresponding  to  them  in  case  that  the given
equation  is  differentiated  on  the  integer  grid,  or  with  pattern
functions on inverse grids  (an inverse  integer grid  is a half-integer
grid, and  vice versa)  compared with  those used for the functions from
the basic term corresponding to  them  in  case  the  given  equation is
differentiated  on  the  half-integer  grid (the discretized term in the
DIFMATCH statement is expressed in the point X(I),  i.e. on  the integer
grid,  and  holds  for  the  discretizing of the equation on the integer
grid; with regard to the substitutions for the pattern index I mentioned

                                        16






later, it is possible to proceed in this way and not necessary to define
the discretization in the points X(I+1/2) too, i.e.  on the half-integer
grid). The program replaces in the thus obtained discretized term:

  1.  The  pattern  coordinate  X  with the particular coordinate s in
  which the discretization is actually performed.

  2. The pattern index I and the grid steps DIM2, DIM1, DI, DIP1, DIP2
  with the expression given in table 2.1 according to the state of the
  CENTERGRID switch and to the  fact  whether  the  given  equation is
  discretized  on  the  integer  or  half-integer grid (i is the index
  corresponding to  the  coordinate  s  according  to  the COORDINATES
  statement, the grid steps were defined in section 2.2)

  3. The  pattern functions  with the corresponding functions from the
  basic  term  and,   possibly,   the   pattern   exponent   with  the
  corresponding exponent from the basic term.

--------------------------------------------------------------------
|                   the equation discretized on                    |
|          the integer grid         |    the half-integer grid     |
|        CENTERGRID      |CENTERGRID|CENTERGRID|     CENTERGRID    |
|           OFF          |   ON     |   OFF    |         ON        |
|------------------------------------------------------------------|
| I  |            i                 |            i+1/2             |
|----|-------------------------------------------------------------|
|DIM2|(Hs(i-2)+Hs(i-1))/2|     Hs(i-1)         |(Hs(i-1)+Hs(i))/2  |
|DIM1|      Hs(i-1)      | (Hs(i-1)+Hs(i))/2   |      Hs(i)        |
|DI  |(Hs(i-1)+Hs(i))/2  |     Hs(i)           |(Hs(i)+Hs(i+1))/2  |
|DIP1|      Hs(i)        | (Hs(i)+Hs(i+1))/2   |      Hs(i+1)      |
|DIP2|(Hs(i)+Hs(i+1))/2  |     Hs(i+1)         |(Hs(i+1)+Hs(i+2))/2|
--------------------------------------------------------------------

        Table 2.1  Values of the pattern index and
                   the pattern grid steps.

     More details  will be  given now to the discretization of the given
functions and  its specification.  The given  function may  occur in the
SAME statement, which makes it bound with some sought function, in other
words it can be discretized only on one grid. This means that  all basic
terms, in  which this  function occurs, must have their pattern terms in
whose discretization definitions by  the DIFMATCH  statement the pattern
function corresponding  to the  mentioned given function has to occur in
the grid specification. If the given function does not occur in the SAME
statement and the TWOGRID switch is OFF, i.e. it can be discretized only
on one grid again, the same holds true. If, however,  the given function
does not  occur in the SAME statement and the TWOGRID switch is ON, i.e.
it can be discretized simultaneously on the integer and the half-integer
grids, then  the basic  terms of  the equations  including this function
have their pattern terms in whose discretization definitions the pattern
function corresponding to the mentioned given function need not occur in
the grid specification.  If, however,  in  spite  of  all,  this pattern

                                        17






function  in  the  discretization  definition  does  occur  in  the grid
specification,  it  is  the  alternative  with   a  smaller   number  of
interpolations occurring  in the DIFMATCH statement that is selected for
each particular basic  term  with  a  corresponding  pattern  (the given
function can be on the integer or half-integer grid).
     Before the  discretization is  executed, it  is necessary to define
using the DIFMATCH statements  the discretization  of all  pattern terms
that are  the patterns  of all basic terms of all equations appearing in
the discretized  system in  all coordinates.  The fact  that the pattern
terms  of  the  basic  terms  of  partial  equations occur repeatedly in
individual systems has made it  possible  to  create  a  library  of the
discretizations  of   the  basic   types  of  pattern  terms  using  the
integro-interpolation method. This library  is a  component part  of the
IIMET module (in its end) and makes work easier for those users who find
the  pattern  matching  mechanism  described  here  too  difficult.  New
DIFMATCH statements  have to  be created  by those  whose equations will
contain a basic term  having no  pattern in  this library,  or those who
need  another  method  to  perform  the  discretization.  The  described
implemented algorithm  of discretizing  the basic  terms is sufficiently
general  to  enable  the  use  of  a  nearly arbitrary discretization on
orthogonal grids.


2.7 Discretization of a system of equations
    ---------------------------------------

     All statements influencing the run of  the discretization  that one
want use  in this  run have  to be executed before the discretization is
initiated. The COORDINATES, DEPENDENCE, and DIFMATCH  statements have to
occur  in  all  applications.  Further,  if necessary, the GRID UNIFORM,
GIVEN, ISGRID,  GRIDEQ, SAME,  and DIFCONST  statements can  be used, or
some of  the CENTREGRID,  TWOGRID, EQFU, and FULLEQ switches can be set.
Only  then  the  discretization  of  a  system  of  partial differential
equations can be started using the IIM statement:

  IIM <array>{,<sought function>,<equation>};
  <array> ::= "identifier" - the name of the array for storing
                             the result
  <sought function> ::= "identifier" - the name of the function
                        whose behavior is described by the
                        equation
  <equation> ::= <left side> = <right side>
  <left side> ::= "algebraic expression" , the derivatives are
                  designated by the DIFF operator
  <right side> ::= "algebraic expression"

Hence, in the IIM statement the name of the array in which the resulting
difference schemes will  be  stored,  and  the  pair  sought  function -
equation, which  describes this  function, are specified. The meaning of
the relation  between the  sought function  and its  equation during the
discretization lies in the fact that the sought function is preferred in
its equation so that the interpolation  is  not,  if  possible,  used in

                                        18






discretizing  the  terms  of  this  equation  that  contain  it.  In the
equations, the functions and the coordinates appear as  identifiers. The
identifiers that  have not  been declared as functions by the DEPENDENCE
statement or as coordinates by the COORDINATES statement  are considered
constants independent  of the  coordinates. The  partial derivatives are
expressed by the DIFF operator that has the same syntax  as the standard
differentiation operator  DF.   The functions and the equations can also
have the vector or tensor character. If these  non-scalar quantities are
applied,  the  EXPRES  module  has  to  be  used together with the IIMET
module, and also non-scalar  differential operators  such as  GRAD, DIV,
etc. can be employed.
     The sequence  performed by the program in the discretization can be
briefly summed up in the following items:

  1. If there are  non-scalar functions  or equations  in a  system of
  equations, they  are automatically  converted into scalar quantities
  by means of the EXPRES module.

  2.  In   each  equation,   the  terms   containing  derivatives  are
  transferred to  the left side, and the other terms to the right side
  of the equation.

  3. For each coordinate, with respect to the  sequence in  which they
  occur in the COORDINATES statement, the following is executed:

  a) It  is determined  on which grids all functions and all equations
  in the actual coordinate will be discretized, and simultaneously the
  limits  are  kept  resulting  from  the  ISGRID,  GRIDEQ,  and  SAME
  statements if they were used. Such  a distribution  of functions and
  equations on  the grids is selected among all possible variants that
  ensures the minimum sum of all numbers of the interpolations  of the
  basic terms  (specified by  the DIFMATCH statement) of all equations
  if the FULLEQ switch is ON, or of all left sides of the equations if
  the FULLEQ  switch is  OFF (after  the loading  the FULLEQ switch is
  ON).

  b) The  discretization  itself  is  executed,  as  specified  by the
  DIFMATCH statements.

  4. If the array name is A, then if there is only one scalar equation
  in the IIM statement, the discretized left side of this  equation is
  stored in  A(0) and  the discretized  right side  in A(1) (after the
  transfer mentioned in item  2), if  there are  more scalar equations
  than one  in the  IIM statement, the discretization of the left side
  of  the  i-th  scalar   equation  is   stored  in   A(i,0)  and  the
  discretization of the right side in A(i,1).

The IIM  statement can  be used  more times  during one program run, and
between its calls, the discretizing process  can be  altered using other
statements of this module.



                                        19






2.8 Error messages
    --------------

     The IIMET  module provides error messages in the case of the user's
errors. Similarly as in the REDUCE system, the error reporting is marked
with five  stars :  "*****" on  the line  start. Some error messages are
identical with those of  the REDUCE  system. Here  are given  some other
error messages that require a more detailed explanation:

***** Matching of X term not found
        - the discretization of the pattern term that is the pattern of
          the basic term printed on the place X has not been
          defined (using the DIFMATCH statement)
***** Variable of type F not defined on grids in DIFMATCH
        - in the definition of the discretizing of the pattern term
          the given functions were not used in the grid
          specification and are needed now
***** X Free vars not yet implemented
        - in the grid specification in the DIFMATCH statement
          more than 3 pattern functions were used
***** All grids not given for term X
        - in the definition of the discretization of the pattern of
          the basic term printed on the place X not all
          necessary combinations of the grid specification
          of the pattern functions were presented




























                                        20















                            3   A P P R O X
                                ===========




             A Module for Determining the Precision Order
                       of the Difference Scheme



     This  module  makes  it  possible  to  determine  the  differential
equation that is solved by the given difference scheme, and to determine
the order  of accuracy  of the solution of this scheme in the grid steps
in individual coordinates. The  discrete  function  values  are expanded
into the Taylor series in the specified point.


3.1 Specification of the coordinates and the indices
    corresponding to them
    ------------------------------------------------

     The COORDINATES  statement, described  in the  IIMET module manual,
specifying the coordinates and the indices corresponding to  them is the
same  for  this  program  module  as  well.  It has the same meaning and
syntax. The  present  module  version  assumes  a  uniform  grid  in all
coordinates. The  grid step  in the  input difference  schemes has to be
designated by an identifier consisting of the character  H and  the name
of the coordinate, e.g. the step of the coordinate X is HX.


3.2 Specification of the Taylor expansion
    -------------------------------------

     In the  determining of the approximation order, all discrete values
of the functions are expanded into the Taylor series in all coordinates.
In order  to determine  the Taylor  expansion, the program needs to know
the point in which it performs this expansion,  and the  number of terms
in the Taylor series in individual coordinates. The center of the Taylor
expansion is specified by the CENTER statement and  the number  of terms
in  the   Taylor  series  in  individual  coordinates  by  the  MAXORDER
statement:


                                        21






  CENTER <center>{,<center>};
  <center> ::= <coordinate> = <increment>
  <increment> ::= "rational number"
  MAXORDER <order>{,<order>};
  <order> ::= <coordinate> = <number of terms>
  <number of terms> ::= "natural number"

The increment in the CENTER statement determines that the center  of the
Taylor expansion  in the given coordinate will be in the point specified
by the index I + <increment>, where I is the index corresponding to this
coordinate, defined  using the COORDINATES statement, e.g. the following
example

  COORDINATE T,X INTO N,J;
  CENTER T = 1/2, X = 1;
  MAXORDER T = 2, X = 3;

specifies that the center of the Taylor expansion  will be  in the point
(t(n+1/2),x(j+1)) and  that until the second derivatives with respect to
t (second powers of ht) and until the third derivatives  with respect to
x (third  powers of  hx) the expansion will be performed. The CENTER and
MAXORDER statements can be placed only after the  COORDINATES statement.
If the center of the Taylor expansion is not defined in some coordinate,
it is supposed to be in the point given by the index  of this coordinate
(i.e.  zero  increment).  If  the  number  of  the  terms  of the Taylor
expansion is not defined in some coordinate, the  expansion is performed
until the third derivatives with respect to this coordinate.


3.3 Function declaration
    --------------------

     All functions  whose discrete  values are  to be  expanded into the
Taylor series must be declared using the FUNCTIONS statement:

  FUNCTIONS <name of function>{,<name of function>};
  <name of function> ::= "identifier"

In the specification of the difference scheme, the functions are used as
operators with one or more arguments, designating the discrete values of
the functions. Each argument is the  sum of  the coordinate  index (from
the  COORDINATES  statement)  and  a  rational  number. If some index is
omitted in  the  arguments  of  a  function,  this  functional  value is
supposed to lie in the point in which the Taylor expansion is performed,
as specified by the CENTER statement. In other words, if the COORDINATES
and CENTER statements, shown in the example in the previous section, are
valid, then it holds that U(N+1) = U(N+1,J+1) and U(J-1) = U(N+1/2,J-1).
The  FUNCTIONS  statement  can  declare  both  the  sought and the known
functions for the expansion.




                                        22






3.4 Order of accuracy determination
    -------------------------------

     The order of accuracy of the difference scheme is determined by the
APPROX statement:

  APPROX (<diff. scheme>);
  <diff. scheme> ::=  <l. side> = <r. side>
  <l. (r.) side> ::= "algebraic expression"

In the  difference scheme  occur the  functions in the form described in
the  preceding  section,  the  coordinate  indices  and  the  grid steps
described  in  section  3.1,  and  the  other symbolic parameters of the
difference scheme. The APPROX statement expands  all discrete  values of
the functions declared in the FUNCTIONS statement into the Taylor series
in all coordinates (the point in which the Taylor expansion is performed
is specified  by the  CENTER statement,  and the number of the expansion
terms by the MAXORDER  statement), substitutes  the expansions  into the
difference  scheme,  which  gives  a modified differential equation. The
modified differential equation, containing the  grid  steps  too,  is an
equation that  is really solved by the difference scheme (into the given
orders in the grid steps).
     The partial differential equation,  whose solution  is approximated
by the  difference scheme,  is determined by replacing the grid steps by
zeros and is displayed after the following message:

  "Difference scheme approximates differential equation"

Then the following message is displayed:

  "with orders of approximation:"

and the lowest powers  (except  for  zero)  of  the  grid  steps  in all
coordinates,  occurring   in  the  modified  differential  equation  are
written. If the PRAPPROX  switch is  ON, then  the rest  of the modified
differential equation is printed. If this rest is added to the left hand
side of the  approximated  differential  equation,  one  obtain modified
equation. By  default the  PRAPPROX switch is OFF. If the grid steps are
found in some denominator in the modified equation, i.e. with a negative
exponent, the  following message  is written, preceding the approximated
differential equation:

  "Reformulate difference scheme, grid steps remain in denominator"

and the approximated differential  equation is  not correctly determined
(one of  its sides is zero). Generally, this message means that there is
a term in the difference scheme that is not a  difference replacement of
the  derivative,  i.e.  the  ratio  of  the  differences of the discrete
function values and the discrete values of the coordinates (the steps of
the difference grid). The user, however, must realize that in some cases
such a term occurs purposefully in  the difference  scheme (e.g.  on the
grid boundary to keep the scheme conservative).

                                        23















                           4   C H A R P O L
                               =============



           A Module for Calculating the Amplification Matrix
               and the Characteristic Polynomial of the
                           Difference Scheme



     This program  module is  used for  the first  step of the stability
analysis  of  the  difference  scheme  using  the  Fourier   method.  It
substitutes   the   Fourier   components  into  the  difference  scheme,
calculates the amplification matrix  of the  scheme for  transition from
one time layer to another, and computes the characteristic polynomial of
this matrix.


4.1 Commands common with the IIMET module
    -------------------------------------

     The COORDINATES and GRID UNIFORM statements, described in the IIMET
module  manual,  are  applied  in  this  module as well, having the same
meaning and syntax. The time coordinate is assumed  to be  designated by
the identifier T. The present module version requires all coordinates to
have uniform grids, i.e. to be declared in  the GRID  UNIFORM statement.
The grid  step in  the input  difference schemes has to be designated by
the identifier consisting  of  the  character  H  and  the  name  of the
coordinate, e.g. the step of the time coordinate T is HT.


4.2 Function declaration
    --------------------

     The  UNFUNC  statement  declares  the names of the sought functions
used in the difference scheme:

  UNFUNC <function>{,<function>}
  <function> ::= "identifier" - the name of the sought function

The functions are used in the difference schemes  as operators  with one
or  more  arguments  for  designating the discrete function values. Each

                                        24






argument is the sum of the index (from the COORDINATES  statement) and a
rational number.  If some  index is  omitted in  the function arguments,
this function value is supposed to  lie in  the point  specified only by
this index,  which means that, with the indices N and J and the function
U, it holds that U(N+1) =  U(N+1,J) and  U(J-1) =  U(N,J-1). As two-step
(in time)  difference schemes may be used only, the time index may occur
either completely alone in the arguments, or in the sum with a one.


4.3 Amplification matrix
    --------------------

     The AMPMAT matrix operator computes the  amplification matrix  of a
two-step difference  scheme. Its argument is an one column matrix of the
dimension  (1,k),  where  k  is  the  number  of  the  equations  of the
difference scheme, that contains the difference equations of this scheme
as algebraic expressions equal to the difference of  the right  and left
sides  of  the  difference  equations.  The  value  of the AMPMAT matrix
operator is the square  amplification  matrix  of  the  dimension (k,k).
During the  computation of the amplification matrix, two new identifiers
are created for each spatial coordinate. The identifier  made up  of the
character K and the name of the coordinate represents the wave number in
this coordinate, and the identifier made up of  the character  A and the
name of  the coordinate  represents the  product of this wave number and
the grid step in this coordinate divided by the least common multiple of
all  denominators  occurring  in  the  scheme  in  the function argument
containing the index of this coordinate.  On the  output an  equation is
displayed defining the latter identifier. For example, if in the case of
function U and index J in the coordinate  X the  expression U(J+1/2) has
been used in the scheme (and, simultaneously, no denominator higher than
2 has occurred in the  arguments  with  J),  the  following  equation is
displayed: AX: = (KX*HX)/2. The definition of these quantities As allows
to express every sum occurring in  the argument  of the  exponentials as
the  sum  of  these  quantities  multiplied by integers, so that after a
transformation, the amplification matrix  will contain  only sin(As) and
cos(As) (for  all spatial  coordinates s).  The AMPMAT operator performs
these transformations  automatically.  If  the  PRFOURMAT  switch  is ON
(after the  loading it is ON), the matrices H0 and H1 (the amplification
matrix is equal to -H1**(-1)*H0) are displayed during  the evaluation of
the AMPMAT  operator. These  matrices can be used for finding a suitable
substitution for the goniometric functions in the next run for a greater
simplification.
     The  TCON  matrix  operator  transforms  the  square  matrix into a
Hermit-conjugate matrix,  i.e. a  transposed and  complex conjugate one.
Its argument  is the  square matrix  and its  value is  Hermit-conjugate
matrix of the argument. The Hermit-conjugate matrix is  used for testing
the  normality   and  unitarity  of  the  amplification  matrix  in  the
determining of the sufficient stability condition.





                                        25






4.4 Characteristic polynomial
    -------------------------

     The CHARPOL operator calculates  the  characteristic  polynomial of
the given  square matrix.  The variable of the characteristic polynomial
is designated by the LAM identifier. The operator has one  argument, the
square matrix, and its value is its characteristic polynomial in LAM.


4.5 Automatic denotation
    --------------------

     Several  statements  and  procedures  are  designed  for  automatic
denotation of some parts of algebraic  expressions by  identifiers. This
denotation is namely useful when we obtain very large expressions, which
cannot fit into the available  memory.  We  can  denote  subparts  of an
expression from the previous step of calculation by identifiers, replace
these  subparts   by  these   identifiers  and   continue  the  analytic
calculation  only  with  these  identifiers.  Every  time  we  use  this
technique we have to explicitly survive  in processed  expressions those
algebraic quantities  which will  be necessary in the following steps of
calculation. The process  of  denotation  and  replacement  is performed
automatically and  the algebraic  values which  are denoted by these new
identifiers can be  written  out  at  any  time.  We  describe  how this
automatic denotation can be used.
     The  statement  DENOTID  defines  the  beginning  letters  of newly
created identifiers. Its syntax is

  DENOTID <id>;
  <id> ::= "identifier"

After this  statement  the  new  identifiers  created  by  the operators
DENOTEPOL and  DENOTEMAT will  begin with  the letters of the identifier
<id> used in this statement. Without using any DENOTID statement all new
identifiers  will  begin  with  one  letter  A.  We  suggest to use this
statement every time before using operators DENOTEPOL or  DENOTEMAT with
some new  identifier and to choose identifiers used in this statement in
such a way that the newly  created  identifiers  are  not  equal  to any
identifiers used in the expressions you are working with.
     The operator  DENOTEPOL has  one argument, a polynomial in LAM, and
denotes  the  real  and  imaginary  part  of  its  coefficients  by  new
identifiers. The  real part of the j-th LAM power coefficient is denoted
by the identifier <id>R0j and the imaginary part by <id>I0j,  where <id>
is the  identifier used in the last DENOTID statement. The denotation is
done only for non-numeric  coefficients. The  value of  this operator is
the  polynomial  in  LAM  with  coefficients  constructed  from  the new
identifiers.  The  algebraic  expressions  which  are  denoted  by these
identifiers are  stored as  LISP data structure standard quotient in the
LISP variable DENOTATION!* (assoc. list).
     The operator DENOTEMAT has one argument, a matrix,  and denotes the
real and  imaginary parts  of its  elements. The  real part of the (j,k)
matrix element is denoted  by the  identifier <id>Rjk  and the imaginary

                                        26






part  by  <id>Ijk.  The  returned  value of the operator is the original
matrix with non-numeric elements replaced by <id>Rjk +  I*<id>Ijk. Other
matters are the same as for the DENOTEPOL operator.
     The statement PRDENOT has the syntax

  PRDENOT;

and writes  from the  variable DENOTATION!*  the definitions  of all new
identifiers introduced  by the  DENOTEPOL and  DENOTEMAT operators since
the last  call of  CLEARDENOT statement (or program start) in the format
defined by  the  present  setting  of  output  control  declarations and
switches. The  definitions are  written   in the same order as they have
been  entered,  so  that  the  definitions  of  the  first  DENOTEPOL or
DENOTEMAT operators  are written  first. This order guarantees that this
statement can be utilized  directly to  generate a  semantically correct
numerical program  (the identifiers from the first denotation can appear
in the second one, etc.).
     The statement CLEARDENOT with the syntax

  CLEARDENOT;

clears the variable DENOTATION!*, so that all denotations  saved earlier
by the  DENOTEPOL and DENOTEMAT operators in this variable are lost. The
PRDENOT statement succeeding this statement writes nothing.





























                                        27















                             5   H U R W P
                                 =========



                A Module for Polynomial Roots Locating



     This module is used  for verifying  the stability  of a polynomial,
i.e. for  verifying if  all roots  of a  polynomial lie in a unit circle
with its center  in  the  origin.  By  investigating  the characteristic
polynomial  of  the  difference  scheme,  the  user  can  determine  the
conditions of the stability of this scheme.


5.1 Conformal mapping
    -----------------

     The HURW  operator  transforms  a  polynomial  using  the conformal
mapping LAM=(z+1)/(z-1).  Its argument  is a  polynomial in  LAM and its
value is a transformed polynomial in LAM (LAM=z).  If P is  a polynomial
in LAM,  then it holds: all roots LAM1i of the polynomial P are in their
absolute values smaller than one, i.e. |LAM1i|<1, iff the real  parts of
all  roots  LAM2i  of  the  HURW(P)  polynomial  are  negative,  i.e. Re
(LAM2i)<0.
     The elimination of the unit polynomial roots (LAM=1),  which has to
occur before  the conformal  transformation is performed, is made by the
TROOT1 operator. The argument of this  operator is  a polynomial  in LAM
and its  value is  a polynomial  in LAM not having its root equal to one
any more. Mostly, the investigated polynomial has some  more parameters.
For some  special values  of those parameters, the polynomial may have a
unit root.  During the evaluation of the TROOT1 operator,  the condition
concerning  the  polynomial  parameters  is  displayed,  and  if  it  is
fulfilled, the resulting polynomial has a unit root.


5.2 Investigation of polynomial roots
    ---------------------------------

     The HURWITZP operator checks  whether a  polynomial is  the Hurwitz
polynomial, i.e.  whether all  its roots  have negative  real parts. The
argument of the HURWITZP operator is  a polynomial  in LAM  with real or

                                        28






complex  coefficients,  and  its  value  is  YES  if the argument is the
Hurwitz polynomial.  It  is  NO  if  the  argument  is  not  the Hurwitz
polynomial, and COND if it is the Hurwitz polynomial when the conditions
displayed by the HURWITZP  operator during  its analysis  are fulfilled.
These  conditions  have  the  form of inequalities and contain algebraic
expressions made up of the polynomial coefficients. The  conditions have
to  be  valid  either  simultaneously,  or  they  are  designated  and a
proposition is created from them by the AND and OR  logic operators that
has  to  be  fulfilled  (it  is  the condition concerning the parameters
occurring in the polynomial  coefficient)  by  a  polynomial  to  be the
Hurwitz one. This proposition is the sufficient condition, the necessary
condition is the fulfillment of all the inequalities displayed.
     If the HURWITZP  operator  is  called  interactively,  the  user is
directly  asked  if  the  inequalities  are  or  are not valid. The user
responds "Y" if the displayed inequality is valid, "N" if it is not, and
"?" if he does not know whether the inequality is true or not.





































                                        29















                           6   L I N B A N D
                               =============



            A Module for Generating the Numeric Program for
            Solving a System of Linear Algebraic Equations
                           with Band Matrix



     The LINBAND  module generates  the numeric  program in  the FORTRAN
language, which solves a system of linear algebraic equations with  band
matrix using  the routine  from the  LINPACK, NAG  ,IMSL or ESSL program
library.  As  input data only  the system of  equations is given  to the
program.   Automatically,  the  statements  of  the FORTRAN language are
generated that fill the band  matrix of the system in  the corresponding
memory mode of chosen library, call the solving routine, and assign  the
chosen variables to the solution of  the system. The module can be  used
for solving linear difference schemes often having the band matrix.


6.1 Program generation
    ------------------

     The  program   in  the   FORTRAN  language   is  generated  by  the
GENLINBANDSOL statement (the  braces  in  this  syntax  definition occur
directly  in  the  program  and  do  not  have  the usual meaning of the
possibility of repetition, they designate REDUCE lists):

  GENLINBANDSOL (<n-lower>,<n-upper>,{<system>});
  <n-lower> ::= "natural number"
  <n-upper> ::= "natural number"
  <system> ::= <part of system> | <part of system>,<system>
  <part of system>::= {<variable>,<equation>} | <loop>
  <variable> ::= "kernel"
  <equation> ::= <left side> = <right side>
  <left side> ::= "algebraic expression"
  <right side> ::= "algebraic expression"
  <loop> ::= {DO,{<parameter>,<from>,<to>,<step>},<c-system>}
  <parameter> ::= "identifier"
  <from> ::= <i-expression>


                                        30






  <to> ::= <i-expression>
  <step> ::= <i-expression>
  <i-expression> ::= "algebraic expression" with natural value
                                        (evaluated in FORTRAN)
  <c-system> ::= <part of c-system> | <part of c-system>,<c-
       system>
  <part of c-system> ::= {<variable>,<equation>}

The first and second  argument of the GENLINBANDSOL  statement specifies
the  number  of  the  lower  (below  the  main  diagonal)  and the upper
diagonals  of  the  band  matrix  of  the  system.  The system of linear
 algebraic equations is specified by means of lists expressed by  braces
{ } in the  REDUCE system. The variables  of the equation system  can be
identifiers, but most  probably they are  operators with an  argument or
with arguments that are analogous to array in FORTRAN. The left side  of
each equation has  to be a  linear combination of  the system variables,
the right side, on the contrary, is not allowed to contain any variables
of the system.  The sequence of  the band matrix  lines is given  by the
sequence  of  the  equations,  and  the  sequence  of the columns by the
sequence of the variables in the list describing the equation system.
     The meaning  of the  loop in  the system list is similar to that of
the DO loop  of  the  FORTRAN  language.  The  individual  variables and
equations described by the loop are obtained as follows:

  1. <parameter> = <from>.
  2.  The  <parameter>  value  is  substituted  into the variables and
  equations of the <c-system> loop,  by  which  further  variables and
  equations of the system are obtained.
  3. <parameter> is increased by <step>.
  4. If <parameter> is less or equal <to>, then go to step 2, else all
  variables and equations described  by  the  loop  have  already been
  obtained.

The variables  and equations  of the system included in the loop usually
contain the loop parameter, which mostly occur in the operator arguments
in the REDUCE language, or in the array indices in the FORTRAN language.
     If NL  = <n-lower>, NU = <n-upper>, and for some loop F = <from>, T
= <to>, S = <step> and  N is  the number  of the  equations in  the loop
<c-system>, it has to be true that

  UP(NL/N) + UP(NU/N) < DOWN((T-F)/S)

where UP  represents the  rounding-off to  a higher  natural number, and
DOWN the rounding-off to a lower natural number. With regard to the fact
that, for  example, the last variable before the loop is not required to
equal the last variable  from  the  loop  system,  into  which  the loop
parameter equal  to F-S  is substituted,  when the  band matrix is being
constructed, from the FORTRAN loop that corresponds to the loop from the
specification   of   the   equation   system,  at  least  the  first  NL
variables-equations have to be moved to precede the FORTRAN loop, and at



                                        31






least the  last NU  variables-equations have  to be moved to follow this
loop in order that the correspondence  of the  system variables  in this
loop  with  the  system  variables  before  and  after this loop will be
secured. And this move requires  the  above  mentioned  condition  to be
fulfilled.   As, in  most cases, NL/N and NU/N are small with respect to
(T-F)/S, this condition does not represent any considerable constrain.
     The loop parameters <from>, <to>, and <step> can be natural numbers
or expressions that must have natural  values in the run of the  FORTRAN
program.


6.2 Choosing the numerical library
    ------------------------------

     The user can choose the routines of which numerical library will be
used in the generated FORTRAN  code.  The supported numerical  libraries
are:   LINPACK,  NAG,  IMSL  and  ESSL  (IBM  Engineering and Scientific
Subroutine Library) . The routines DGBFA, DGBSL (band solver) and  DGTSL
(tridiagonal solver)  are used  from the  LINPACK library,  the routines
F01LBF, F04LDF (band solver) and F01LEF, F04LEF (tridiagonal solver) are
used from  the NAG  library, the  routine LEQT1B  is used  from the IMSL
library  and  the  routines  DGBF,  DGBS  (band  solver)  and DGTF, DGTS
(tridiagonal solver)  are used  from the  ESSL library.   By default the
 LINPACK library  routines are  used. The  using of  other libraries  is
controlled by the switches NAG,IMSL and ESSL. All these switches are  by
default OFF. If the switch IMSL  is ON then the IMSL library  routine is
used. If  the switch  IMSL is  OFF and  the switch  NAG is  ON then  NAG
library routines are used. If the switches IMSL and NAG are OFF and  the
switch ESSL is ON then the  ESSL library is used. During generating  the
code using LINPACK, NAG or  ESSL libraries the special routines  are use
for systems with tridiagonal  matrices, because tridiagonal solvers  are
faster than the band matrix solvers.


6.3 Completion of the generated code
    --------------------------------

     The GENLINBANDSOL statement generates a block  of FORTRAN  code ( a
block of  statements of the FORTRAN language) that performs the solution
of the given system of linear algebraic equations. In order  to be used,
this  block  of  code  has  to  be  completed with some declarations and
statements, thus getting  a  certain  envelope  that  enables  it  to be
integrated into the main program.
     In order  to be able to work, the generated block of code has to be
preceded by:

  1. The declaration  of arrays as  described by the  comments generated
  into the FORTRAN code (near the calling of library routines)

  2. The assigning  the values to  the integer variables  describing the
  real  dimensions  of  used  arrays  (again  as  described in generated
  FORTRAN comments)

                                        32






  3. The filling of the variables that can occur in the loop parameters.

  4. The filling or declaration of all variables and arrays occurring in
  the system equations, except for the variables of the system of linear
  equations.

  5. The definition of subroutine ERROUT the call to which is  generated
  after some routines found that the matrix is algorithmically singular

     The  mentioned  envelope  for  the  generated  block can be created
manually, or directly using  the GENTRAN program package  for generating
numeric programs. The  LINBAND module itself  uses the GENTRAN  package,
and the  GENLINBANDSOL statement  can be  applied directly  in the input
files of the GENTRAN package (template processing). The GENTRAN  package
has to be loaded prior to loading of the LINBAND module.
     The  generated  block  of  FORTRAN  code  has to be linked with the
routines from chosen numerical library.







References
----------

[1] R. Liska:  Numerical Code Generation  for Finite Difference  Schemes
     Solving.   In  IMACS  World  Congress  on  Computation  and Applied
     Mathematics. Dublin, July 22-26, 1991, Dublin,(In press).























                                        33