@@ -1,1983 +1,1983 @@ - - - - - - - - - - - - - - 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 ,,...,,,...,; - ::= 2 | 3 coordinate system dimension - ::= "algebraic expression" the expression of the i-th - Cartesian coordinate in new - coordinates - ::= "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: - - ,,...,; - ::= VECTORS | DYADS | TENSOR - ::= "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(,) expresses the derivative of -with respect to the coordinate . 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 is described as follows (the -conditions (d=i) represent the limitation on the dimension of the -coordinate system equalling i): - - ::= | | - ::= "algebraic expression, can contain " | - "tensor variable with scalar value" | - . | # | - (d=2)? | DIV | - LAPL | (d=2) ROT | - DIRDF(,) - ::= "identifier declared by VECTORS statement" | - "tensor variable with vector value" | - VECT(,...,) | - | - + | - | - * | / | - . | . | (d=3) - ? | (d=2) ? | - (d=2) ? | GRAD | - - 6 - - - - - - - DIV | LAPL | (d=3) ROT | - DIRDF(,) | DF(,"usual - further arguments") - ::= "identifier declared by DYADS statement" | - "tensor variable with dyadic value" | - DYAD((,...,),...,(, - ...,)) | - | + | - - | * | / - | . | & | - (d=3) ? | (d=3) ? | - GRAD | DF(,"usual further arguments") - - -1.6 Assigning statement - ------------------- - - The assigning statement for tensor variables has a usual syntax, -namely: - - := - ::= "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 {,} [ INTO - {,}]; - ::= "identifier" - the name of the coordinate - ::= "identifier" - the name of the index - -This statement specifies that the will correspond to the -. 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,{,}; - -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 {,}; - ::= ([],{, - }) - ::= "identifier" - the name of the function - ::= 1|2 tensor order of the function (the value of - the function is 1 - vector, 2 - dyad (two- - dimensional tensor)) - -Every in the statement determines on which -the depends. If the tensor of the function occurs in -the , the is declared as a vector or a dyad. If, -however, the has been declared by the VECTORS and DYADS -statements of the EXPRES module, the user need not present the tensor -. 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 {,}; - - 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 {,}; - -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 {,}; - ::= ([,]{,}) - ::= .. , - - 11 - - - - - - - ::= ONE | HALF designation of the integer - (ONE) and half-integer (HALF) - grids - ::= | for the vector - , for the dyadic - it is not presented for the - scalar - ::= *| "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 {,}; - ::= ({,}) - -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 ,,{{,} - , }; - ::= ALL | "identifier" - the coordinate name from - the COORDINATES statement - ::= | - | - | * - | / | - ** | - DIFF(,[,])| - ({,}) - ::= X - ::= U | V | W - ::= F | G - ::= N | "integer greater than 1" - ::= "integer greater than 2" - ::= = - ::= | - - - 13 - - - - - - - ::= "non-negative integer" - ::= ()| - "natural number"|DI|DIM1|DIP1|DIM2|DIP2| - | - | - + | - * | - / | - () | - ** - ::= X | U | V | W | F | G - ::= | - + | - - - ::= I - = "rational number" - DIFCONST {,}; - ::= "identifier" - the constant parameter of - the difference scheme. - DIFFUNC {,}; - ::= "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 - - =1/(X(I+1/2)-X(I-1/2))*DINT(X(I-1/2),X(I+1/2), - ,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 {,,}; - ::= "identifier" - the name of the array for storing - the result - ::= "identifier" - the name of the function - whose behavior is described by the - equation - ::= = - ::= "algebraic expression" , the derivatives are - designated by the DIFF operator - ::= "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
{,
}; -
::= = - ::= "rational number" - MAXORDER {,}; - ::= = - ::= "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 + , 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 {,}; - ::= "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 (); - ::= = - ::= "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 {,} - ::= "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 ; - ::= "identifier" - -After this statement the new identifiers created by the operators -DENOTEPOL and DENOTEMAT will begin with the letters of the identifier - 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 R0j and the imaginary part by I0j, where -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 Rjk and the imaginary - - 26 - - - - - - -part by Ijk. The returned value of the operator is the original -matrix with non-numeric elements replaced by Rjk + I*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 (,,{}); - ::= "natural number" - ::= "natural number" - ::= | , - ::= {,} | - ::= "kernel" - ::= = - ::= "algebraic expression" - ::= "algebraic expression" - ::= {DO,{,,,},} - ::= "identifier" - ::= - - - 30 - - - - - - - ::= - ::= - ::= "algebraic expression" with natural value - (evaluated in FORTRAN) - ::= | , - ::= {,} - -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. = . - 2. The value is substituted into the variables and - equations of the loop, by which further variables and - equations of the system are obtained. - 3. is increased by . - 4. If is less or equal , 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 = , NU = , and for some loop F = , T -= , S = and N is the number of the equations in the loop -, 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 , , and 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 - + + + + + + + + + + + + + + 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 ,,...,,,...,; + ::= 2 | 3 coordinate system dimension + ::= "algebraic expression" the expression of the i-th + Cartesian coordinate in new + coordinates + ::= "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: + + ,,...,; + ::= VECTORS | DYADS | TENSOR + ::= "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(,) expresses the derivative of +with respect to the coordinate . 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 is described as follows (the +conditions (d=i) represent the limitation on the dimension of the +coordinate system equalling i): + + ::= | | + ::= "algebraic expression, can contain " | + "tensor variable with scalar value" | + . | # | + (d=2)? | DIV | + LAPL | (d=2) ROT | + DIRDF(,) + ::= "identifier declared by VECTORS statement" | + "tensor variable with vector value" | + VECT(,...,) | - | + + | - | + * | / | + . | . | (d=3) + ? | (d=2) ? | + (d=2) ? | GRAD | + + 6 + + + + + + + DIV | LAPL | (d=3) ROT | + DIRDF(,) | DF(,"usual + further arguments") + ::= "identifier declared by DYADS statement" | + "tensor variable with dyadic value" | + DYAD((,...,),...,(, + ...,)) | - | + | + - | * | / + | . | & | + (d=3) ? | (d=3) ? | + GRAD | DF(,"usual further arguments") + + +1.6 Assigning statement + ------------------- + + The assigning statement for tensor variables has a usual syntax, +namely: + + := + ::= "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 {,} [ INTO + {,}]; + ::= "identifier" - the name of the coordinate + ::= "identifier" - the name of the index + +This statement specifies that the will correspond to the +. 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,{,}; + +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 {,}; + ::= ([],{, + }) + ::= "identifier" - the name of the function + ::= 1|2 tensor order of the function (the value of + the function is 1 - vector, 2 - dyad (two- + dimensional tensor)) + +Every in the statement determines on which +the depends. If the tensor of the function occurs in +the , the is declared as a vector or a dyad. If, +however, the has been declared by the VECTORS and DYADS +statements of the EXPRES module, the user need not present the tensor +. 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 {,}; + + 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 {,}; + +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 {,}; + ::= ([,]{,}) + ::= .. , + + 11 + + + + + + + ::= ONE | HALF designation of the integer + (ONE) and half-integer (HALF) + grids + ::= | for the vector + , for the dyadic + it is not presented for the + scalar + ::= *| "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 {,}; + ::= ({,}) + +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 ,,{{,} + , }; + ::= ALL | "identifier" - the coordinate name from + the COORDINATES statement + ::= | + | + | * + | / | + ** | + DIFF(,[,])| + ({,}) + ::= X + ::= U | V | W + ::= F | G + ::= N | "integer greater than 1" + ::= "integer greater than 2" + ::= = + ::= | + + + 13 + + + + + + + ::= "non-negative integer" + ::= ()| + "natural number"|DI|DIM1|DIP1|DIM2|DIP2| + | - | + + | + * | + / | + () | + ** + ::= X | U | V | W | F | G + ::= | + + | + - + ::= I + = "rational number" + DIFCONST {,}; + ::= "identifier" - the constant parameter of + the difference scheme. + DIFFUNC {,}; + ::= "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 + + =1/(X(I+1/2)-X(I-1/2))*DINT(X(I-1/2),X(I+1/2), + ,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 {,,}; + ::= "identifier" - the name of the array for storing + the result + ::= "identifier" - the name of the function + whose behavior is described by the + equation + ::= = + ::= "algebraic expression" , the derivatives are + designated by the DIFF operator + ::= "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
{,
}; +
::= = + ::= "rational number" + MAXORDER {,}; + ::= = + ::= "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 + , 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 {,}; + ::= "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 (); + ::= = + ::= "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 {,} + ::= "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 ; + ::= "identifier" + +After this statement the new identifiers created by the operators +DENOTEPOL and DENOTEMAT will begin with the letters of the identifier + 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 R0j and the imaginary part by I0j, where +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 Rjk and the imaginary + + 26 + + + + + + +part by Ijk. The returned value of the operator is the original +matrix with non-numeric elements replaced by Rjk + I*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 (,,{}); + ::= "natural number" + ::= "natural number" + ::= | , + ::= {,} | + ::= "kernel" + ::= = + ::= "algebraic expression" + ::= "algebraic expression" + ::= {DO,{,,,},} + ::= "identifier" + ::= + + + 30 + + + + + + + ::= + ::= + ::= "algebraic expression" with natural value + (evaluated in FORTRAN) + ::= | , + ::= {,} + +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. = . + 2. The value is substituted into the variables and + equations of the loop, by which further variables and + equations of the system are obtained. + 3. is increased by . + 4. If is less or equal , 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 = , NU = , and for some loop F = , T += , S = and N is the number of the equations in the loop +, 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 , , and 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 +