ON THE ECONOMICAL SOLUTION METHOD FOR A SYSTEM OF LINEAR ALGEBRAIC EQUATIONS

The present work proposes a novel optimal and exact method of solving large systems of linear algebraic equations. In the approach under consideration, the solution of a system of algebraic linear equations is found as a point of intersection of hyperplanes, which needs a minimal amount of computer operating storage. Two examples are given. In the first example, the boundary value problem for a three-dimensional stationary heat transfer equation in a parallelepiped in R3 is considered, where boundary value problems of first, second, or third order, or their combinations, are taken into account. The governing differential equations are reduced to algebraic ones with the help of the finite element and boundary element methods for different meshes applied. The obtained results are compared with known analytical solutions. The second example concerns computation of a nonhomogeneous shallow physically and geometrically nonlinear shell subject to transversal uniformly distributed load. The partial differential equations are reduced to a system of nonlinear algebraic equations with the error of O(hx1 + h 2 x2 ). The linearization process is realized through either Newton method or differentiation with respect to a parameter. In consequence, the relations of the boundary condition variations along the shell side and the conditions for the solution matching are reported.


Introduction
It is obvious that a vast number of problems in physics, mechanics, and technology is modelled through linear and nonlinear partial differential equations (PDEs, equations of mathematical physics).In the next step, PDEs are usually reduced to a linear or nonlinear set of algebraic equations applying either finite element or finite difference methods.In what follows, only the method enabling solution of a large number of algebraic equations is considered.To solve the mentioned problem, there are two classical approaches that are direct and iterational methods.Perhaps the most popular are the Gauss reduction for a general system, and the relaxation or relaxation matrix methods usually applied for triple-diagonal or block-diagonal matrices.The mentioned methods belong to the direct ones.Their efficiency depends on the system equation order and on the matrix structure.
From the point of view of iterational methods, a system of equations is treated as the first-order operator equations of the form Ax = b, where A = (a ik ) is a square matrix of dimension n × n, x=(x 1 ,x 2 ,...,x n ) stands for the sought vector, and b = (b 1 ,b 2 ,...,b 4 ) is a given vector (right-hand side).In the case of the mentioned iterational methods, the system of algebraic equations can be treated as the first-order operator equations with the operator acting in an n-dimensional space H n (A : It is worth noticing that the application of the general theory makes it possible to prove convergence of iterations for Seidel's and upper relaxation methods with the minimal constraints used by operator A. Usually, two types of approaches are applied: (i) for known boundaries γ 1 > 0 and γ 1 ≥ γ 2 for a spectrum of operator A lying in a certain energy space H p and (ii) for the case when boundaries γ 1 and γ 2 are unknown.However, the triangle variational method seems to be more effective in application.
For all numerical methods, the reduction of the infinite-dimensional problem to that of the finite dimension plays a key role.It is expressed by the fact that a computational algorithm should yield a solution of the initial problem with a given accuracy ε > 0 through a finite number Q(k) of actions.Moreover, the algorithm should be practically realized, that is, should yield a solution to the problem within the required computer time.Furthermore, the number of actions (and hence the time of solution) Q(ε) should be minimal for the considered problem.Algorithms with the mentioned properties are called economical.It is obvious that the choice of a numerical method of solving a system of linear algebraic equations depends on many circumstances, that is, on the matrix A properties, on the computational type applied, and so forth.
A computational type is one of the possible formulations of the problem, like finding a solution to one special problem of A x = b, or finding solutions to a few variants of problem A x = b with the same matrix A and different right-hand sides B. It may happen that a nonoptimal choice of the problem with one matrix A x = b can be suitable for multivariant computations.
Note that for multivariant computations, one may decrease the average number of operations for one variant if some quantities are conserved and not computed once more for each variant.It depends on a computer and its operating storage.
The choice of a computational algorithm should depend on the computational type, on the volume of the operating storage, and on the considered system order.
In this work, the presented method of solution of large algebraic equation systems offers reduction of the infinite-dimensional problem of mathematical physics to the finitedimensional one through the finite difference approximation.
The used approach yields a band matrix possessing only a few diagonals with nonzero elements.Direct methods of solving those systems matched with the transformation of input matrix A do not enable arrival at similar results, that is, the possibility of solution of the system of high-order algebraic equations offered by them is rather doubtful.

Elimination method for equations
Consider a direct method [6], further referred to as the elimination method for equations, making it possible to find the system solution as a point of intersection of hyperplanes and Jan Awrejcewicz et al. 379 requiring a minimal operating storage volume.The latter property is rarely considered by the standard direct methods [3].
The idea of the proposed method will be illustrated and clarified when applied to a system of linear algebraic equations (SLAE).
Consider first SLAE I of the general form n k=1 a ik x k = b i i = 1,2,...,n. (2.1) The process of the proposed elimination consists of the following few steps.
Step 2. Compute i = j − 1 and construct set M of points that are intersections of straight lines going through points {X j ,X k }, k = 1,2,...,i, and the hyperplane defined by the ith system of algebraic equations of the form n l=1 a il x l = b i , and the coordinates of the kth point of set M i are computed in the following way: The set M i of points {X k }, k = 1,2,...,i, belongs to the hyperplane defined through the ith equation of system (2.1).The coordinates of the points of the mentioned set satisfy the last n − i equations of system (2.1).In Table 2.1, coordinates of set M i for i = n − 1 are reported.
Step 3. Assume j = i and go to Step 2. Repeating computations n − 1 times, one gets a point in an n-dimensional space, which is a solution to system (2.1).
The above considerations allow the conclusion that the above algorithm is suitable for solving SLAE through computations of the matrix rows.
Observe also that for a chosen parametrization of set M n and a given method applied for the construction of the straight lines, only the coordinates It is worth noticing that the proposed novel approach becomes four times smaller, which enables a considerable decrease in the volume of the computer operating storage used to keep coordinates of points X k .
Consider SLAE II with a band matrix: where l is the number of nonzero diagonals and I is a set consisting of l elements, that is, ordered numbers of nonzero band diagonals with an account of sign with respect to the main diagonal.
The following parameters are further introduced: l 1 , m 1 are the numbers of nonzero diagonals and width of the band part lying under the main diagonal, respectively; l 2 , m 2 are the numbers of nonzero diagonals and width of the band part lying over the main diagonal; and m = m 1 + m 2 + 1 is the bandwidth.
Observe that (2.7) The algorithm for solving SLAE (2.6) runs as follows.
Step 1. Assume j = n and solve the jth equation of (2.6) m 1 + 1 times with respect to X n introducing variables x n−m1 ,x n−m1+1 ,...,x n−1 .As a result, the set M j of points X k in Jan Awrejcewicz et al. 381 the (m 1 + 1)th-dimensional space of variables x n−m1 ,x n−m1+1 ,...,x n is obtained: (2.8) Note that it is advisable to take c = 1, p = −1.
Step 2. Compute i = j − 1.If i > m 1 , then the set M j is modified by introducing the variable and the points (2.11) Construct a set of points attained through intersection of the hyperplane defined by the ( j − m)th equation of system (2.6) with the straight lines going through points (2.12) are solved.System (2.12) yields (2.13) Then coordinates of the kth intersection point are computed: Step 3. Assume j = i and return to Step 2. In Table 2.2, coordinates of the points of set M j for m 1 < j n − m 2 are reported.Notice that for the repeated index combination, the following notation is introduced: Remark.The already mentioned advantages with respect to the chosen parametrization of the input points and straight line construction methods hold also for this case.
Step 4. By repeating the computations through Steps 2 and 3 n − 1 times, a point in the space with coordinates x 1 ,x 2 ,...,x m2 is obtained.The space with obtained coordinates overlaps with the values of the unknown x 1 ,x 2 ,...,x m2 occurring in the solution of system (2.6).
Step 5.The order of the initial system is decreased first by neglecting m 2 equations and by putting the found values of the m 2 unknowns into the remaining equations.Notice that during this operation, the parameters of the truncated band matrix do not change.
Step 6. Parameters of the next group are chosen from the m 2 unknowns when the computation through Steps 1, 2, and 3 of the truncated system is carried out.
The described procedure of the successive decrease in the system order is repeated until all values of n unknowns are found.Generally speaking, one may require to get n and m 2 as integers.Notice that by adding a zero diagonal to the compressed band matrix, parameter m 2 may be introduced arbitrarily.
Using the presented method (see [6]), it is possible to apply other variants of SLAE solution to the band matrix as well.
The proposed novel algorithm is stable with respect to the errors introduced by roundings.In this sense, as practical computations show, this method is close to that of the Gauss elimination technique for unknowns with a partial choice of the main element.For instance, solving the system of equations considered in [6, page 61], for the matrix conditioned by ω = 4.7 • 10 5 , the proposed method gives exact decimal digits, that is, of one order less than the Gauss method.This result is obtained owing to a particular symmetry of the fundamental computational formulas (2.13), (2.14) and a homogeneity in computations of all unknowns.
Owing to the algorithm, one may conclude that in order to store coordinates of the points of set M j , one may isolate independently of the order of the sought system, that is, the set consisting of (m 1 + 2)(m 2 + 3) elements.From this point of view, the considered direct method is similar to iterational methods owing to the required volume of operational storage.
Furthermore, the proposed method is characterized by a real computational process cycling and a weak coupling of the algorithm due to constant and full refreshments of transitional results.
Indeed, giving the corresponding coordinates of the points of manifold M for a certain j < n, one may begin the process of solution of the system of equations from the equation with the number i = j − 1.
The mentioned property accounts for an essential increase in the solution process on a computer owing to the elimination of computation repetitions for truncated systems.Namely, while solving an input system of order n, we remember the coordinates of the points of the sets M j , j = n + 1 − lm 2 , l = 1,2,...,n/m 2 − 2, j m 2 + 1Q.
In order to store the values of coordinates, it is advisable to use external computer memory.Next, after a successive truncated system of order n − k m2 , k = 1,2,..., has been formulated, the corresponding set M j of points j = 2m 2 + 1 + (k − 1)m 2 is created.
The discussed algorithm for solving SLAE with a band matrix is also suitable for solving boundary value problems using the method of finite differences in a rectangular parallelepiped since the matrix of SLAE possesses a band with a regular structure.The latter property allows the use of one and only one set I for all equations.
For the space with a solution to a boundary value problem slightly differing from the canonical one, an artificial way of introducing fictitious equations for the mesh nodes to SLAE is as follows: which supplements the given space to a canonical one.
We dwell now for a while on the computation of elements of set I ordered by the numbers of nonzero diagonals of a band of the SLAE matrix with respect to its main diagonal.In what follows, a three-dimensional case is considered.
Let n 1 , n 2 , n 3 be the numbers of mesh nodes in directions ox 1 , ox 2 , ox 3 , respectively.For approximation of partial derivatives of second order in the node {i, j,k}, the seventh point pattern is used: {i, j,k;i ± 1, j,k;i, j ± 1,k;i, j,k ± 1}, where 1 i n 1 , 1 j n 2 , 1 k n 3 .Observe that the number of nonzero band diagonals is equal to that of nodes in the pattern l = 7.The position of nonzero diagonals with respect to the main diagonal of the band matrix is defined using ordered numbers of the unknowns corresponding to the nodes appearing in the pattern.
Numbering of unknowns is defined by the number of mash nodes during SLAE formulation.Although various variants are possible, only those which formulate a band matrix of a regular structure are interesting, that is, one set I gives the position of nonzero coefficients in all equations.
The variant associated with the successive picking of nodes in directions ox 1 , ox 2 , ox 3 and corresponding to the formula for the computation of unknown number n in the node {i, j,k} reads (2.16) Formula (2.16) yields (2.17) For the last band matrix parameters, one gets (2.18) The formula used for the successive picking of nodes along directions ox 3 , ox 1 , ox 2 , that is, the formula applied to compute an unknown number in the node {k, i, j}, is where elements i s of set I read and the band parameter matrix follows: Finally, the last variant of the successive node picking in directions ox 2 , ox 3 , ox 1 corresponds to the formula for computation of the unknown number in the node { j,k,i}, that is, where elements i s of set I read and the band element matrix has the form It follows from the algorithm description that both the parameter m-band SLAE matrix width and the system of order n determine the dimension of the operating storage and the machine time required to solve the equations on a computer.
Therefore, if values n 1 , n 2 , n 3 are different, then the variant of mesh nodes preserving the minimal matrix bandwidth should be chosen.It is clear that the last picking procedure should be associated with the maximal number of nodes in direction n x .
In practice, when the SLAE matrix in the form of a compressed band is built (only nonzero elements are considered), only one picking variant of the shell nodes (independently of the values of n 1 , n 2 , n 3 ) is taken.Therefore, minimization of the bandwidth is introduced, which is reduced to the exchange of elements in rows, exchange of rows, and exchange of column elements of free terms.
In the considered example, the correspondence between elements of the band row is given in Table 2.3 for three variants of picking of mesh nodes.Relations (2.16), (2.19), and (2.22) are used for the exchange of rows and of column elements of free terms, which define simultaneously the ordering numbers of the system equations for the various variants of picking of the shell nodes.
If one considers the boundary value problem for a system of differential equations with respect to a few sought functions, then with the use of finite difference method, the structure of the SLAE matrix becomes dependent on the choice of the numbering of unknowns.
In the case of successive numeration of unknowns into groups corresponding to sought functions, the matrix has a block structure.The bandwidth exceeds the system order, which makes the considered algorithm noneffective.
Note that numeration of similar unknowns with a step equal to the number of sought functions corresponds to successive numeration of all (related to a mesh node) unknowns, and the SLAE matrix possesses a regular band-type structure.The bandwidth is essentially smaller than the system order, and hence the proposed method can be applied.
In practice, both numeration methods mentioned above are applied.The first one is used in the computation of nonzero elements of the equation matrix.The second one is applied during transformation of the matrix to a form suitable for the application of the direct solution method.
It is worth noticing that using the described approach, the general properties of the matrix transformation algorithm are applied independently of the number of the sought functions.The procedure consists of the following steps.
Step 1. Elements of set I, that is, ordering numbers of nonzero diagonals of the matrix with respect to the main diagonals, are computed.For this purpose (analogously to (2.16), (2.17)), ordering numbers of the unknowns corresponding to the nodes of the difference pattern scheme are computed first with respect to each group of the same unknowns in accordance with the first numbering order.
Step 2. Elements of set I k are computed, that is, the relative numbers of these unknowns for the cases when the unknowns lie on the main diagonal associated with the first, the second,..., the kth group of the unknowns named the same.
Step 3. A nonorder set I matching all sets I k is constructed.Further, an increasing ordering of the elements of set I is carried out and the number of elements of set I with parameters Step 4. The reconstruction of the initial block-band compressed matrix is initiated through changes of row elements, rows, and free column elements with respect to the second numbering way for unknowns: first elements of all blocks, second elements of all blocks, and so on; first rows of blocks corresponding to the equations of the boundary value problem, second rows of these blocks, and so forth.
Elements of the column of free terms are transformed in a way similar to that applied to the rows.
Step 5. Then a band-type matrix in the compressed form is built in the following way: ordering numbers of the row elements from the block of the input matrix are defined through the condition of equality of the corresponding elements of sets I k , I.
Notice that in its essential part, the resulting band is filled with zero elements.Consider an example of the SLAE matrix transformation occurring during the solving of the boundary value problem for the system of two differential equations in partial derivatives of fourth order with variable coefficients in two-dimensional space using the method of finite differences.Denote by n 1 , n 2 the numbers of mesh nodes in directions ox 1 , ox 2 , respectively.
The derivatives in mesh node {i, j} are approximated through difference relations on the pattern consisting of 13 nodes (see Figure 2.1).
The following coordinates of the local nodes {i s , j s }, s = 1,2,...,13, are introduced: i 7 = i, j 7 = j; i 8 = i, j 8 = j + 1; i 9 = i, j 9 = j + 2; (2.25) The input matrix in the compressed form (nonzero diagonals) is built owing to the first numbering way for unknowns, and it describes a two-dimensional packet of dimension k n • k l , where k = 2 is the number of sought functions, n = n 1 n 2 is the number of mesh nodes for which equations are constructed, and l = 13 is the pattern dimension.
The discussed matrix is a block-diagonal one consisting of k 2 blocks and possessing 39 nonzero diagonals.
Local index s of the pattern node defines a position in the row of a coefficient standing by the corresponding unknown in the equation for node {i, j}.The computation is carried out with respect to the beginning of the corresponding block.
Then, elements of the following sets are computed: For this purpose, the relative members of all unknowns corresponding to pattern nodes and appearing in equations for node {i, j} are computed: where even k : The set I = I 1 ∪ I 2 is constructed and its elements are ordered in an increasing manner.The corresponding results for n 2 = 5 are reported in the first four columns of Table 2.4.
As seen from this table, the number of nonzero diagonals of the band matrix is equal to 31.This value does not depend on n 2 but is defined only through a type of the pattern chosen for approximation of derivatives and by the number of the functions sought in the boundary value problem.
Comparing elements of sets I 1 , I 2 , and I, one obtains a rule of distribution of the row elements in the modified matrix during construction of the band matrix in the compact form.
The correspondence between the ordering numbers of nonzero coefficients for the even and odd rows of the modified matrix and rows of the band matrix is reported in the On the economical solution of algebraic equations  As known, when solving numerically boundary value problems for PDEs using finite element methods, there appears SLAE associated with matrices including a relatively small number of nonzero elements, most of which remain on the main diagonal, that is, in this case, the equations dealt with are associated with the band matrix of irregular structure (the so-called cutting-off matrix).
However, there exists a modification of the elimination method for equations due to which it is possible to solve SLAE of a similar form.The cutting-off matrix will be introduced through two sets: the set of nonzero matrix elements and the corresponding set of indices, that is, relative numbers of nonzero elements with respect to the elements lying on the main diagonal, which are associated with the indices.
In what follows, contrary to the band matrix of a regular structure requiring only one row of indices for nonzero elements (set I), in our case, this type of row is formulated for each of the equations (set I i ).Here lies the fundamental difference in the application of the elimination method for the equations of SLAE using the cutting-off matrix.
We comment on parameters m 1 , m 2 .Being the same for all equations, they are computed initially as where l 1 denotes the number of nonzero coefficients in the ith equation.Therefore, the row of nonzero coefficients and the row of indices for each equation are supplemented, if necessary, by zeros in the case of elements, and by in the case of indices, where l i is the value of I i satisfying the above requirement.As a result, SLAE with a band matrix of nonregular structure is obtained.Note that the formulation of the cutting-off matrix with the help of two packets is optimal with respect to the memory storage volume since for the index conservation, it is sufficient that the machine word be 2 bytes long.
The corresponding computation shows, for instance, that for the conservation of a band-type matrix of the 27th order and a regular structure obtained through approximation of the Laplace operator in a three-dimensional space by finite difference equations, 5832 bytes are required in the general case, 1512 bytes are required in the case of its compressed form (only nonzero elements are considered), and finally 1350 bytes are needed when the matrix is represented by two packets.
It is worth noticing that the method of elimination extends essentially the possibility of computer computation needed to solve large-order SLAE through reaching a match between the advantages of both direct (universality, finiteness of the computational process, exactness) and iterational (minimal requirement of the storage volume) methods.
We finally emphasize particularities of the SLAE matrix in each separate case, which can also be treated as an advantage of the method.

Numerical solution of a three-dimensional equation of elliptic type
Many stationary processes of different physical properties lead to PDEs of elliptic type, to mention only the stationary problem of current distribution in a medium, problems of electrostatics and magnetostatics of the theory of plates and shells, problems of the theory of elasticity or the theory of filtration, and so forth.
Exact solution to the boundary value problems for elliptic equations can be found in rare cases only.Therefore, numerical methods are usually applied to solve differential equations of elliptic type (linear or nonlinear).In the latter case, equations are linearized through differentiation along the parameter [7] or using the "setup" method [2].Hence, the problem is finally reduced to the solution of the SLAE.
Moreover, a majority of methods for nonlinear problem solving are reduced to a series of linear systems.Application of the method of finite differences to elliptic equations results in SLAE with a band-type matrix.In the case of PDEs, the matrix (SLAE) possesses the cutting-off band, with only a few nonzero elements (see Section 2).
We briefly model the problem of stationary heat distribution in a certain volume G with surface Γ of three-dimensional space x = (x 1 ,x 2 ,x 3 ).The heat transfer equation is governed by the Fourier principle.A vector of heat stream density W is proportional to the temperature gradient V = V (x) such that where K = K(x) is the heat transfer coefficient.The density of the heat stream is equal to the amount of the heat stream, a unit of time passing through a unit area of an isotermic surface.
To derive the equation governing the heat balance for certain volume U ∈ G and surface S, let the distributed heat sources exist inside volume U with density ϕ(x), where ϕ(x)dU is the heat amount occurring in volume dU.
Let W n describe the vector W projection onto external normal n to surface S. The heat balance equation is governed by the known rule, that is, the total heat stream passing through surface S, S W n dS, (3.2) should be equal to the heat amount of Using the Gauss formula S W n dS = U div, balance equation (3.4) is transformed to the form As volume U is arbitrary, therefore, if ϕ(x) and divW are continuous functions of point x = (x 1 ,x 2 ,x 3 ), then (3.5) yields Jan Awrejcewicz et al. 391 Substituting here expression (3.1) for the vector of heat stream W, the following equation for the stationary temperature V = V (x) is obtained: or in the equivalent form, where K is the function of point x = (x 1 ,x 2 ,x 3 ).Heat transfer equation (3.8) is obtained under an assumption of isotropy of the heat transfer process.If a heat coefficient depends on direction and is a tensor (an isotropic medium), then, instead of (3.8), the following equation is taken: (3.9) (3.10) Equation (3.8) holds for all internal points of space G.Additional conditions of the following form are attached to boundary Γ: (1) temperature is given: V (x) = g(x) for x ∈ Γ; (2) heat stream is given: K(∂V /∂n) = g(x) for x ∈ Γ; (3) heat transfer of Newton's rule is applied: Consider the boundary value problem for the equation of the stationary heat transfer in a rectangular parallelepiped of a three-dimensional space: (3.12) In this case, the method of finite differences is used; the construction of the difference equation on the mesh with a variable step is realized through the method of functional approximation [5].Consider the boundary value conditions (1), (2), and (3) in various states, and, if necessary, symmetry q of solution is accounted: where K α (x) is the heat transfer coefficient in direction x α , x ±α (x) is the heat exchange coefficient with surrounding medium on walls x α = 0,l α , and V α (x) is the temperature of surrounding medium.
It is known that a solution to the boundary value problem (3.13), (3.14) satisfies the functional minimum: where g ±α (x) = ae ±α (x)V 0 (x).Indeed, variating functional (3.15) with respect to V , V ,xα , integrating by parts the terms with δV ,xα , and comparing variations of the functional to zero, relations (3.4) and (3.14) are obtained as the necessary conditions to minimize functional (3.15).
Further, instead of V and K, v and k will be used.We introduce the following fundamental mesh: where ωα = {x αi , i = 0,1,2,...,N α − 1,N α }-mesh on interval [0, l α ].The set of interval mesh nodes is denoted by where ω α = {x αi , i = 1,2,...,N α − 1}.We introduce the following notations: In what follows, the function with n variables v i, j,k is obtained: As known, the minimum of this function is achieved at the point where its partial derivatives are equal to zero with respect to v i, j,k .We construct a system of difference equations of order n to define values v i, j,k realizing the minimum of the functional (3.15) and being a solution to the boundary value problem (3.13), (3.14).The following four cases are considered.
(1) Internal mesh nodes of space Ḡ.Terms (3.19) including v i, j,k are reported below: (3.20) By differentiating S 1 with respect to v (i, j,k) and comparing its derivative to zero, the following difference equation is obtained at node {i, j,k} of the mesh owing to the division by 2h 1i h 2 j h 3k : (2) Node of the parallelepiped wall.
By differentiating S 2 with respect to v i, j,k and comparing the corresponding derivative to zero, the following difference equation (after division by 2 h1i h2j h3k ) is obtained at the parallelepiped node x 1 = 0: . (3.24) Difference equations at other parallelepiped wall nodes are constructed analogously.
Jan Awrejcewicz et al. 395 (3) Mesh nodes belong to the parallelepiped By proceeding in a way similar to that applied in the previous case, the following difference equation is obtained: .
(3.26) (4) A node of the mash coincides with the parallelepiped vertex Proceeding again as in the previous case, one arrives at the following corresponding equations: (3.28) (3.29) Thus, all possible cases occurring while formulating the difference equations on the mesh with variable step for the third boundary value problem have been discussed.We discuss briefly the peculiarities involved in the formulation of difference equations in the case when the following Dirichlet conditions are given on some parallelepiped walls: where v u (x) is the given function.
One of the possible ways to include a condition of such a type is to formulate the following corresponding equation associated with nodes: (3.31) The latter one allows for construction of both an algorithm and computer programs universal with respect to the boundary condition type.Indeed, the drawback of the approach is associated with conservation of the order of the difference equation system.Furthermore, if it is known a priori that the solution to problem (3.13), (3.14) possesses one, two, or three plane symmetries x α = l α /2, then these conditions are recommended to be included while difference equations are being formulated, since they essentially decrease the order of the investigated system.
If the solution is symmetric with respect to, say, the plane x 1 = 1 /2, then we take In what follows, in the difference equations formulated for the nodes of the plane x α = α /2, accounting for the symmetry property the coefficient standing by the unknown v i−1, j,k is doubled.The wall node is treated as an internal one in this case.The solution to the boundary value problem (3.13), (3.14) consists of two steps: formulation of the system of difference equations with a band-type compressed matrix having a regular structure (in the operating storage, only nonlinear matrix diagonals are conserved), and solution to the obtained system using the elimination method for equations [3].
We discuss the computations of the first stage.In order to identify a particular problem, the following items are required: node numbers N α with respect to each direction x α , α = 1,2,3; packets of steps of the mesh h 1i , h 2 j , h 3k ; formulas for k α (x), g ±α (x), ϕ(x) computations, and v u (x) (if necessary); values of symmetry indicators IS, JS, KS; packet of six elements to formulate the type of the boundary condition on each of the six parallelepiped walls; and values of χ ±α (x).
Having this information, the operator coefficients k i, j,k 1 , k i−0,5, j,k 1 are computed and the matrix of difference equations is formulated due to formulas (3.21)-(3.29).Analogous operations are carried out with k 2 (x), k 3 (x).
Jan Awrejcewicz et al. 397 In what follows, the matrix of coefficients of the difference equations, which corresponds to the third boundary value problem, is obtained.Then free terms of the system of equations are computed, and (if necessary) the equations involving the Dirichlet conditions are formulated.
Recall that when solving numerically the boundary value problems, the reality estimation of the computer-yielded results is in the focus of attention.Therefore, the correct derivation of appropriate relations together with suitable formulation of the variational problem should be emphasized, and then an application of the exact method of solution to the system of difference equations with double precision should be carried out.
During analysis, a series of computations to verify the mentioned reliability has been carried out.
(1) Model problem.The possibilities of both the developed algorithm and the system subroutines are checked through solution to problem (3.13), (3.14), with , where The exact solution reads and hence We introduce the Dirichlet condition on the wall x 1 = 0.The symmetry of solution is taken into account with respect to the plane x 1 = 1/2, x 2 = 1/2.The solution is defined in the space Ḡ = {0 x α 1, α = 1,2,3}.In Table 3.1, the results of the solution to this problem for c = 2048 are reported for the mesh with a constant step.An analysis of the given results yields the conclusion that during the change of the mesh step, the approximated solution converges to the exact one.Variation of the operator coefficients in the considered case is high, that is, Notice that for c = 0, the numerical solution with accuracy greater than five meaningful digits overlaps with the exact solution.
(2) Problem without a source.The convergence of the solution to the problem without a source (owing to the decreasing mesh step) is testified through the results included in Table 3.2.The following parameters are taken: and the following boundary conditions are applied: The following symmetry solution with respect to the planes is taken into account: Owing to the analysis of the results reported in Table 3.3, it is clear that without sources, the numerical solution is obtained with an efficient accuracy on the mesh 5 × 5.
(3) Problem with a point source.The convergence of the solution to the problem with a point source obtained during the decrease of the mesh step is considered on the basis of the following problem: where the attached boundary conditions are as follows: The following symmetry condition with respect to the planes is applied.In Tables 3.4 and 3.5, the results of the solution to this problem on four meshes, 3 × 3 × 5, 5 × 5 × 5, 9 × 9 × 5, and 17 × 17 × 5, are given.For the three partitions, the mesh step (in the vicinity of the node with the source) is taken to be constant and equal to the step of partition 17 × 17 × 5.In all four cases, the constant volume power is achieved.The changeable step allows application of the mentioned actions in a relatively simple manner by introduction of additional nodes with respect to x 1 , x 2 in the first three cases, that is, the solution is obtained for the shells 4 × 4 × 5, 6 × 6 × 5, and 10 × 10 × 5.
Owing to the given results, stable convergence is achieved during the mesh step variation.Maximal differences are observed at the node with the source, but the difference is less than 1% for two neighboring partitions.Finally, we present the results of computation through the finite difference method for two problems for which both analytical and boundary value solutions are known.
In Table 3.6, the results of the solution to the first boundary value problem for the Laplace equation for the cube with ribs of unit length are reported.
The boundary conditions are as follows: (3.45) In Table 3.7, the results of the solution to the Laplace equation with a hybrid boundary condition for the rectangular parallelepiped are reported.The boundary conditions are Notice that when numerical solutions are constructed, the symmetry with respect to two planes, x 2 = 0, x 3 = 0, is taken into account.
The reported results show that the solution obtained through the finite difference method during the mesh step decrease converges to the analytical solution, and its error is small for the appropriate partition.The results of the numerical experiments yield the conclusion that our proposed method and the associated subroutines allow efficient numerical solution of three-dimensional boundary problems for the stationary heat transfer equation without any essential limitations.The same conclusion holds when the problems of plates and shells are considered.

Computation of geometrically nonlinear nonhomogenous shallow shells with mixed boundary condition along their sides
One of the important computation problems connected with shallow nonhomogeneous and geometrically nonlinear shells concerns conditions of the boundary value problems variations along a supporting side.
In order to derive appropriate relations, a system consisting of two identical shells stiffened by ribs along two sides linked by a hinge is considered (see Figure 4.1).
The shell material is assumed to be isotropic but nonhomogeneous, that is, shear modulus G and Poisson coefficient µ are functions of point coordinates x = (x 1 ,x 2 ,x 3 ), or, according to the theory of small elastic-plastic deformations, they depend on the stressstrain material state at the considered point.
Following the Kirchhoff-Love hypothesis, the components of shell and rib deformation read where , ε 11 , ε 22 , ε 12 -length and stream deformations of the mean surface; where χ 11 , χ 22 , χ 12 are variations of curvature and torsion.
In accordance with the Hook principle for a flat-strain state and (4.1), the coupling between stressed σ ik and deformation e ik is presented in the following form: where, here and further, (x ←→ y) denotes the change of index x to y. where The deformation energy of the considered system reads where 1 -rib stiffness on elongation (compressing), bending in a vertical plane, bending in a horizontal plane, and torsion, respectively; α ± 1 , α-stiffness of hinge ribshell, rib-rib interaction, respectively, β ± 1 , β-elasticity of support under the rib; [w ,x1 ]jump of derivative w ,x1 on the line x 1 = 0.
We introduce the following functional: Conditions for the coupling of solutions on the line x 1 = 0 are derived from the following expression: where This means that the following conditions should be satisfied on the line of the solution coupling: as well as the conditions on the continuity of displacements ), the conditions for the solution coupling for a plate are obtained [5].
We determine agreement conditions in the shell angle x 1 = 0, x 2 = 0.For this purpose, the forms related to an angle are considered in (4.15), that is,

.20)
One may derive now, for instance, the following compatibility conditions in the angle: 406 On the economical solution of algebraic equations Notice that if the ribs are excluded, then only one compatibility condition M 12 = 0 remains.
We derive compatibility conditions for the point x 1 = 0, where the change of boundary conditions occurs.They are obtained from the terms of (4.15) with respect to the point x 1 = 0, x 2 = 0, and the coefficients with index 2 are equal to zero (homogeneous shell):

.22)
As an example, we consider the compatibility condition for a hinged support of the shell side x 2 = 0 with different stiffness: Owing to relations (4.22), (4.23), the necessity of application of the compatibility condition appears only for the shells supported by ribs.
We introduce the following force function: and we derive the equations in a hybrid form with respect to force F and deflection w.
For this purpose, the deformation continuity equation and the third equation of (4.15) are applied to yield By solving (4.6) with respect to ε ik , one gets where where In what follows, the algorithm for solving the stability problem of physically nonlinear flexible shallow shells with mixed boundary conditions and without ribs is discussed.
Step 1.In the space occupied by the shell, a rectangular mash with equal step {x 1i ,x 2 j ,x 3k } is introduced.Assuming in the beginning that the shell is unloaded, modulus G i jk of the shell nodes is introduced.
Step 2. At the nodes of mesh surface x 3 = 0, coefficients C α β (4.11) are computed.Then, substituting the partial derivatives by difference expressions with error Q(h 2 x1 + h 2 x2 ), instead of the system of equations (4.30), the system of nonlinear algebraic equations with respect to w i j , F i j , that is, deflection and stress functions formulated at nodes of the mean surface, is obtained.Boundary conditions are considered while out-contour nodes are being eliminated from the equations system.If the mesh node overlaps with the point where the boundary conditions change, for instance, in the case of hinge clamping, then the clamping is understood as the value of derivative with respect to a normal solution to a system of algebraic equations for a given load q(x 1 x 2 ) defined with the help of the Newton or differentiation along the parameter methods [7], and then the elimination method is applied.
Step 3. New values of G i jk and µ i jk are computed at the mesh nodes in accordance with the theory of small elastic-plastic deformation [1,4].Modulus of elongation and modulus of volume deformation K are found using formulas of the theory of small elastic-plastic deformation

.31)
Notice that in this theory, it is assumed that K does not depend on the deformable state at a point since in the case of a body which is nonhomogeneous before deformation, K = K(x 1 ,x 2 ,x 3 ), and for a body homogeneous before deformation, K = K 0 = const.The shear modulus is defined through the formula σ i e i e i (x) , (4.32) and G is called the cutting modulus.In order to compute the shell, dependence σ i (e i ) should be explicitly given (intensity of stress σ i versus strain (e i )).Some of the diagrams σ i (e i ) are reported in [4].Then the procedure goes back to Step 2.
The computations are repeated until the obtained solution overlaps with the previous one keeping the assumed accuracy.
Step 4. The load q = q + h q is increased and the computations start at Step 2.
As a result, dependence q(w) yielding the critical values is obtained.On a basis of the introduced algorithm, the program in Fortran has been developed.
As an example, a square shell made of aluminium and with parameters k 1 = k 2 = 18, λ = a/h = 50, and q = const [4] has been considered.
We apply the Mises flow condition (σ i = σ s ) and the dependence σ i (e i ) taken in the form where e s = σ s /3σ 0 is the artificial intensity of the flow deformation.
The problem is considered symmetric with respect to x 1 , x 2 for four types of the boundary conditions, namely, (1) hinged support on the contour; (2) hinged support in an angle, and clamping in the middle of the side over the interval 1/4a; (3) clamping in an angle, and hinged support in the middle of the side over the interval 1/4a; (4) clamping along a contour.The stress function on the side x 2 = 0 is F ,x1x1 = 0, F ,x1x2 = 0, whereas for the deflection, w = 0, w x2 = 0 or M 22 = 0.
Consider a few computational results.In Figure 4.2, dependencies of load versus center deflection represented by curves 1-4 corresponding to the considered boundary conditions are reported.Analyzing the behavior of the curves, one may conclude that the value of the upper critical load depends essentially on the position of the clamped part.

Figure 4 .
3 shows influence of the boundary conditions on the distribution of plasticity zones on the upper shell surface x 3 = −h/2 for the loads corresponding to deflection w = 0,4 in the shell center.The largely developed plasticity zones are achieved for hinged support.Owing to clamping of the shell point in the middle of its side, the plastic zones are sharply decreased and their position is subject to change.The minimal areas of the plasticity zones are achieved in the case of clamping along the whole length.JanAwrejcewicz et al. 409

Figure 4 . 3 .
Figure 4.3.Influence of boundary conditions on plasticity zone distribution (see text for more details).

Table 2 .
1. Coordinates of points in set M i

Table 2 .
2. Coordinates of points M j (for m 1

Table 2
.4.The free positions in the row of the matrix are filled out by zeros.Notice that for arbitrary n 1 , n 2 , this correspondence is the same.
.18) JanAwrejcewicz et al. 393All integrals occurring in functional (3.15) are substituted by quadratures.If the underintegral function includes a partial derivative, then the center triangle formula is applied with respect to the corresponding variable (in other cases, the trapezoid rule is used).Derivatives are approximated through difference approximations, for instance, by v

Table 3 .
2. Numerical solution of the model problems at point:x 1k = x 2k = x 3k .

Table 3 .
3. Numerical solution to the problem without source.

Table 3 .
4. Numerical results of the problem with point source.

Table 3 .
5. Numerical results of the problem with point source.

Table 3 .
6. Solution to the first boundary value problem for Laplace equation.