# On the economical solution method for a system of linear algebraic equations

On the economical solution of algebraic equations From the point of view of iterational methods, a system of equations is treated as the ï¬rst-order operator equations of the form Ax = b, where A = (aik ) is a square matrix of dimension n Ã n, x=(x1 ,...,xn ) stands for the sought vector, and b = (b1 ,b2 ,...,b4 ) 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 ï¬rst-order operator equations with the operator acting in an n-dimensional space Hn (A : Hn â Hn ), x,b â Hn . 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 eï¬ective in application. For all numerical methods, the reduction of the inï¬nite-dimensional problem to that of the ï¬nite 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 ï¬nite 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 ï¬nding a solution to one special problem of Ax = b, or ï¬nding solutions to a few variants of problem Ax = b with the same matrix A and diï¬erent right-hand sides B. It may happen that a nonoptimal choice of the problem with one matrix Ax = 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 inï¬nite-dimensional problem of mathematical physics to the ï¬nitedimensional one through the ï¬nite diï¬erence 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 oï¬ered by them is rather doubtful. 2. Elimination method for equations Consider a direct method [6], further referred to as the elimination method for equations, making it possible to ï¬nd 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 clariï¬ed when applied to a system of linear algebraic equations (SLAE). Consider ï¬rst SLAE I of the general form aik xk = bi k=1 i = 1,2,...,n. (2.1) The process of the proposed elimination consists of the following few steps. Step 1. Assume j = n. Solve the jth equation of system (2.1) j times with respect to xn by keeping variables x1 ,... ,xni ï¬xed. As a result, the set M j consisting of points {X k }, k = k = 1,2,..., j, with the coordinates 1 X 1 c,c,... ,c,xn , 2 X 2 p,c,... ,c,xn , 3 X 3 c, p,... ,c,xn , (2.2) X j c,c,... ,c, p,xn , is obtained, where c, p are arbitrary numbers (it is suitable to take c = 1, p = 1). 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 deï¬ned by the ith system of algebraic equations of the form n j ail xl = bi , l=1 j xl = xl (2.3) xlk j xl Î»k , l = 1,2,...,n, is solved. System (2.3) yields Î»k = bi j n l=1 ail xl j n k l=1 ail xl xl (2.4) and the coordinates of the kth point of set Mi are computed in the following way: Ë xlk = xl + xlk xl Î»k , l = 1,2,...,n. (2.5) The set Mi of points {X k }, k = 1,2,...,i, belongs to the hyperplane deï¬ned 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 Mi for i = n 1 are reported. On the economical solution of algebraic equations Table 2.1. Coordinates of points in set Mi (for i = n 1). x1 x2 xn3 xn2 xn1 xn X1 c c c c 1 xn1 1 xn X2 2 x1 c c c 2 xn1 2 xn X3 c 3 x2 c c 3 xn3 3 xn Â·Â·Â· Â·Â·Â· Â·Â·Â· Â·Â·Â· Â·Â·Â· Â·Â·Â· Â·Â·Â· X i1 c c i xn13 c i xn11 i xn1 Xi c c c i xn2 i xn1 i xn 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 Mn and a given method applied k k k for the construction of the straight lines, only the coordinates xk1 , xik , xi+1 ,...,xn should be determined; other coordinates preserve values equal to c. 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: a +i ,x j+i = b j , s=1 j = 1,2,...,n, 1 j + is n, is â I, (2.6) 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: l1 , m1 are the numbers of nonzero diagonals and width of the band part lying under the main diagonal, respectively; l2 , m2 are the numbers of nonzero diagonals and width of the band part lying over the main diagonal; and m = m1 + m2 + 1 is the bandwidth. Observe that m1 = I(1), m2 = I(l). (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) m1 + 1 times with respect to Xn introducing variables xnm1 ,xnm1 +1 ,...,xn1 . As a result, the set M j of points X k in Jan Awrejcewicz et al. 381 the (m1 + 1)th-dimensional space of variables xnm1 ,xnm1 +1 ,...,xn is obtained: X j m1 c,c,... ,c,c,xn X j m1 +1 j m1 , , (2.8) j m +1 p,c,... ,c,c,xn 1 X j 1 c,c,... , p,c,xn j X j c,c,... ,c, p,xn j 1 Note that it is advisable to take c = 1, p = 1. Step 2. Compute i = j 1. If i > m1 , then the set M j is modiï¬ed by introducing the variable x j m1 1 = c and the points X j m1 1 = X j m1 , X j m1 p,c,c,... ,c,c,xn j m1 (2.9) (2.10) Ë A new set M j of points {X k }, k = j m1 1, j m1 ,..., j, is obtained in the space of variables xnm1 1 ,xnm1 ,...,xn of the form X j m1 1 c,c,c,... ,c,c,xn X j m1 X j 1 c,c,c,... , p,c,xn j X j c,c,c,... ,c, p,xn j 1 j m1 1 j m p,c,c,... ,c,c,xn 1 , (2.11) , Construct a set of points attained through intersection of the hyperplane deï¬ned by the ( j m)th equation of system (2.6) with the straight lines going through points {X j ,X k }, k = j m1 1, j m1 ,..., j 1. For this purpose, i, or m1 + 1 if i > m1 , systems of equations aii+is xi+is = bi , j xi+is = xi+is s=1 j k xi+is xi+is (2.12) i + is n, Î»k , s = 1,2,...,l, 1 are solved. System (2.12) yields Î»k = bi j n l=1 aii+is xi+is j n k l=1 aii+is xi+is xi+is (2.13) On the economical solution of algebraic equations Table 2.2. Coordinates of points M j (for m1 < j n m2 ). xk1 xk x j 2 x j 1 xj x j+1 xm1 xm X k1 c c c c xk1 j xk1 j+1 k xm11 k xm1 Xk p c c c xk j xk j+1 k xm1 k xm X k+1 c k+1 xk c c xk+1 j xk+1 j+1 k+1 xm1 k+1 xm Â·Â·Â· Â·Â·Â· Â·Â·Â· Â·Â·Â· Â·Â·Â· Â·Â·Â· Â·Â·Â· X j 1 c c x j 2 c xj x j+1 j 1 j 1 Xj c c c j x j 1 xj j 1 j 1 Â·Â·Â· Â·Â·Â· x j+1 xm1 xm xm1 xm j 1 Then coordinates of the kth intersection point are computed: k Ëk xr = xr + xr xr Î» k , r = k,i,i + 1,...,i + m2 , r n. (2.14) Step 3. Assume j = i and return to Step 2. n m2 are reported. In Table 2.2, coordinates of the points of set M j for m1 < j Notice that for the repeated index combination, the following notation is introduced: k = j m 1 , m = j + m2 . 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 x1 ,...,xm2 is obtained. The space with obtained coordinates overlaps with the values of the unknown x1 ,...,xm2 occurring in the solution of system (2.6). Step 5. The order of the initial system is decreased ï¬rst by neglecting m2 equations and by putting the found values of the m2 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 m2 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 m2 as integers. Notice that by adding a zero diagonal to the compressed band matrix, parameter m2 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 Jan Awrejcewicz et al. 383 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 Â· 105 , 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 (m1 + 2)(m2 + 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 lm2 , l = 1,2,...,n/m2 2, j m2 + 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 km2 , k = 1,2,..., has been formulated, the corresponding set M j of points j = 2m2 + 1 + (k 1)m2 is created. The discussed algorithm for solving SLAE with a band matrix is also suitable for solving boundary value problems using the method of ï¬nite diï¬erences 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 diï¬ering from the canonical one, an artiï¬cial way of introducing ï¬ctitious equations for the mesh nodes to SLAE is as follows: xi (2.15) 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 n1 , n2 , n3 be the numbers of mesh nodes in directions ox1 , ox2 , ox3 , 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 n1 , 1 j n2 , 1 k n3 . 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 deï¬ned using ordered numbers of the unknowns corresponding to the nodes appearing in the pattern. On the economical solution of algebraic equations Numbering of unknowns is deï¬ned 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 coeï¬cients in all equations. The variant associated with the successive picking of nodes in directions ox1 , ox2 , ox3 and corresponding to the formula for the computation of unknown number n in the node {i, j,k} reads si, j,k = i + ( j 1)n1 + (k 1)n1 n2 . Formula (2.16) yields i1 = si, j,k1 si, j,k = n1 n2 , i2 = si, j 1,k si, j,k = n1 , i3 = si1, j,k si, j,k = 1, i4 = si, j,k si, j,k i5 = si+1, j,k si, j,k = 1, i6 = si, j+1,k si, j,k = n1 , i7 = si, j,k+1 si, j,k = n1 n2 . For the last band matrix parameters, one gets m1 = n 1 n 2 , m2 = m 1 , m = 2n1 n2 + 1 , n = n1 n2 n3 . (2.18) (2.17) (2.16) The formula used for the successive picking of nodes along directions ox3 , ox1 , ox2 , that is, the formula applied to compute an unknown number in the node {k,i, j }, is sk,i, j = j + (k 1)n3 + ( j 1)n1 n3 , where elements is of set I read i1 = n1 n3 , i4 i5 = 1, i 2 = n 3 , i6 = n 3 , i3 = 1, i7 = n 2 n 3 , (2.20) (2.19) and the band parameter matrix follows: m1 = n 1 n 3 , m2 = m 1 , m = 2n1 n3 + 1. (2.21) Finally, the last variant of the successive node picking in directions ox2 , ox3 , ox1 corresponds to the formula for computation of the unknown number in the node { j,k,i}, that is, s j,k,i = j + (k 1)n2 + (i 1)n2 n3 , where elements is of set I read i1 = n2 n3 , i4 i5 = 1, i 2 = n 2 , i6 = n 2 , i3 = 1, i7 = n 2 n 3 , (2.23) (2.22) Jan Awrejcewicz et al. 385 Table 2.3. Exchange of mesh nodes. Pattern Node i, j,k 1 i, j 1,k i 1, j,k i, j,k i + 1, j,k i, j + 1,k i, j,k + 1 Row Element 1 2 3 4 5 6 7 i, j,k 1 2 3 4 5 6 7 Variant of the nodes exchange k,i, j 3 1 2 4 6 7 5 j,k,i 2 3 1 4 7 5 6 and the band element matrix has the form m1 = n 2 n 3 , m2 = m 1 , m = 2n2 n3 + 1. (2.24) 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 n1 , n2 , n3 are diï¬erent, 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 nx . 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 n1 , n2 , n3 ) 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 deï¬ne 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 diï¬erential equations with respect to a few sought functions, then with the use of ï¬nite diï¬erence 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 noneï¬ective. 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. On the economical solution of algebraic equations In practice, both numeration methods mentioned above are applied. The ï¬rst 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 diï¬erence pattern scheme are computed ï¬rst with respect to each group of the same unknowns in accordance with the ï¬rst numbering order. Step 2. Elements of set Ik are computed, that is, the relative numbers of these unknowns for the cases when the unknowns lie on the main diagonal associated with the ï¬rst, the second,..., the kth group of the unknowns named the same. Step 3. A nonorder set I matching all sets Ik 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 l1 , l2 , m1 , m2 , m is computed. 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: ï¬rst elements of all blocks, second elements of all blocks, and so on; ï¬rst 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 deï¬ned through the condition of equality of the corresponding elements of sets Ik , I. Notice that in its essential part, the resulting band is ï¬lled 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 diï¬erential equations in partial derivatives of fourth order with variable coeï¬cients in two-dimensional space using the method of ï¬nite diï¬erences. Denote by n1 , n2 the numbers of mesh nodes in directions ox1 , ox2 , respectively. The derivatives in mesh node {i, j } are approximated through diï¬erence relations on the pattern consisting of 13 nodes (see Figure 2.1). The following coordinates of the local nodes {is , js }, s = 1,2,...,13, are introduced: i1 = i 2, j1 = j; i2 = i 1, j2 = j 1; i3 = i 1, j3 = j; i4 = i 1, j4 = j + 1; i5 = i, j5 = j 2; i6 = i, j6 = j 1; i7 = i, j7 = j; i8 = i, j8 = j + 1; i9 = i, j9 = j + 2; i10 = i + 1, j10 = j 1; i11 = i + 1, j11 = j; i12 = i + 1, j12 = j + 1; i13 = i + 2, j13 = j. (2.25) Jan Awrejcewicz et al. 387 i, j + 2 i 1, j + 1 i, j + 1 i + 1, j + 1 i 2, j i 1, j i, j i + 1, j i + 2, j i 1, j 1 i, j 1 i + 1, j 1 i, j 2 Figure 2.1. The pattern of 13 nodes. The input matrix in the compressed form (nonzero diagonals) is built owing to the ï¬rst numbering way for unknowns, and it describes a two-dimensional packet of dimension kn Â· kl , where k = 2 is the number of sought functions, n = n1 n2 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 blocks and possessing 39 nonzero diagonals. Local index s of the pattern node deï¬nes a position in the row of a coeï¬cient 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: i1,k â I1 , i2,k â I2 , k = 1,2,...,26. (2.26) For this purpose, the relative members of all unknowns corresponding to pattern nodes and appearing in equations for node {i, j } are computed: i1,k = lk l13 , i2,k = lk l14 , (2.27) where even k : lk = 2((is 1)n2 + js ), odd k : lk1 = 2((is 1)n2 + js ) 1, s = 1,2,...,13. The set I = I1 âª I2 is constructed and its elements are ordered in an increasing manner. The corresponding results for n2 = 5 are reported in the ï¬rst 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 n2 but is deï¬ned 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 I1 , I2 , and I, one obtains a rule of distribution of the row elements in the modiï¬ed matrix during construction of the band matrix in the compact form. The correspondence between the ordering numbers of nonzero coeï¬cients for the even and odd rows of the modiï¬ed matrix and rows of the band matrix is reported in the On the economical solution of algebraic equations Table 2.4. Comparison of obtained results (see text). N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 I1 20 19 12 11 10 9 8 7 4 3 2 1 I2 21 20 13 12 1 10 9 8 5 4 3 2 1 21 20 19 13 12 11 10 9 8 7 5 4 3 2 1 0 1 2 3 4 5 8 9 10 11 12 13 20 21 â â â â â 0 1 2 3 4 7 8 9 10 11 12 19 20 â â â â â k=1 2 3 5 6 7 8 9 10 12 13 14 15 16 17 18 19 20 21 23 24 25 26 27 28 30 31 â â â â â k=2 1 2 4 5 6 7 8 9 11 12 13 14 15 16 17 18 19 20 22 23 24 25 26 27 29 30 â â â â â last two columns of Table 2.4. The free positions in the row of the matrix are ï¬lled out by zeros. Notice that for arbitrary n1 , n2 , this correspondence is the same. As known, when solving numerically boundary value problems for PDEs using ï¬nite 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 Jan Awrejcewicz et al. 389 is, in this case, the equations dealt with are associated with the band matrix of irregular structure (the so-called cutting-oï¬ matrix). However, there exists a modiï¬cation of the elimination method for equations due to which it is possible to solve SLAE of a similar form. The cutting-oï¬ 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 Ii ). Here lies the fundamental diï¬erence in the application of the elimination method for the equations of SLAE using the cutting-oï¬ matrix. We comment on parameters m1 , m2 . Being the same for all equations, they are computed initially as m1 = max Ii (1) , m2 = max Ii li , (2.28) where l1 denotes the number of nonzero coeï¬cients in the ith equation. Therefore, the row of nonzero coeï¬cients and the row of indices for each equation are supplemented, if necessary, by zeros in the case of elements, and by Ii (1) = m1 , Ii li = m2 (2.29) in the case of indices, where li is the value of Ii satisfying the above requirement. As a result, SLAE with a band matrix of nonregular structure is obtained. Note that the formulation of the cutting-oï¬ matrix with the help of two packets is optimal with respect to the memory storage volume since for the index conservation, it is suï¬cient 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 ï¬nite diï¬erence 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 ï¬nally 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, ï¬niteness of the computational process, exactness) and iterational (minimal requirement of the storage volume) methods. We ï¬nally emphasize particularities of the SLAE matrix in each separate case, which can also be treated as an advantage of the method. 3. Numerical solution of a three-dimensional equation of elliptic type Many stationary processes of diï¬erent 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 ï¬ltration, and so forth. On the economical solution of algebraic equations 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 diï¬erential equations of elliptic type (linear or nonlinear). In the latter case, equations are linearized through diï¬erentiation along the parameter [7] or using the âsetupâ method [2]. Hence, the problem is ï¬nally 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 ï¬nite diï¬erences to elliptic equations results in SLAE with a band-type matrix. In the case of PDEs, the matrix (SLAE) possesses the cutting-oï¬ band, with only a few nonzero elements (see Section 2). We brieï¬y model the problem of stationary heat distribution in a certain volume G with surface Î of three-dimensional space x = (x1 ,x3 ). 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 W = K grad V , (3.1) where K = K(x) is the heat transfer coeï¬cient. 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 Wn 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, Wn dS, (3.2) should be equal to the heat amount of Ï(x)dU (3.3) appearing in volume U, that is, Wn dS = Ï(x)dU. (3.4) Using the Gauss formula form S Wn dS = U div, balance equation (3.4) is transformed to the div W Ï(x) dU = 0. (3.5) As volume U is arbitrary, therefore, if Ï(x) and div W are continuous functions of point x = (x1 ,x3 ), then (3.5) yields W = Ï(x). (3.6) 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: LV = div(K grad V ) = Ï(x), or in the equivalent form, (3.7) K(x)V (x) = Ï(x), (3.8) where K is the function of point x = (x1 ,x3 ). Heat transfer equation (3.8) is obtained under an assumption of isotropy of the heat transfer process. If a heat coeï¬cient depends on direction and is a tensor (an isotropic medium), then, instead of (3.8), the following equation is taken: KÎ±Î² V (x),xÎ² Î±1 Î²=1 = Ï(x). (3.9) If KÎ±Î² â¡ O and Î± = Î², then (3.9) reads KÎ± V (x) = Ï(x). (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: K dV = Ã¦V + g(x), dn x â Î,Ã¦ = Ã¦(x) > 0. (3.11) Consider the boundary value problem for the equation of the stationary heat transfer in a rectangular parallelepiped of a three-dimensional space: Â¯ G= 0 xÎ± lÎ± , Î± = 1,2,3 . (3.12) In this case, the method of ï¬nite diï¬erences is used; the construction of the diï¬erence 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: KÎ± (X)V (x) = Ï(x), Â¯ x â G; xÎ± xÎ± = l Î± , (3.13) (3.14) KÎ± (x)V (x) = Ã¦ (x) V (x) V0 (x) , KÎ± (x)V (x) = Ã¦+Î± (x) V (x) V0 (x) , On the economical solution of algebraic equations where KÎ± (x) is the heat transfer coeï¬cient in direction xÎ± , xÂ±Î± (x) is the heat exchange coeï¬cient with surrounding medium on walls xÎ± 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) satisï¬es the functional minimum: I(u) = + G 3 lÎ² 2 kÎ± V 2V Ï dx lÎ³ 0 0 V 2 2g V xÎ± =0 + V 2 2g+Î± V xÎ± =lÎ± dxÎ² dxÎ³ , Î± = Î² = Î³, (3.15) where gÂ±Î± (x) = Ã¦Â±Î± (x)V0 (x). Indeed, variating functional (3.15) with respect to V , V , integrating by parts the terms with Î´V , 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: Â¯ Ï= Â¯ ÏÎ± , (3.16) Â¯ where ÏÎ± = {xÎ±i , i 1,2,...,NÎ± 1,NÎ± }âmesh on interval [0, lÎ± ]. The set of interval mesh nodes is denoted by Ï= ÏÎ± , (3.17) where ÏÎ± = {xÎ±i , i = 1,2,...,NÎ± 1}. We introduce the following notations: hÎ±i = xÎ±i xÎ±i1 , Ë hÎ±i = ï£´0, ï£± ï£´0, ï£´ ï£´ ï£´ ï£² ï£´ ï£´ ï£´ ï£³0, i = 1,2,...,NÎ± , 5Î±1 , i 5hÎ±i1 + hÎ±i , i = 1,2,...,NÎ± 1, 5hÎ±NÎ± , i = NÎ± , vi, j,k = v x1i, x2 j, x3k , (3.18) , , . i, j,k i, j,k i, j,k iÂ±1, j,k i, j,Â±1,k i, j,kÂ±1 k1 iÂ±0,5, j,k 5 k1 5 + k1 + i, j,Â±0,5,k i, j,kÂ±0,5 k3 5 k3 + k3 Jan Awrejcewicz et al. 393 All 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 diï¬erence approximations, for instance, by v (x1i 0,5h1i x2 j, x3k ) = (vi, j,k vi1, j,k )/h1i . In what follows, the function with n variables vi, j,k is obtained: N 1 N 2 N3 In (v) = i=1 j =0 k=0 k1 i0,5, j,k 2Ë Ë vi, j,k vi1, j,k h2 j h3k h1i 2Ë Ë vi, j,k vi, j 1,k h1i h3k h2 j 2Ë vi, j,k vi, j,k1 h1i h2i h3k N 1 N2 N 3 i=0 j =1 k=0 N 1 N2 N 3 k3 i, j 0,5,k i, j,k0,5 i=0 j =0 k=1 N 1 N2 N 3 i=0 j =0 k=0 3 NÎ² NÎ³ Ë Ë Ë Ïi, j,k vi, j,k h1i h2 j h3k Ã¦ v xÎ±0 ,xÎ² j ,xÎ³k j =0 k=0 2g v xÎ±0 ,xÎ² j ,xÎ³k + Ã¦+Î± v xÎ±NÎ± ,xÎ² j ,xÎ³k Ë Ë 2g+Î± v xÎ±NÎ± ,xÎ² j ,xÎ³k hÎ² j hÎ³k , Î± = Î² = Î³, n= NÎ± + 1 . (3.19) As known, the minimum of this function is achieved at the point where its partial derivatives are equal to zero with respect to vi, j,k . We construct a system of diï¬erence equations of order n to deï¬ne values vi, 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 G. Terms (3.19) including vi, j,k are reported below: S1 = + k1 i0,5, j,k vi, j,k vi1, j,k h1i k1 vi+1, j,k vi, j,k h1i+1 Ë Ë h2 j h3k i, j 0,5,k vi, j,k vi, j 1,k h2 j vi, j+1,k vi, j,k h2 j+1 Ë Ë h1i h3k (3.20) vi, j,k vi, j,k1 h3k Ë 1i h2 j h3k . Ë Ë 2Ïi, j,k vi, j,k h + i, j,k0,5 k3 vi, j,k+1 vi, j,k h3k+1 Ë Ë h1i h2 j On the economical solution of algebraic equations By diï¬erentiating S1 with respect to v(i, j,k) and comparing its derivative to zero, the following diï¬erence equation is obtained at node {i, j,k} of the mesh owing to the division by 2h1i h2 j h3k : k1 vi+1, j,k vi, j,k /h1i+1 k1 Ë h1i i0,5 j,k vi, j,k vi1, j,k /h1i i, j 0,5,k vi, j+1,k vi, j,k /h2 j+1 Ë h2 j vi, j,k vi, j 1,k /h2 j ) (3.21) vi, j,k vi, j,k1 /h3k k N3 . = Ïi, j,k , N1 1, 1 Ë h3k j N2 1, 1 (2) Node of the parallelepiped wall. i 1 j N2 1, 1 k N3 1. (3.22) Terms (3.19) (including variable vi, j,k for i = 0) follow: S2 = k1 2Ë Ë vi+1, j,k vi, j,k h2 j h3k h1i+1 i, j 0,5,k vi, j,k vi, j 1,k h2 j vi, j+1,k vi, j,k h2 j+1 Ë Ë h1i h3k (3.23) vi, j,k vi, j,k1 vi, j,k+1 vi, j,k + h3k h3k+1 2 Ë Ë Ë Ë Ë 2Ïi, j,k vi, j,k h1i h2 j h3k + Ã¦1 vi, j,k 2g1 vi, j,k h2 j h3k . i, j,k0,5 k3 Ë Ë h1i h2 j By diï¬erentiating S2 with respect to vi, j,k and comparing the corresponding derivative Ë Ë Ë to zero, the following diï¬erence equation (after division by 2h1i h2 j h3k ) is obtained at the parallelepiped node x1 = 0: k1 vi+1, j,k vi, j,k /h1i+1 Ã¦1 vi, j,k ) Ë h1i vi, j+1,k vi, j,k /h2 j+1 Ë h2 j Ë h3k i, j 0,5,k vi, j,k vi, j 1,k /h2 j (3.24) vi, j,k vi, j,k1 /h3k = Ïi, j,k Ë . h1i g1 Diï¬erence equations at other parallelepiped wall nodes are constructed analogously. Jan Awrejcewicz et al. 395 (3) Mesh nodes belong to the parallelepiped i = j 1 k N3 1. (3.25) By proceeding in a way similar to that applied in the previous case, the following diï¬erence equation is obtained: k1 vi+1, j,k vi, j,k Ë Ã¦1 vi, j,k h1i h1i+1 i, j+0,5k vi, j+1,k vi, j,k /h2 j+1 Ã¦2 vi, j,k Ë h2 j (3.26) vi, j,k vi, j,k1 /h3k = Ïi, j,k Ë Ë . h1i h2 j g1 g2 Ë h3k (4) A node of the mash coincides with the parallelepiped vertex i N1 , j N2 , k N3 . (3.27) Proceeding again as in the previous case, one arrives at the following corresponding equations: for i = j = k k1 vi+1, j,k vi, j,k /h1i+1 Ã¦1 vi, j,k Ë h1i k3 vi, j+1,k vi, j,k /h2 j+1 Ã¦2 vi, j,k Ë h2 j (3.28) vi, j,k+1 vi, j,k /h3k+1 Ã¦3 vi, j,k + Ë h3k g1 g2 g3 = Ïi, j,k Ë Ë Ë ; h1i h2 j h3k for i = N1 , j = N2 , and k = N3 , Ã¦1 vi, j,k k1 i0,5, j,k i, j,k+0,5 vi, j,k vi1, j,k /h1i ) i, j 0.5,k Ë h1i Ã¦2 vi, j,k vi, j,k vi, j 1,k /h2 j ) (3.29) Ë h2 j i, j,k0,5 Ã¦3 vi, j,k k3 vi, j,k vi, j,k1 /h3k Ë h3k g+1 g+2 g+3 = Ïi, j,k Ë Ë Ë . h1i h2 j h3k On the economical solution of algebraic equations Thus, all possible cases occurring while formulating the diï¬erence equations on the mesh with variable step for the third boundary value problem have been discussed. We discuss brieï¬y the peculiarities involved in the formulation of diï¬erence equations in the case when the following Dirichlet conditions are given on some parallelepiped walls: v(x)|xÎ± =0,IÎ± = vu (x), (3.30) where vu (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: vi, j,k = vu . (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 diï¬erence 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 diï¬erence 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 x1 = 1 /2, then we take Ë h1N1 = h1N1 , k1 1 N +0,5, j,k = k1 N1 0,5, j,k (3.32) In what follows, in the diï¬erence equations formulated for the nodes of the plane xÎ± = Î± /2, accounting for the symmetry property v + h1N1 ,x3 = v h1N1 ,x3 , (3.33) the coeï¬cient standing by the unknown vi1, 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 diï¬erence 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 ï¬rst 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 h1i , h2 j , h3k ; formulas for kÎ± (x), gÂ±Î± (x), Ï(x) computations, and vu (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). i, j,k i0,5, j,k Having this information, the operator coeï¬cients k1 , k1 are computed and the matrix of diï¬erence equations is formulated due to formulas (3.21)â(3.29). Analogous operations are carried out with (x), k3 (x). Jan Awrejcewicz et al. 397 In what follows, the matrix of coeï¬cients of the diï¬erence 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 diï¬erence 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 k1 = 1 + 2 , Ë Ë Ë 3c x1 + x2 + x3 2 , = 1 + Ë Ë Ë 3c 3/4 x1 x2 x3 2 , k3 = 1 + Ë Ë Ë 3c 1/2 x1 x2 + x3 (3.34) where Ë xÎ± = xÎ± The exact solution reads vu (x) = and hence 2 Ï(x) = x3 2 + Î± = 1,2,3. (3.35) Ë Ë 2 1/2 x1 x2 x3 , 2 (3.36) 4 Ë Ë 3c 3/8 + x1 x2 1+ 2 Ë Ë Ë 3c 1/2 x1 x2 + x3 + 2x3 x3 1/2 , (3.37) 1 Ë Ë x1 x2 2 2 3 Ë g2 x1 ,0,x3 = x1 + 2 , 4 Ë 3c 1/2 x1 x2 x3 /2 g3 x1 ,0,x3 g+3 x1 ,1 = 1 Ë Ë x1 x2 2 3 3 + Ë Ë 2 2c 3/4 x1 x2 assuming that Ã¦Â±Î± (x) = 1. We introduce the Dirichlet condition on the wall x1 = 0. The symmetry of solution is taken into account with respect to the plane x1 = 1/2, x2 = 1/2. The solution is deï¬ned in Â¯ the space G = {0 xÎ± 1, Î± = 1,2,3}. On the economical solution of algebraic equations Table 3.1. Numerical solution of the model problems at point: x1 = 0.5, x2 = 0.5, x3 = x3k . k 1 2 3 4 5 x3k 0 0.25 0.5 0.75 1.0 3Ã3Ã5 0 0.0112 0.504 0.120 0.223 Mesh 5Ã5Ã5 0 0.0120 0.0523 0.122 0.0224 9Ã9Ã9 0 0.0147 0.0599 0.136 0.243 Exact solution 0 0.0156 0.0625 0.141 0.250 Table 3.2. Numerical solution of the model problems at point: x1k = x2k = x3k . k 1 2 3 4 x3k 0 0.25 0.5 0.75 3Ã3Ã3 0 â 0.378 â Mesh 5Ã5Ã5 0 0.0038 0.0375 0.114 9Ã9Ã5 0 0.0060 0.0444 0.127 Exact solution 0 0.0068 0.0469 0.132 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 coeï¬cients in the considered case is high, that is, 1 kÎ± (x) 1024, Î± = 1,2,3. (3.38) Notice that for c the numerical solution with accuracy greater than ï¬ve 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 testiï¬ed through the results included in Table 3.2. The following parameters are taken: k1 = = k3 = 207, and the following boundary conditions are applied: v 0 ,x3 k3 v,x3 x1 ,0 = 10(v 100), v x1 ,0,x3 k3 v,x3 x1 ,l3 = 10v. (3.40) Ï(x) (3.39) The following symmetry solution with respect to the planes is taken into account: x1 = l1 , 2 x2 = l2 ; 2 l1 = l2 = 1m, l3 = 0.02m. (3.41) Jan Awrejcewicz et al. 399 Table 3.3. Numerical solution to the problem without source. Mesh k 1 2 3 4 5 x1k 0 0.125 0.25 0.375 0.5 5Ã5Ã5 0 9.17 15.3 18.8 19.9 9Ã9Ã5 0 9.2 15.3 18.8 20.0 Recipes â x1k = x2k x3 = 3 â â 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 eï¬cient 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: ï£± ï£´ ï£²4 Â· 106 , k1 = = k3 = 1, Ï(x) = ï£´ ï£³0, for x1 = x2 = otherwise, , x3 = 3 , (3.42) where the attached boundary conditions are as follows: v 0 ,x3 v,x3 x1 ,0 v x1 ,0,x3 = 10v, v,x3 x1 ,l3 = 0; (3.43) = 1, = 0.02. The following symmetry condition with respect to the planes x1 = x2 = (3.44) 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 x1 , x2 in the ï¬rst 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 diï¬erences are observed at the node with the source, but the diï¬erence is less than 1% for two neighboring partitions. On the economical solution of algebraic equations Table 3.4. Numerical results of the problem with point source. k k3k 1 0 2 0.25 3 0.5 4 0.75 5 1 3Ã3Ã5 264.5 269.5 284.4 310.9 351.7 5Ã5Ã5 300.9 305.7 320.6 347.2 388.0 Mesh 9Ã9Ã5 321.0 325.9 340.9 367.5 408.3 17 Ã 17 Ã 5 324.3 329.2 344.2 370.8 411.6 Recipes â x1 = 0.5 x2 = 0.5 â â Table 3.5. Numerical results of the problem with point source. k k1k 1 0 2 0.25 3 0.5 4 0.75 3Ã3Ã5 0 â 36.6 â 5Ã5Ã5 0 12.05 38.9 90.0 Mesh 9Ã9Ã5 0 12.25 39.6 92.0 17 Ã 17 Ã 5 0 12.31 39.8 92.5 Recipes â x2k = x1k x3k = 3 â (4) Comparison of solutions obtained through various methods. Finally, we present the results of computation through the ï¬nite diï¬erence 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 ï¬rst boundary value problem for the Laplace equation for the cube with ribs of unit length are reported. The boundary conditions are as follows: v 0.5 x3 = 2, v x1 , 0.5,x3 v x1 , 0.5 v 0.5 ,x3 = 1, v x1 ,0.5,x3 v x1 ,0.5 = 0. (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 v 0,5 ,x3 = 10, v + 5v x2 =Â±1 v + 5v v,x3 + 5v x1 =0,5 x3 =Â±1 = 0. (3.46) Notice that when numerical solutions are constructed, the symmetry with respect to two planes, x2 x3 is taken into account. The reported results show that the solution obtained through the ï¬nite diï¬erence method during the mesh step decrease converges to the analytical solution, and its error is small for the appropriate partition. Jan Awrejcewicz et al. 401 Table 3.6. Solution to the ï¬rst boundary value problem for Laplace equation. x1k 0.375 0.250 0.125 0 0.125 0.250 0.375 Boundary element method N = 12 N = 24 1.637 1.472 1.044 0.979 0.678 0.661 0.5 0.5 0.478 0.472 0.597 0.566 0.770 0.770 Finite diï¬erence method 5 Ã 3 Ã 3 9 Ã 5 Ã 5 17 Ã 9 Ã 9 â 1.413 1.426 0.926 0.953 0.963 â 0.653 0.657 0.5 0.5 0.5 â 0.471 0.472 0.545 0.555 0.558 â 0.740 0.746 Analytical solution 1.430 0.967 0.659 0.5 0.472 0.560 0.748 Table 3.7. Solution of the Laplace equation with hybrid boundary conditions (x3k = x2k ). Boundary element method x1k 0.25 Finite element method 5Ã3Ã3 7.155 4.724 3.821 2.778 â 2.105 â 9Ã5Ã5 7.234 4.804 3.835 2.825 2.645 2.094 1.202 17 Ã 9 Ã 9 7.257 4.829 3.841 2.839 2.655 2.091 1.186 Analytical solution 7.259 4.837 3.843 2.844 2.658 2.089 1.180 x2k 0 0 0.5 0 0.25 0 0.75 N = 24 7.387 4.827 3.745 2.816 2.612 2.0 1.050 N = 48 7.282 4.84 3.843 2.843 2.658 2.073 1.144 The results of the numerical experiments yield the conclusion that our proposed method and the associated subroutines allow eï¬cient 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. 4. 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 stiï¬ened 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 coeï¬cient Âµ are functions of point coordinates x = (x1 ,x3 ), or, On the economical solution of algebraic equations x2 x3 x1 Figure 4.1. Investigated shell. according to the theory of small elastic-plastic deformations, they depend on the stressstrain material state at the considered point. We introduce the following notation: x1 , x2 , x3 ârectangular system of coordinates (Figure 4.1); Î±, b, hâshell dimension and thickness, respectively; k1 , âcurvatures; u, v, wâcomponents of displacements of mean surface point; qâtransversal load intensity. Following the Kirchhoï¬-Love hypothesis, the components of shell and rib deformation read eii = Îµii + Ïii x3 where Îµ11 = u k1 w + 2 w,x Îµ22 = v w + 2 , 2 2 w , 2 (i = 1,2), e12 = Îµ12 + Ï12 x3 , (4.1) (4.2) Îµ12 = U + v + w w , Îµ11 , Îµ22 , Îµ12 âlength and stream deformations of the mean surface; Ï11 = wx1 x1 , Ï22 = wx2 x2 , Ï12 = wx1 x2 , (4.3) where Ï11 , Ï22 , Ï12 are variations of curvature and torsion. In accordance with the Hook principle for a ï¬at-strain state and (4.1), the coupling between stressed Ïik and deformation eik is presented in the following form: Ï11 = 2G Îµ11 + Ï11 x3 + Âµ Îµ22 + Ï22 x3 1v Ï12 = G Îµ12 + 2Ï12 x3 , (1 ââ 2), (4.4) where, here and further, (x ââ y) denotes the change of index x to y. Jan Awrejcewicz et al. 403 Integration of (4.4) with respect to x3 yields the following expressions for forces and movements: T11 = c00 Îµ11 + c10 Îµ22 + c01 Ï11 + c11 Ï22 , T22 = c00 Îµ11 + c00 Îµ22 + c11 Ï11 + c01 Ï22 , T12 = 0.5 c00 c10 Îµ12 + c01 c11 Ï12 , M11 = c01 Îµ11 + c11 Îµ22 + c02 Ï11 + c12 Ï22 , M22 = c11 Îµ11 + c01 Îµ22 + c12 Ï11 + c02 Ï22 , M12 = 0.5 c01 c11 Îµ12 + c02 c12 Ï12 , where k 2G(x)Âµi (x)x3 3 dx 1 Âµ(x) h/2 h/2 (4.5) (4.6) (4.7) (4.8) (4.9) (4.10) cik x1 = (i 1; k 1,2). (4.11) The deformation energy of the considered system reads â = 0.5 Î± Î± Adx1 dx2 + H1 + P 1 Adx1 dx2 + H1 + P 1 x2 =0 dx1 (4.12) x1 =0 dx2 x2 =0 dx1 + + H2 + H2 + P where A = T11 Îµ11 + T22 Îµ22 + T12 Îµ12 + M11 Ï11 + M22 Ï22 + 2M12 Ï12 , Â± Â± 2 2 Â± 2 Â± 2 H1 = D1 Îµ11 + AÂ± w x1 + B1 v x1 + C1 w x2 , 1 Â± Â± 2 2 Â± Â± 2 H2 = D2 Îµ22 + AÂ± w + B2 u2 2 x2 + C2 w x2 , ,x 2 Â± 2 Â± P1 = Î±Â± w + Î²1 w2 , 1 (4.13) P = Î±[w ]2 + Î²w2 ; Â± Â± Â± D1 , AÂ± , B1 , C1 ârib stiï¬ness on elongation (compressing), bending in a vertical plane, 1 bending in a horizontal plane, and torsion, respectively; Î±Â± , Î±âstiï¬ness of hinge rib1 Â± shell, rib-rib interaction, respectively, Î²1 , Î²âelasticity of support under the rib; [w ]â jump of derivative w on the line x1 = 0. We introduce the following functional: I(u,v,w) = â qw dx1 dx2 Î± 0 qw dx1 dx2 . (4.14) On the economical solution of algebraic equations Variating (4.14) with respect to u, v, w and carrying out transformations with the help of integration by parts, one gets Î´I = 0 Ë A dx1 dx2 + Î± 0 Ë A dx1 dx2 x2 =0 R1 |x2 =0,b + S + Î²1 wÎ´w + Î± w Î´w 1 1 dx1 Î± 0 b + R1 |x2 =0,b + S Î²1 wÎ´w + Î±+ w Î´w 1 1 x2 =0 dx1 x1 =0 (4.15) dx2 + Q1 R2 |x1 =,0 + R2 |x1 =0,Î± + S + S Î± w Î´ w + Î²wÎ´w 2 2 x1 =,0 + x2 =0 + Q2 + Q2 + x1 =0 + Q1 x1 =0,Î± x2 =0,b x2 =b 2M12 Î´w x1 =,0 2M12 Î´w x1 =0,Î± x2 =0,b x2 =0,b where Ë A = T11 + T12 Î´u T12 + T22 Î´v M11 x1 + M22 + 2M12 x2 + k1 T11 + T22 + T11 w + T12 w + T12 w + T22 w q Î´w, R1 = T12 Î´u + T22 Î´v + M22 + 2M12 + T22 w Î´w M22 Î´w , Â± SÂ± = D1 Îµ11 1 Î´u + Â± B1 v x1 x1 Î´v + Â± Â± k1 D1 Îµ11 D1 Îµ11 w (4.16) Â± + AÂ± w x1 Î´w C1 w x2 Î´wx2 , 1 Â± Â± Â± Q1 = D1 Îµ11 Î´ B1 v x1 Î´v + Â± B1 v x1 Î´v Â± Î´w + AÂ± w x1 Î´w + C1 w x2 Î´w 1 Â± + D1 Îµ11 w AÂ± w x1 1 1 ââ 2, u ââ v, x1 ââ x2 . We consider some of the relations that follow from the previous considerations. Diï¬erential equations governing the equilibrium state of a nonhomogeneous shell with respect to displacements are obtained by comparison to zero of the coeï¬cients standing by variations of u, v, w in doubled integral (4.15), and by taking into account relations (4.2), (4.6), and (4.9). However, they are not explicitly given owing to their complexity. Boundary conditions for the side x2 = 0 follow from the examinations of the under-integral expressions. If the displacements on the contour are not given, then the boundary conditions are deï¬ned by comparison to zero of the coeï¬cients standing by variations of u, v, w, w . Jan Awrejcewicz et al. 405 Conditions for the coupling of solutions on the line x1 = 0 are derived from the following expression: R2 |x1 =0 R2 |x1 =0 + S + S Î± w Î´ w + Î²wÎ´w |x1 =0 2 2 = T11 + B2 u Î´u + T12 D2 Îµ22 Î´v + T13 D2 Îµ22 D2 Îµ22 w + M11 C2 w x2 + where T1i = T1i + M11 C2 w x2 + A2 w + Î²w Î´w (4.17) Î± w Î´w Î± w Î´w x1 =0 T1i x1 =0 (i = 1,2,3), + B2 = B 2 + B 2 , T13 = M11 + 2M12 + T11 w , + D2 = D2 + D2 , (4.18) A2 = A2 + A+ . 2 This means that the following conditions should be satisï¬ed on the line of the solution coupling: T11 = B2 u , T12 = D2 Îµ22 , T13 = D2 Îµ22 D2 Îµ22 w Î± + M11 = C2 w x2 w = M11 M11 + + A2 w , + Î²w, C2 w x2 + C2 w x2 (4.19) + C2 w x2 , as well as the conditions on the continuity of displacements [u] = [v] = [w] = 0. Taking = Îµ22 = u = v = T11 w = 0 in (4.19), the conditions for the solution coupling for a plate are obtained [5]. We determine agreement conditions in the shell angle x1 x2 = 0. For this purpose, the forms related to an angle are considered in (4.15), that is, Q1 Q2 2M12 Î´w x1 =0 x2 =0 = D1 Îµ11 + B2 u Î´u + + B1 v x1 + D2 Îµ22 Î´v 2M12 D1 Îµ11 w + A+ w x1 1 + D2 Îµ22 w + A+ w 2 Î´w (4.20) + B1 v x1 Î´v B2 u Î´u A+ w x1 C2 w x2 Î´w 1 + A+ wx2 x2 + C1 w x2 Î´w = 0. 2 One may derive now, for instance, the following compatibility conditions in the angle: 2M12 + A1 w x1 + A2 w D1 Îµ11 w D2 Îµ22 w = 0; B2 u B1 v x1 A+ w + C1 w x2 = 0. A+ w x1 + C2 w x2 1 2 D1 Îµ11 + B2 u D2 Îµ22 + (B1 v x1 (4.21) On the economical solution of algebraic equations Notice that if the ribs are excluded, then only one compatibility condition M12 = 0 remains. We derive compatibility conditions for the point x1 where the change of boundary conditions occurs. They are obtained from the terms of (4.15) with respect to the point x1 x2 and the coeï¬cients with index 2 are equal to zero (homogeneous shell): + Q1 Q1 x=0 x2 =0 = D1 Îµ11 D1 Îµ11 Î´u + B1 v x1 Î´v + D1 Îµ11 w + D1 Îµ11 w + A+ w x1 1 Î´v B1 v x1 Î´v + + A w x1 Î´w A+ w x1 Î´w 1 1 A w x1 1 Î´w + C1 w x2 C1 w x2 Î´w (4.22) = 0. As an example, we consider the compatibility condition for a hinged support of the shell side x2 = 0 with diï¬erent stiï¬ness: D1 Îµ11 D1 Îµ11 B1 v x1 C1 w x2 C1 w x2 = 0. (4.23) 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: F = T11 , F,z1 x2 = T12 , F x1 = T22 , (4.24) and we derive the equations in a hybrid form with respect to force F and deï¬ection w. For this purpose, the deformation continuity equation and the third equation of (4.15) are applied to yield Îµ11 + Îµ22 x1 Îµ12 x2 + â2 w + 0,5L(w,w) k M11 x1 + M22 + 2M12 x2 + â2 F + L(w,F) q = 0. k By solving (4.6) with respect to Îµik , one gets Îµ11 = AF + A F x1 + Bw x1 = B w , Îµ22 = A F + AF x1 + B w x1 + Bw , where A = A1 + A2 , A = A1 A2 , 0.5 A1 = , c00 + c10 B1 = c01 + c11 A1 , B = B 1 + B2 , B = B1 B2 , 0.5 A2 = , c00 c10 B2 = c01 c11 A2 . Îµ12 = 4A2 F x2 + 4B2 w x2 , (4.26) (4.25) (4.27) Jan Awrejcewicz et al. 407 Taking into account (4.27) in (4.9), one obtains M11 = BF + B F x1 Cw x1 C w , M22 = B F + BF x1 C w x1 Cw , M12 = 2B2 F x2 (C C)w x2 , where C = c02 B1 B2 , B1 = c01 + c11 B1 , C = c12 B1 B2 , B2 = c01 c11 B2 . (4.28) (4.29) Equations (4.25) with an account of (4.27), (4.28) read A F x1 + AF + Bw x1 + B w + AF x1 + A F + B w x1 + Bw 4A2 F x2 + 4B2 w x2 x1 2 x2 + âk w + 0.5 L(w,w) = 0; BF x1 + B F C w x1 Cw + B F x1 + BF Cw x1 C Ï 2 2B2 F x2 (C C)w x2 x2 + â2 F + L(w,F) q = 0. k (4.30) x1 In what follows, the algorithm for solving the stability problem of physically nonlinear ï¬exible 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 {x1i j ,x3k } is introduced. Assuming in the beginning that the shell is unloaded, modulus Gi jk of the shell nodes is introduced. Step 2. At the nodes of mesh surface x3 coeï¬cients CÎ± Î² (4.11) are computed. Then, substituting the partial derivatives by diï¬erence expressions with error Q(h21 + h22 ), inx x stead of the system of equations (4.30), the system of nonlinear algebraic equations with respect to wi j , Fi j , that is, deï¬ection 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(x1 x2 ) deï¬ned with the help of the Newton or diï¬erentiation along the parameter methods [7], and then the elimination method is applied. Step 3. New values of Gi 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 E= 9KG , 3K + G Âµ= 1 3K 2G . 2 3K + G (4.31) On the economical solution of algebraic equations 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(x1 ,x3 ), and for a body homogeneous before deformation, K = K0 = const. The shear modulus is deï¬ned through the formula G= 1 Ïi ei , 3 ei (x) (4.32) and G is called the cutting modulus. In order to compute the shell, dependence Ïi (ei ) should be explicitly given (intensity of stress Ïi versus strain (ei )). Some of the diagrams Ïi (ei ) 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 + hq 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 k1 = = 18, Î» = a/h = 50, and q = const [4] has been considered. We apply the Mises ï¬ow condition (Ïi = Ïs ) and the dependence Ïi (ei ) taken in the form Ïi = ÏS 1 exp ei eS (4.33) where es = Ïs /3Ï0 is the artiï¬cial intensity of the ï¬ow deformation. The problem is considered symmetric with respect to x1 , x2 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 x2 = 0 is F x1 F x2 whereas for the deï¬ection, w wx2 = 0 or M22 = 0. Consider a few computational results. In Figure 4.2, dependencies of load versus center deï¬ection 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 inï¬uence of the boundary conditions on the distribution of plasticity zones on the upper shell surface x3 = h/2 for the loads corresponding to deï¬ection w 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. Jan Awrejcewicz et al. 409 q 50 1.0 w Figure 4.2. Load versus shell center deï¬ection (see text for more details). Figure 4.3. Inï¬uence of boundary conditions on plasticity zone distribution (see text for more details). M11 1 0 0.25 x2 Figure 4.4. Moment M11 versus x2 on shell side x1 = 0 (see text). Figure 4.4 displays graphs of moments M on the side x1 = 0 for the following cases: 2âsolid curve; 3âdashed curve; 4âdash-dotted curve. On the economical solution of algebraic equations Deï¬ection in the shell center achieved in all cases is w = 0.4. In Figures 4.2 and 4.3, digits 1â4 should be considered simultaneously. For the given computational model, contrary to [8], a jump of the moments is observed, which is close to real observations. The results of Section 4 make the elimination method for equations eï¬ective and also for solving the complicated PDEs like nonlinear equations governing elastic-plastic problems of shells with ï¬nite deï¬ections. Here, on each loading step, one has to solve many times a large SLAE, which, without any doubt, exerts signiï¬cant impact on the chosen computational algorithm. References [1] [2] [3] [4] [5] [6] [7] [8] J. Awrejcewicz and V. A. Krysko, Nonclassical Thermoelastic Problems in Nonlinear Dynamics of Shells, Scientiï¬c Computation, Springer-Verlag, Berlin, 2003. V. I. Feodosâev, On a method of solution of the nonlinear problems of stability of deformable systems, J. Appl. Math. Mech. 27 (1963), 392â404. G. E. Forsythe and C. B. Moler, Computer Solution of Linear Algebraic Systems, Prentice-Hall, New Jersey, 1967. V. A. Krysko, Nonlinear Statics and Dynamics of Non-Homogeneous Shells, Saratov State University Press, Saratov, 1976. A. A. Samarskii and V. B. Andreev, Diï¬erence Methods for Elliptic Equations, Nauka, Moscow, 1976. A. S. Saninskiy, Direct Method of Solution for a System of Linear Algebraic Equations with Calculated Matrix, Saratov State University Press, Saratov, 1877. V. I. Shalashvili and E. B. Kusnetsov, Application of Solution along the Parameter Method and on the Best Parametrization, Editorial URSS, Moscow, 1999. L. F. Vakhlaeva and V. A. Krysko, Stability of ï¬exible shallow rectangular shells with various boundary conditions on their sides, Building and Architecture. IV, VUZ Press, Kazan, 1984, pp. 21â25. Jan Awrejcewicz: Department of Automatics and Biomechanics, Technical University of Lodz, 90-924 Lodz, Poland E-mail address: awrejcew@mail.p.lodz.pl Vadim A. Krysko: Department of Mathematics, Saratov state University, 410054 Saratov, Russia E-mail address: tek@sun.ru Anton V. Krysko: Department of Mathematics, Saratov state University, 410054 Saratov, Russia
http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png
http://www.deepdyve.com/lp/hindawi-publishing-corporation/on-the-economical-solution-method-for-a-system-of-linear-algebraic-qgk5kpOazQ
/lp/hindawi-publishing-corporation/on-the-economical-solution-method-for-a-system-of-linear-algebraic-qgk5kpOazQ