“Woah! It's like Spotify but for academic articles.”

Instant Access to Thousands of Journals for just $40/month

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

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 first-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 first-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 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 Ax = b, or finding solutions to a few variants of problem Ax = 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 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 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. 2. 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 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 fixed. 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 defined 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 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 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 modified 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 defined 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 first 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 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: 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 defined 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 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 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 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 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 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. On the economical solution of algebraic equations 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 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 first, 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: 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 Ik , 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 n1 , n2 the numbers of mesh nodes in directions ox1 , ox2 , 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 {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 first 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 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: 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 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 n2 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 I1 , I2 , 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 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 filled 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 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 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-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 Ii ). 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 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 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 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-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. 3. 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. 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 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 = (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 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 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 coefficient 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 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: 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 coefficient in direction xα , x±α (x) is the heat exchange coefficient 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) satisfies 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 difference 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 difference equations of order n to define 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 differentiating S1 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 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 differentiating S2 with respect to vi, j,k and comparing the corresponding derivative ˜ ˜ ˜ to zero, the following difference 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 Difference 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 difference 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 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: 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 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 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 difference equations formulated for the nodes of the plane xα = α /2, accounting for the symmetry property v + h1N1 ,x3 = v h1N1 ,x3 , (3.33) the coefficient 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 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 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 coefficients k1 , k1 are computed and the matrix of difference 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 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 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 defined 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 coefficients 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 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: 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 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:   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 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. 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 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: 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 finite difference 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 first 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 difference 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 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. 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 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 = (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 Kirchhoff-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 flat-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 stiffness on elongation (compressing), bending in a vertical plane, 1 bending in a horizontal plane, and torsion, respectively; α± , α—stiffness 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. Differential equations governing the equilibrium state of a nonhomogeneous shell with respect to displacements are obtained by comparison to zero of the coefficients 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 defined by comparison to zero of the coefficients 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 satisfied 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 coefficients 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 different stiffness: 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 deflection 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 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 {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 coefficients Cα β (4.11) are computed. Then, substituting the partial derivatives by difference 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, 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(x1 x2 ) 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 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 defined 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 flow 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 artificial intensity of the flow 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 deflection, w wx2 = 0 or M22 = 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 x3 = h/2 for the loads corresponding to deflection 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 deflection (see text for more details). Figure 4.3. Influence 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 Deflection 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 effective and also for solving the complicated PDEs like nonlinear equations governing elastic-plastic problems of shells with finite deflections. Here, on each loading step, one has to solve many times a large SLAE, which, without any doubt, exerts significant 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, Scientific 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, Difference 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 flexible 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

Loading next page...
/lp/hindawi-publishing-corporation/on-the-economical-solution-method-for-a-system-of-linear-algebraic-qgk5kpOazQ

You're reading a free preview. Subscribe to read the entire article.

And millions more from thousands of peer-reviewed journals, for just $40/month

To be the best researcher, you need access to the best research

  • With DeepDyve, you can stop worrying about how much articles cost, or if it's too much hassle to order — it's all at your fingertips. Your research is important and deserves the top content.
  • Read from thousands of the leading scholarly journals from Springer, Elsevier, Nature, IEEE, Wiley-Blackwell and more.
  • All the latest content is available, no embargo periods.

Stop missing out on the latest updates in your field

  • We’ll send you automatic email updates on the keywords and journals you tell us are most important to you.
  • There is a lot of content out there, so we help you sift through it and stay organized.