김경종
24 KiB
Hence, we see that the measurement was not correctly performed for the force in clamp 3 when U_{3} = 1 and U_{1} = U_{4} = 0 .
So far we have assumed that no loads are applied to the structure since the operations on the stiffness matrix are independent of the operations on the load vector. When the load vector is not a null vector, the elimination of degrees of freedom proceeds as described above, but in this case the equation used to eliminate a displacement variable from the remaining equations involves a load term. In the elimination the effect of this load is carried over to the remaining degrees of freedom. Therefore, in summary, the physical process of the Gauss elimination procedure is that n stiffness matrices of order n, n - 1, \ldots, 1 , corresponding to the n - i last degrees of freedom, i = 0, 1, 2, \ldots, n - 1 , respectively, of the same physical system are established. In addition, the appropriate load vectors corresponding to the n stiffness matrices are calculated. These load vectors are such that the corresponding displacements at the unreleased degrees of freedom are the displacements calculated for the system when described by all n degrees of freedom. The unknown displacements are then obtained by considering in succession the systems with only one, two, \ldots , degrees of freedom (these correspond to the last, two last, \ldots , of the original number of degrees of freedom).
We can now explain why, because of physical reasons, all diagonal elements in the Gauss elimination procedure must remain positive. This follows because the final ith diagonal element is the stiffness at degree of freedom i when the first i - 1 degrees of freedom of the system have been released, and this stiffness should be positive. If a zero (or negative) diagonal element occurs in the Gauss elimination, the structure is not stable. An example of such a case is shown in Fig. 8.7, where after the release of degrees of freedom U_{1} , U_{2} , and U_{3} , the last diagonal element is zero.
text_image
Hinge U₃ U₄ U₁ U₂
Flexural rigidity EI
Figure 8.7 Example of an unstable structure
So far we have considered that the Gauss elimination proceeds in succession from the first to the (n - 1) st degree of freedom. However, we may in the same way perform the elimination backward (i.e., from the last to the second degree of freedom), or we may choose any desirable order, as we already indicated in the discussion of the physical laboratory experiment shown in Fig. 8.2.
EXAMPLE 8.2: Obtain the solution to the equilibrium equations of the beam in Fig. 8.1 by eliminating the displacement variables in the order U_{3} , U_{2} , U_{4} .
We may either write down the individual equations during the elimination process or directly perform Gauss elimination in the prescribed order. In the latter case we obtain,
eliminating U_{3} ,
\left[ \begin{array}{c c c c} \frac {2 9}{6} & - \frac {1 0}{3} & 0 & \frac {7}{3} \\ - \frac {1 0}{3} & \frac {1 0}{3} & 0 & - \frac {5}{3} \\ 1 & - 4 & 6 & - 4 \\ \frac {2}{3} & - \frac {5}{3} & 0 & \frac {7}{3} \end{array} \right] \left[ \begin{array}{l} U _ {1} \\ U _ {2} \\ U _ {3} \\ U _ {4} \end{array} \right] = \left[ \begin{array}{l} 0 \\ 1 \\ 0 \\ 0 \end{array} \right]
Next, we eliminate U_{2} and obtain
\left[ \begin{array}{c c c c} \frac {3}{2} & 0 & 0 & - 1 \\ - \frac {1 0}{3} & \frac {1 0}{3} & 0 & - \frac {5}{3} \\ 1 & - 4 & 6 & - 4 \\ - 1 & 0 & 0 & \frac {3}{2} \end{array} \right] \left[ \begin{array}{c} U _ {1} \\ U _ {2} \\ U _ {3} \\ U _ {4} \end{array} \right] = \left[ \begin{array}{c} 1 \\ 1 \\ 0 \\ \frac {1}{2} \end{array} \right]
and finally we eliminate U_{4} ,
\left[ \begin{array}{c c c c} \frac {5}{6} & 0 & 0 & 0 \\ - \frac {1 0}{3} & \frac {1 0}{3} & 0 & - \frac {5}{3} \\ 1 & - 4 & 6 & - 4 \\ - 1 & 0 & 0 & \frac {3}{2} \end{array} \right] \left[ \begin{array}{l} U _ {1} \\ U _ {2} \\ U _ {3} \\ U _ {4} \end{array} \right] = \left[ \begin{array}{l} \frac {4}{3} \\ 1 \\ 0 \\ \frac {1}{2} \end{array} \right]
The solution for the displacements is now obtained as follows:
\begin{array}{l} U _ {1} = \frac {8}{5} \\ U _ {4} = \frac {\frac {1}{2} - (- 1) \frac {8}{5}}{\frac {3}{2}} = \frac {7}{5} \\ U _ {2} = \frac {1 - (- \frac {1 0}{3}) \frac {8}{5} - (- \frac {5}{3}) \frac {7}{5}}{\frac {1 0}{3}} = \frac {1 3}{5} \\ U _ {3} = \frac {0 - (1) \frac {8}{5} - (- 4) \frac {7}{5} - (- 4) \frac {1 3}{5}}{6} = \frac {1 2}{5} \\ \end{array}
which is the solution obtained earlier.
8.2.2 The LDL ^{T} Solution
We have seen in the preceding section that the basic procedure of the Gauss elimination solution is to reduce the equations to correspond to an upper triangular coefficient matrix from which the unknown displacements U can be calculated by a back-substitution. We now want to formalize the solution procedure using appropriate matrix operations. An additional important purpose of the discussion is to introduce a notation that can be used throughout the following presentations. The actual computer implementation is given in the next section.
Considering the operations performed in the Gauss elimination solution presented in the preceding section, the reduction of the stiffness matrix K to upper triangular form can be written
\mathbf {L} _ {n - 1} ^ {- 1} \dots \mathbf {L} _ {2} ^ {- 1} \mathbf {L} _ {1} ^ {- 1} \mathbf {K} = \mathbf {S} \tag {8.10}
where S is the final upper triangular matrix and
\mathbf {L} _ {i} ^ {- 1} = \left[ \begin{array}{c c c c c c c c} 1 & & & & \text {elements} & & \\ & \cdot & & & \text {not shown} & & \\ & & \cdot & & & \text {are zeros} ^ {1} & & \\ & & & 1 & & & \\ & & & - l _ {i + 1, i} & . & & \\ & & & - l _ {i + 2, i} & & & \\ & & & \cdot & & & \\ & & & \cdot & & & \\ & & & - l _ {n i} & & & 1 \end{array} \right]; \quad l _ {i + j, i} = \frac {k _ {i + j , i} ^ {(i)}}{k _ {i i} ^ {(i)}} \tag {8.11}
The elements l_{i+j,i} are the Gauss multiplying factors, and the right superscript (i) indicates that an element of the matrix L_{i-1}^{-1} \ldots L_{2}^{-1}L_{1}^{-1}K is used.
We now note that L_{i} is obtained by simply reversing the signs of the off-diagonal elements in L_{i}^{-1} . Therefore, we obtain
\mathbf {K} = \mathbf {L} _ {1} \mathbf {L} _ {2} \dots \mathbf {L} _ {n - 1} \mathbf {S} \tag {8.12}
where
\mathbf {L} _ {i} = \left[ \begin{array}{c c c c c c c c c} 1 & & & & & & & \\ & \cdot & & & & & & \\ & & \cdot & & & & & \\ & & & \cdot & & & & \\ & & & & 1 & & & \\ & & & & l _ {i + 1, i} & & & & \\ & & & & l _ {i + 2, i} & & & & \\ & & & & \cdot & & \cdot & & \\ & & & & \cdot & & & \cdot & \\ & & & & l _ {n i} & & & & 1 \end{array} \right] \tag {8.13}
Hence, we can write
\mathbf {K} = \mathbf {L S} \tag {8.14}
where \mathbf{L} = \mathbf{L}_1\mathbf{L}_2\ldots \mathbf{L}_{n - 1} ; i.e., \mathbf{L} is a lower unit triangular matrix,
\mathbf {L} = \left[ \begin{array}{c c c c c c c c c} 1 & & & & & & & & \\ l _ {2 1} & 1 & & & & & & & \\ l _ {3 1} & l _ {3 2} & 1 & & & & & & \\ l _ {4 1} & l _ {4 2} & & 1 & & & & & \\ \cdot & \cdot & & & \cdot & & & & \\ \cdot & \cdot & & & & \cdot & & & \\ \cdot & \cdot & & & & & \cdot & & \\ & & & & & & & 1 & \\ l _ {n 1} & \cdot & & \cdot & & \cdot & & l _ {n, n - 1} & 1 \end{array} \right] \tag {8.15}
Since S is an upper triangular matrix and the diagonal elements are the pivots in the Gauss elimination, we can write S = D\bar{S} , where D is a diagonal matrix storing the diagonal
^{1} Throughout this book, elements not shown in matrices are zero.
elements of \mathbf{S} ; i.e., d_{ii} = s_{ii} . Substituting for \mathbf{S} into (8.14) and noting that \mathbf{K} is symmetric and the decomposition is unique, we obtain \tilde{\mathbf{S}} = \mathbf{L}^T , and hence,
\mathbf {K} = \mathbf {L D L} ^ {T} \tag {8.16}
It is this LDL^{T} decomposition of K that can be used effectively to obtain the solution of the equations in (8.1) in the following two steps:
\mathbf {L V} = \mathbf {R} \tag {8.17}
\mathbf {D L} ^ {T} \mathbf {U} = \mathbf {V} \tag {8.18}
where in (8.17) the load vector R is reduced to obtain V,
\mathbf {V} = \mathbf {L} _ {n - 1} ^ {- 1} \dots \mathbf {L} _ {2} ^ {- 1} \mathbf {L} _ {1} ^ {- 1} \mathbf {R} \tag {8.19}
and in (8.18) the solution \mathbf{U} is obtained by a back-substitution,
\mathbf {L} ^ {T} \mathbf {U} = \mathbf {D} ^ {- 1} \mathbf {V} \tag {8.20}
In the implementation the vector V is frequently calculated at the same time as the matrices L_{i}^{-1} are established. This was done in the example solution of the simply supported beam in Section 8.2.1.
It should be noted that in practice the matrix multiplications to obtain L in (8.15) and V in (8.19) are not formally carried out, but that L and V are established by directly modifying K and R. This is discussed further in the next section, in which the computer implementation of the solution procedure is presented. However, before proceeding, consider the example in Section 8.2.1 for the derivation of the matrices defined above.
EXAMPLE 8.3: Establish the matrices L_{i}^{-1} , i = 1, 2, 3, S, L, and D and the vector V corresponding to the stiffness matrix and the load vector of the simply supported beam treated in Section 8.2.1.
Using the information given in Section 8.2.1, we can directly write down the required matrices:
\begin{array}{l} \mathbf {L} _ {1} ^ {- 1} = \left[ \begin{array}{c c c c} 1 & & & \\ \frac {4}{5} & 1 & & \\ - \frac {1}{5} & 0 & 1 & \\ 0 & 0 & 0 & 1 \end{array} \right]; \\ \mathbf {L} _ {2} ^ {- 1} = \left[ \begin{array}{c c c c} 1 & & & \\ 0 & 1 & & \\ 0 & \frac {8}{7} & 1 & \\ 0 & - \frac {5}{1 4} & 0 & 1 \end{array} \right] \\ \mathbf {L} _ {3} ^ {- 1} = \left[ \begin{array}{c c c c} 1 & & & \\ 0 & 1 & & \\ 0 & 0 & 1 & \\ 0 & 0 & \frac {4}{3} & 1 \end{array} \right]; \\ \mathbf {S} = \left[ \begin{array}{c c c c} 5 & - 4 & 1 & 0 \\ & \frac {1 4}{5} & - \frac {1 6}{5} & 1 \\ & & \frac {1 5}{7} & - \frac {2 0}{7} \\ & & & \frac {5}{6} \end{array} \right] \\ \end{array}
where the matrix L_{i}^{-1} stores in the ith column the multipliers that were used in the elimination of the ith equation and the matrix S is the upper triangular matrix obtained in (8.5). The matrix D is a diagonal matrix with the pivot elements on its diagonal. In this case,
\mathbf {D} = \left[ \begin{array}{c c c c} 5 & & & \\ & \frac {1 4}{5} & & \\ & & \frac {1 5}{7} & \\ & & & \frac {5}{6} \end{array} \right]
To obtain L we use (8.15); hence,
\mathbf {L} = \left[ \begin{array}{c c c c} 1 & & & \\ - \frac {4}{5} & 1 & & \\ \frac {1}{5} & - \frac {8}{7} & 1 & \\ 0 & \frac {5}{1 4} & - \frac {4}{3} & 1 \end{array} \right]
and we can check that \mathbf{S} = \mathbf{DL}^T .
The vector V was obtained in (8.5):
\mathbf {V} = \left[ \begin{array}{l} 0 \\ 1 \\ \frac {8}{7} \\ \frac {7}{6} \end{array} \right]
8.2.3 Computer Implementation of Gauss Elimination—The Active Column Solution
The aim in the computer implementation of the Gauss solution procedure is to use a small solution time. In addition, the high-speed storage requirements should be small to avoid the use of backup storage. However, for large systems it will nevertheless be necessary to use backup storage, and for this reason it should also be possible to modify the solution algorithm for effective out-of-core solution.
An advantage of finite element analysis is that the stiffness matrix of the element assemblage is not only symmetric and positive definite but also banded; i.e., k_{ij} = 0 for j > i + m_{K} , where m_{K} is the half-bandwidth of the system (see Fig. 2.1). The fact that in finite element analysis all nonzero elements are clustered around the diagonal of the system matrices greatly reduces the total number of operations and the high-speed storage required in the equation solution. However, this property depends on the nodal point numbering in the finite element mesh, and the analyst must take care to obtain an effective nodal point numbering (see Chapter 12).
Assume that for a given finite element assemblage a specific nodal point numbering has been determined and the corresponding column heights and the stiffness matrix K have been calculated (see Section 12.2.3 for details). The LDL^{T} decomposition of K can be obtained effectively by considering each column in turn; i.e., although the Gauss elimination is carried out by rows, the final elements of D and L are calculated by columns. Using d_{11} = k_{11} , the algorithm for the calculation of the elements l_{ij} and d_{jj} in the jth column is, for j = 2, \ldots, n ,
\left. \begin{array}{l} g _ {m _ {j}, j} = k _ {m _ {j}, j} \\ g _ {i j} = k _ {i j} - \sum_ {r = m _ {m}} ^ {i - 1} l _ {r i} g _ {r j} \quad i = m _ {j} + 1, \dots , j - 1 \end{array} \right\} \tag {8.21}
where m_{j} is the row number of the first nonzero element in column j and m_{m} = \max\{m_{i}, m_{j}\} (see Fig. 12.2). The variables m_{i}, i = 1, \ldots, n , define the skyline of the matrix; also, the values i - m_{i} are the column heights and the maximum column height is the half-bandwidth m_{K} . The elements g_{ij} in (8.21) are defined only as intermediate quantities, and
the calculation is completed using
l _ {i j} = \frac {g _ {i j}}{d _ {i i}} \quad i = m _ {j}, \dots , j - 1 \tag {8.22}
d _ {i j} = k _ {i j} - \sum_ {r = m _ {j}} ^ {j - 1} l _ {r j} g _ {r j} \tag {8.23}
It should be noted that the summations in (8.21) and (8.23) do not involve multiplications with zero elements outside the skyline of the matrix and that the l_{ij} are elements of the matrix L^{T} rather than of L. We refer to the solution algorithm given in (8.21) to (8.23) [actually with (8.24) and (8.25)] as the active column solution or the skyline (or column) reduction method.
Considering the storage arrangements in the reduction, the element l_{ij} when calculated for use in (8.23) immediately replaces g_{ij} , and d_{jj} replaces k_{jj} . Therefore, at the end of the reduction we have elements d_{jj} in the storage locations previously used by k_{jj} , and l_{rj} is stored in the locations of k_{rj} , j > r.
In order to get familiar with the solution algorithm, we consider the following examples.
EXAMPLE 8.4: Use the solution algorithm given in (8.21) to (8.23) to calculate the triangular factors D and L^{T} of the stiffness matrix of the beam considered in Example 8.3.
The initial elements considered are, when written in their respective matrix locations,
\left[ \begin{array}{c c c c} 5 & - 4 & 1 & \\ & 6 & - 4 & 1 \\ & & 6 & - 4 \\ & & & 5 \end{array} \right]
with m_1 = 1 , m_2 = 1 , m_3 = 1 , and m_4 = 2 . Using (8.21) to (8.23), we obtain, for j = 2 ,
d _ {1 1} = k _ {1 1} = 5
g _ {1 2} = k _ {1 2} = - 4
l _ {1 2} = \frac {g _ {1 2}}{d _ {1 1}} = - \frac {4}{5}
d _ {2 2} = k _ {2 2} - l _ {1 2} g _ {1 2} = 6 - (- 4) \left(- \frac {4}{5}\right) = \frac {1 4}{5}
and thus the resulting matrix elements are now, using a dotted line to separate the reduced from the unreduced columns,
\left[ \begin{array}{c c c c} 5 & - \frac {4}{5} & 1 & \\ & \frac {1 4}{5} & - 4 & 1 \\ & & 6 & - 4 \\ & & & 5 \end{array} \right]
Next, we obtain, for j = 3,
g _ {1 3} = k _ {1 3} = 1
g _ {2 3} = k _ {2 3} - l _ {1 2} g _ {1 3} = - 4 - (- \frac {4}{5}) (1) = - \frac {1 6}{5}
l _ {1 3} = \frac {g _ {1 3}}{d _ {1 1}} = \frac {1}{5}
l _ {2 3} = \frac {g _ {2 3}}{d _ {2 2}} = \frac {- \frac {1 6}{5}}{\frac {1 4}{5}} = - \frac {8}{7}
d _ {3 3} = k _ {3 3} - l _ {1 3} g _ {1 3} - l _ {2 3} g _ {2 3} = 6 - (\frac {1}{5}) (1) - (- \frac {8}{7}) (- \frac {1 6}{5}) = \frac {1 5}{7}
and the resulting matrix elements are
\left[ \begin{array}{c c c c} 5 & - \frac {4}{5} & \frac {1}{5} & \\ & \frac {1 4}{5} & - \frac {8}{7} & 1 \\ & & \frac {1 5}{7} & - 4 \\ & & & 5 \end{array} \right]
Finally, we have, for j = 4,
g _ {2 4} = k _ {2 4} = 1
g _ {3 4} = k _ {3 4} - l _ {2 3} g _ {2 4} = - 4 - (- \frac {8}{7}) (1) = - \frac {2 0}{7}
l _ {2 4} = \frac {g _ {2 4}}{d _ {2 2}} = \frac {1}{\frac {1 4}{5}} = \frac {5}{1 4}
l _ {3 4} = \frac {g _ {3 4}}{d _ {3 3}} = \frac {- \frac {2 0}{7}}{\frac {1 5}{7}} = - \frac {4}{3}
d _ {4 4} = k _ {4 4} - l _ {2 4} g _ {2 4} - l _ {3 4} g _ {3 4} = 5 - (\frac {5}{1 4}) (1) - (- \frac {4}{3}) (- \frac {2 0}{7}) = \frac {5}{6}
and the final elements stored are
\left[ \begin{array}{c c c c} 5 & - \frac {4}{5} & \frac {1}{5} & \\ & \frac {1 4}{5} & - \frac {8}{7} & \frac {5}{1 4} \\ & & \frac {1 5}{7} & - \frac {4}{3} \\ & & & \frac {5}{6} \end{array} \right]
We should note that the elements of \mathbf{D} are stored on the diagonal and the elements l_{ij} have replaced the elements k_{ij}, j > i .
Although the details of the solution procedure have already been demonstrated in Example 8.4, the importance of using the column reduction method could not be shown, since the skyline coincides with the band. The effectiveness of the skyline reduction scheme is more apparent in the factorization of the following matrix.
EXAMPLE 8.5: Use the solution algorithm given in (8.21) to (8.23) to evaluate the triangular factors D and L^{T} of the stiffness matrix K, where
\mathbf {K} = \left[ \begin{array}{c c c c c} 2 & - 2 & & & - 1 \\ & 3 & - 2 & & 0 \\ & & 5 & - 3 & 0 \\ & \text { Symmetric } & & 1 0 & 4 \\ & & & & 1 0 \end{array} \right]
For this matrix we have m_1 = 1 , m_2 = 1 , m_3 = 2 , m_4 = 3 , and m_5 = 1 .
The algorithm gives, in this case with d_{11} = 2 , for j = 2 ,
g _ {1 2} = k _ {1 2} = - 2
l _ {1 2} = \frac {g _ {1 2}}{d _ {1 1}} = \frac {- 2}{2} = - 1
d _ {2 2} = k _ {2 2} - l _ {1 2} g _ {1 2} = 3 - (- 1) (- 2) = 1
and thus the resulting matrix elements are
\left[ \begin{array}{c c c c c c} 2 & - 1 & & & & - 1 \\ & 1 & - 2 & & & 0 \\ & & 5 & - 3 & & 0 \\ & & & 1 0 & & 4 \\ & & & & & 1 0 \end{array} \right]
For j = 3 , g_{23} = k_{23} = -2
l _ {2 3} = \frac {g _ {2 3}}{d _ {2 2}} = \frac {- 2}{1} = - 2
d _ {3 3} = k _ {3 3} - l _ {2 3} g _ {2 3} = 5 - (- 2) (- 2) = 1
and the coefficient array is now
\left[ \begin{array}{c c c c c} 2 & - 1 & & & - 1 \\ & 1 & - 2 & & 0 \\ & & 1 & - 3 & 0 \\ & & & 1 0 & 4 \\ & & & & 1 0 \end{array} \right]
For j = 4 , g_{34} = k_{34} = -3
l _ {3 4} = \frac {g _ {3 4}}{d _ {3 3}} = \frac {- 3}{1} = - 3
d _ {4 4} = k _ {4 4} - l _ {3 4} g _ {3 4} = 1 0 - (- 3) (- 3) = 1
and the resulting matrix elements are
\left[ \begin{array}{c c c c c c} 2 & - 1 & & & & - 1 \\ & 1 & - 2 & & & 0 \\ & & 1 & - 3 & & 0 \\ & & & 1 & & 4 \\ & & & & & 1 0 \end{array} \right]
Finally we have, for j = 5,
g _ {1 5} = k _ {1 5} = - 1
g _ {2 5} = k _ {2 5} - l _ {1 2} g _ {1 5} = 0 - (- 1) (- 1) = - 1
g _ {3 5} = k _ {3 5} - l _ {2 3} g _ {2 5} = 0 - (- 2) (- 1) = - 2
g _ {4 5} = k _ {4 5} - l _ {3 4} g _ {3 5} = + 4 - (- 3) (- 2) = - 2
\begin{array}{l} l _ {1 5} = \frac {g _ {1 5}}{d _ {1 1}} = \frac {- 1}{2} = - \frac {1}{2} \\ l _ {2 5} = \frac {g _ {2 5}}{d _ {2 2}} = \frac {- 1}{1} = - 1 \\ l _ {3 5} = \frac {g _ {3 5}}{d _ {3 3}} = \frac {- 2}{1} = - 2 \\ l _ {4 5} = \frac {g _ {4 5}}{d _ {4 4}} = \frac {- 2}{1} = - 2 \\ d _ {5 5} = k _ {5 5} - l _ {1 5} g _ {1 5} - l _ {2 5} g _ {2 5} - l _ {3 5} g _ {3 5} - l _ {4 5} g _ {4 5} \\ = 1 0 - \left(- \frac {1}{2}\right) (- 1) - (- 1) (- 1) - (- 2) (- 2) - (- 2) (- 2) = \frac {1}{2} \\ \end{array}
\begin{array}{l} d _ {5 5} = k _ {5 5} - l _ {1 5} g _ {1 5} - l _ {2 5} g _ {2 5} - l _ {3 5} g _ {3 5} - l _ {4 5} g _ {4 5} \\ = 1 0 - \left(- \frac {1}{2}\right) (- 1) - (- 1) (- 1) - (- 2) (- 2) - (- 2) (- 2) = \frac {1}{2} \\ \end{array}
and the final matrix elements are
\left[ \begin{array}{c c c c c} 2 & - 1 & & & - \frac {1}{2} \\ & 1 & - 2 & & - 1 \\ & & 1 & - 3 & - 2 \\ & & & 1 & - 2 \\ & & & & \frac {1}{2} \end{array} \right]
As in Example 8.4, we have the elements of \mathbf{D} and \mathbf{L}^T replacing the elements k_{ii} and k_{ij}, j > i of the original matrix \mathbf{K} , respectively.
In the preceding discussion we considered only the decomposition of the stiffness matrix K, which constitutes the main part of the equation solution. Once the L, D factors of K have been obtained, the solution for U is calculated using (8.19) and (8.20), where it may be noted that the reduction of R in (8.19) can be performed at the same time as the stiffness matrix K is decomposed or may be carried out separately afterward. The equation to be used is similar to (8.23); i.e., we have V_{1} = R_{1} and calculate for i = 2, \ldots, n ,
V _ {i} = R _ {i} - \sum_ {r = m _ {i}} ^ {i - 1} l _ {r i} V _ {r} \tag {8.24}
where R_{i} and V_{i} are the i th elements of \mathbf{R} and \mathbf{V} . Considering the storage arrangements, the element V_{i} replaces R_{i} .
The back-substitution in (8.20) is performed by evaluating successively U_{n} , U_{n-1}, \ldots, U_{1} . This is achieved by first calculating \overline{V} , where \overline{V} = D^{-1}V . Then using \overline{\mathbf{V}}^{(n)} = \overline{\mathbf{V}} , we have U_{n} = \overline{V}_{n}^{(n)} and we calculate for i = n, \ldots, 2 ,
\left. \begin{array}{c} \overline {{V}} _ {r} ^ {(i - 1)} = \overline {{V}} _ {r} ^ {(i)} - l _ {r i} U _ {i}; \quad r = m _ {i}, \dots , i - 1 \\ U _ {i - 1} = \overline {{V}} _ {i - 1} ^ {(i - 1)} \end{array} \right\} \tag {8.25}
where the superscript (i - 1) indicates that the element is calculated in the evaluation of U_{i-1} . It should be noted that \overline{V}_{k}^{(j)} for all j is stored in the storage location of V_{k} , i.e., the original storage location of R_{k} .
EXAMPLE 8.6: Use the algorithm given in (8.24) and (8.25) to calculate the solution to the problem \mathbf{KU} = \mathbf{R} , when \mathbf{K} is the stiffness matrix considered in Example 8.5 and
\mathbf {R} = \left[ \begin{array}{l} 0 \\ 1 \\ 0 \\ 0 \\ 0 \end{array} \right]
In the solution we employ the D, L^{T} factors of K calculated in Example 8.5. Using (8.24) for the forward reduction, we obtain
\begin{array}{l} V _ {1} = R _ {1} = 0 \\ V _ {2} = R _ {2} - l _ {1 2} V _ {1} = 1 - 0 = 1 \\ V _ {3} = R _ {3} - l _ {2 3} V _ {2} = 0 - (- 2) (1) = 2 \\ V _ {4} = R _ {4} - l _ {3 4} V _ {3} = 0 - (- 3) (2) = 6 \\ V _ {5} = R _ {5} - l _ {1 5} V _ {1} - l _ {2 5} V _ {2} - l _ {3 5} V _ {3} - l _ {4 5} V _ {4} \\ = 0 - 0 - (- 1) (1) - (- 2) (2) - (- 2) (6) = 1 7 \\ \end{array}
\begin{array}{l} V _ {5} = R _ {5} - l _ {1 5} V _ {1} - l _ {2 5} V _ {2} - l _ {3 5} V _ {3} - l _ {4 5} V _ {4} \\ = 0 - 0 - (- 1) (1) - (- 2) (2) - (- 2) (6) = 1 7 \\ \end{array}
Immediately after calculation of V_{i} , the element replaces R_{i} . Thus, we now have in the vector that initially stored the loads,
\mathbf {V} = \left[ \begin{array}{l} 0 \\ 1 \\ 2 \\ 6 \\ 1 7 \end{array} \right]
The first step in the back-substitution is to evaluate \overline{V} , where \overline{V} = D^{-1}V . Here we obtain
\overline {{{\mathbf {V}}}} = \left[ \begin{array}{l} 0 \\ 1 \\ 2 \\ 6 \\ 3 4 \end{array} \right] \tag {a}
and thus, U_{5} = \overline{V}_{5} = 34
Now we use (8.25) with \overline{\mathbf{V}}^{(5)} = \overline{\mathbf{V}} of (a). Hence we obtain, for i = 5 ,
\begin{array}{l} \overline {{{V}}} _ {1} ^ {(4)} = \overline {{{V}}} _ {1} ^ {(5)} - l _ {1 5} U _ {5} = 0 - (- \frac {1}{2}) (3 4) = 1 7 \\ \overline {{{V}}} _ {2} ^ {(4)} = \overline {{{V}}} _ {2} ^ {(5)} - l _ {2 5} U _ {5} = 1 - (- 1) (3 4) = 3 5 \\ \overline {{{V}}} _ {3} ^ {(4)} = \overline {{{V}}} _ {3} ^ {(5)} - l _ {3 5} U _ {5} = 2 - (- 2) (3 4) = 7 0 \\ \overline {{{V}}} _ {4} ^ {(4)} = \overline {{{V}}} _ {4} ^ {(5)} - l _ {4 5} U _ {5} = 6 - (- 2) (3 4) = 7 4 \\ \end{array}
and U_{4} = \overline{V}_{4}^{(4)} = 74
For i = 4 , \overline{V}_3^{(3)} = \overline{V}_3^{(4)} - l_{34}U_4 = 70 - (-3)(74) = 292