Artifact ab87ff8129108f11b1446afaac102b5229140e1732c5269f17d5d6aba43045d2:
- Executable file
r36/doc/FIDE.DOC
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 78156) [annotate] [blame] [check-ins using] [more...]
- Executable file
r37/packages/fide/fide.txt
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 78156) [annotate] [blame] [check-ins using]
- Executable file
r38/packages/fide/fide.txt
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 78156) [annotate] [blame] [check-ins using]
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