24.2: Gaussian Elimination and Sparsity

If we were to consider the most obvious tactic - find \(K^{-1}\) and then compute \(K^{-1} f\) - the result would be disastrous: days of computing (if the calculation even completed). And indeed even if we were to apply a method of choice - Gaussian elimination (or LU decomposition) - without any regard to the actual structure of the matrix \(K\), the result would still be disastrous. Modern computational solution strategies must and do take advantage of a key attribute of \(K\) - sparseness. \({ }^{2}\) In short, there is no reason to perform operations which involve zero operands and will yield zero for a result. In MechE systems sparseness is not an exceptional attribute but rather, and very fortunately, the rule: the force in a (Hookean) spring is just determined by the deformations of nearest neighbors, not by distant neighbors; similar arguments apply in heat transfer, fluid mechanics, and acoustics. (Note this is not to say that the equilibrium displacement at one node does not depend on the applied forces at the other nodes; quite to the contrary, a force applied at one node will yield nonzero displacements at all the nodes of the system. We explore this nuance more deeply when we explain why formation of the inverse of a (sparse) matrix is a very poor idea.)

We present in Figure \(24.2\) the structure of the matrix \(K\) : the dots indicate nonzero entries in the matrix. We observe that most of the matrix elements are in fact zero. Indeed, of the \(3,603,600,900\) entries of \(K, 3,601,164,194\) entries are zero; put differently, there are only \(2,436,706\)

\({ }^{2}\) Note in this unit we shall consider only direct solution methods; equally important, but more advanced in terms of formulation, analysis, and robust implementation (at least if we consider the more efficient varieties), are iterative solution methods. In actual practice, most state-of-the-art solvers include some combination of direct and iterative aspects. And in all cases, sparsity plays a key role.

Figure 24.2: Structure of stiffness matrix \(K\).

nonzero entries of \(K\) - only \(0.06 \%\) of the entries of \(K\) are nonzero. If we exploit these zeros, both in our numerical approach and in the implementation/programming, we can now solve our system in reasonable time: about 230 seconds on a Mac laptop (performing a particular version of sparse Gaussian elimination based on Cholesky factorization). However, we can do still better.

In particular, although some operations "preserve" all sparsity, other operations - in particular, Gaussian elimination - result in "fill-in": zero entries become nonzero entries which thus must be included in subsequent calculations. The extent to which fill-in occurs depends on the way in which we order the equations and unknowns (which in turn determines the structure of the matrix). There is no unique way that we must choose to order the unknowns and equations: a particular node say near the elbow of the robot arm could be node (column) " 1 " - or node (column) " 2,345 "; similarly, the equilibrium equation for this node could be row "2" - or row " 58,901 " \(.3\) We can thus strive to find a best ordering which minimizes fill-in and thus maximally exploits sparsity. In fact, this optimization problem is very difficult, but there are efficient heuristic procedures which yield very good results. The application of such a heuristic to our matrix \(K\) yields the new (that is, reordered) matrix \(K^{\prime}\) shown in Figure \(24.3\). If we now reapply our sparse Cholesky approach the computational time is now very modest - only 7 seconds. Hence proper choice of algorithm and an appropriate implementation of the algorithm can reduce the computational effort from days to several seconds.