The matrix that replaced the hand methods

Almost every analytical method in this collection was invented to be done by hand. Moment distribution passes moments round a frame until they stop moving; virtual work computes one deflection at a time from a unit load; the area-moment theorems read rotations off a diagram. All of them are exact, all of them are ingenious, and all of them stop scaling somewhere in the low tens of members.

What replaced them adds no physics whatever. It writes the same equilibrium and the same compatibility as one matrix equation,

$$\mathbf{K}\,\mathbf{d} = \mathbf{f}$$

assembled by scattering each member’s own six-by-six stiffness into the rows and columns of the freedoms it touches, and then solves it. The whole of the invention is the bookkeeping, and it is worth looking at what that bookkeeping produces, because three properties of the resulting matrix are structural statements in their own right.

Which free body produced the number #

The member is the free body, and it is a small one. Take a plane-frame element with three freedoms at each end — two translations and a rotation — and ask what forces are needed at its ends to produce a unit value of each freedom with the other five held at zero. Six questions, six answers each, and the six-by-six matrix of them is the element stiffness.

Every entry of it comes from an elementary beam solution. A unit end rotation with everything else fixed needs $4EI/L$ at that end and $2EI/L$ at the other; a unit end translation needs $12EI/L^3$ and $6EI/L^2$; a unit axial extension needs $EA/L$. Nothing there is new — those are the fixed-end moments of the slope-deflection method, written out as a table instead of as equations.

Assembly is two operations and no physics. Rotate the element matrix into the structure’s axes, then add each of its thirty-six entries into the row and column of the global freedom that entry belongs to. Every member does the same; the sum is the structure.

The equilibrium is in the addition — the forces at a node from every member meeting there must sum to whatever is applied — and the compatibility is in the sharing of an index, because two members attached to the same node are being given the same displacement by construction. Those are the two conditions every method in this field enforces, and here they are enforced by an array index rather than by an argument.

The check, which has to come from outside #

A method that assembles its own answer from its own parts can be internally consistent and wrong, so the assembly needs to be held against something it does not contain.

The propped cantilever is the right test because its answer is classical and short: a cantilever of span $L$ under a uniform load $w$, propped at its free end, has a prop reaction of $3wL/8$. That number comes from compatibility — the prop deflection cancels the free cantilever’s — and nothing in a stiffness assembly derives it.

Assembled from sixteen elements and solved, the prop comes out at $3wL/8$ to a part in a billion, which is the residual of the elimination and not of the physics. The method is exact for a structure made of members whose behaviour is the element matrix’s, and the only approximation anywhere is the arithmetic’s.

Singular until the supports go in #

The first structural property of the matrix is one nobody sets out to produce.

Assemble a frame and do not restrain it. The resulting $\mathbf{K}$ is singular, and its null space has exactly three dimensions: the two rigid-body translations and the rigid-body rotation that a plane structure can undergo without straining anything. A vector in that null space is a displacement pattern that costs no energy, so $\mathbf{K}$ maps it to zero.

That is the same statement a truss that is a mechanism makes, arrived at from the other end. There, the equilibrium matrix has a null vector and it is a mechanism of the structure; here, the stiffness matrix has three of them and they are mechanisms of the structure with respect to the ground. An unsupported structure is a mechanism, and the linear algebra says so before any structural reasoning does.

The practical value of that is diagnostic. A frame with a member accidentally disconnected, or a joint released in both rotation and translation, has a fourth null vector — and the solver reports a singular matrix, which is a true and unhelpful message. Computing the null space instead returns the displacement pattern that costs nothing, and that pattern points straight at the missing restraint. Nearly every finite-element package has that facility and it is under-used.

Counting restraints against equations is the hand version of the same question, and the comparison is instructive about what each instrument can do. The count returns a single integer and knows nothing about where the deficiency is; the null space returns a displacement pattern and points at it. Both detect the same defect, and only one of them can be drawn on the frame.

Sparse, banded, and the numbering decides which #

The second property is the one that made the method practical rather than merely correct.

A member reaches only its own two nodes, so it can contribute nothing to any row belonging to a node it does not touch. Every such entry is zero exactly rather than nearly, and for a frame of any size the great majority of the matrix is zeros: 16 per cent non-zero for the small frame above, 5 per cent for a twenty-storey one, and falling.

More usefully, the non-zeros sit in a band about the diagonal, whose width is the largest difference between the global freedom numbers at the two ends of any member. And that number is a property of the labelling rather than of the frame — the same eight-storey matrix drawn above has a second version of itself with the same entries in different places and a band a third as wide.

Take a one-bay, twenty-storey frame. Number the nodes column by column — all the way up one column, then all the way up the other — and every beam connects nodes twenty-one apart, giving a bandwidth of 65. Number them storey by storey and every member connects nodes at most two apart, giving 8.

A banded elimination costs about $n b^2$ operations against a dense $n^3/3$. So the two numberings of the same frame differ by a factor of sixty-six in the work of solving it, for answers identical to the last bit. That is why every solver in existence renumbers before it factorises, using a bandwidth-reducing permutation, and why nothing about the physics changes when it does.

There is a general point in that worth extracting. The structure of a matrix is not a property of the problem, it is a property of the ordering — and a great many computational costs in this subject are decided by choices that are invisible in the engineering.

Symmetry, and what it is a statement about #

There is a third property, and it is the one that halves the cost of everything above.

$\mathbf{K}$ is symmetric: the force at freedom $i$ due to a unit displacement at $j$ equals the force at $j$ due to a unit displacement at $i$. That is not a convenience of the algebra, it is Maxwell’s reciprocal theorem restated — a statement about linear elastic bodies proved in 1864, three-quarters of a century before anybody assembled a matrix.

Every solver relies on it. Symmetric storage halves the memory; a Cholesky factorisation instead of a general elimination halves the operations again; and the symmetry is what allows the factorisation to be done without pivoting, which is what makes the band survive the solve. A non-symmetric formulation — which follower forces and some contact problems produce — loses all three at once, and the cost of analysing such a problem is several times that of an ordinary one for reasons that are entirely about the theorem no longer holding.

The physical content of the symmetry is worth restating because it is easy to lose. It says that a structure has a strain energy — a scalar function of its displacements — and that the stiffness matrix is that function’s second derivative. Symmetry follows because mixed partial derivatives commute. Anything that breaks the symmetry is something that has no potential, which is to say a load that does work round a closed cycle, and such loads are the ones that produce instabilities of a kind no eigenvalue of a symmetric matrix can find.

Where the reactions come from, and why that check is weak #

The supports enter by deleting things. A restrained freedom is known — usually zero — so its equation is struck out and its column with it, and what is left is a smaller system in the free freedoms alone. That is the operation the previous section’s singularity was waiting for: the null space vanishes the moment three independent restraints are removed.

What happens to the deleted equations is worth following, because it decides how much a common check is worth. They are not discarded. After the free displacements are solved, each struck-out row is re-evaluated — multiplied through by the displacement vector — and the number that comes out is that support’s reaction. The reactions are the residuals of the equations that were removed.

Which means the standard verification, summing the reactions and comparing with the applied load, is not the independent check it appears to be. It re-evaluates rows of the same $\mathbf{K}$ against the same $\mathbf{d}$, so it will close whenever the elimination was arithmetically sound — whatever the element stiffnesses were. A member entered with ten times its true second moment produces a wrong displacement field, wrong member forces, and reactions that sum perfectly. That is exactly the class of error a global equilibrium check cannot see, and this is the mechanism: the check and the answer come out of one matrix.

The fourth bar is the whole of the argument above, drawn. Equilibrium of the whole body is one equation for each of its own three freedoms, and a stiffness error lives entirely in the many equations underneath that — so a model can satisfy every equilibrium check ever devised and still be a model of a different structure.

What would be independent is a free body drawn by hand across the whole structure — total load against total reaction, computed from the loading rather than from the model. It takes a minute, it uses no element matrix at all, and it is the only version of that check with anything outside the solve in it.

What it asserts, and what it does not #

It is worth being precise about the method’s claims, because it is easy to credit it with more than it has.

It asserts nothing about the members. The element stiffness is a beam solution: plane sections, small displacements, linear elasticity, no shear deformation unless it is put in. Every limitation of beam theory is inherited unchanged, and the matrix does not improve any of them.

It asserts nothing about the joints. A shared node means a rigid connection, and a real connection is neither pinned nor rigid. Modelling a semi-rigid joint means putting a rotational spring in, which is a decision somebody has to make and a stiffness somebody has to know.

And it gives displacements. Forces come out afterwards by multiplying the element matrix by its own end displacements, and stresses come out after that by a section calculation the method has not performed. The quantity that decides a design is nearly always a stress at a point that no element in the model contains — which is why the output of an analysis is the input to a design and not a substitute for it.

What was lost, which is not nothing #

Moment distribution was not merely slower. It was legible: every cycle showed a moment moving from one joint to the next, the sequence made the load path visible, and an engineer developed a feel for which joints mattered by watching them converge.

A matrix solve shows none of that. It produces a vector of displacements in one step, and the intermediate states have no physical meaning at all — the partial results of a Gaussian elimination are not partial structures. So the method that scales is also the method that hides its own reasoning, and the compensating discipline is to know what the answer should roughly be before the answer arrives.

There is a related loss worth naming. Hand methods forced a choice of which quantity to compute. Virtual work computes the deflection asked for and nothing else; the flexibility method computes the redundant asked for. A stiffness solve computes everything, whether or not anybody wanted it, and the temptation to read all of it as meaningful is the origin of a good deal of bad engineering. A number the analysis printed is not the same as a number the analysis is about.

The order it arrived in, which is backwards #

The history is worth a paragraph because it explains why the method looks the way it does rather than the way its name suggests.

The flexibility method came first and is the natural hand method: choose redundants, release them, compute the deflections at the released points under the applied load and under unit values of each redundant, and solve a small system for the redundants. It is small — the system is as large as the degree of indeterminacy, which for a hand-analysable frame is a handful — and it needs judgement, because the released structure has to be chosen and a bad choice gives an ill-conditioned system.

The stiffness method is the opposite in every respect. Its system is as large as the number of freedoms, which is far larger than the degree of indeterminacy — thirty-six against nine for the small frame above. It needs no judgement at all: the freedoms are whatever the nodes are, the assembly is mechanical, and there is nothing to choose.

The choice it needs is a real one, and its consequence is visible in the shape of the matrix it produces.

That pair is the flexibility method’s whole difficulty in one picture. Nothing about the beam changed between the two matrices, and one of them can be solved by hand while the other cannot. A stiffness assembly has no equivalent picture, because it has no equivalent choice: the freedoms are whatever the nodes are, and the only thing left to decide is the order they are written in.

By hand that trade is absurd, which is why the flexibility method dominated for a century. On a machine it inverts completely: a large system with no choices in it is exactly what a computer is good at, and a small system requiring judgement is exactly what it is not. The method that won is the one that was worse for humans, and it won on the property that made it worse.

The transition is visible in the literature of the 1950s. The aircraft industry needed to analyse structures with hundreds of redundancies, which is beyond any hand method, and matrix formulations of both approaches were written there — the flexibility one first, because it was familiar, and the stiffness one shortly after, because it was the one that automated. By the mid-1960s the stiffness method had absorbed both structural analysis and the finite-element method, which is the same assembly applied to elements that represent pieces of continuum instead of members.

Where the model stops #

The condition number is a physical statement. A frame with one member a thousand times stiffer than its neighbours produces a matrix whose eigenvalues span the same range, and the digits lost in the solve are the logarithm of that. Modelling a rigid link as a very stiff member — a common expedient — is a good way of destroying the precision of everything else in the model, and the symptom is not an error message but an answer that is wrong in its last several figures.

The element is the model. A frame element cannot represent a member that is warping, a panel zone that is shearing, or a joint region that is not a point. Every one of those needs either a different element or a spring, and choosing them is engineering rather than computation.

And the method is linear. Second-order effects, material yielding, contact and cracking all require the matrix to be reassembled as the structure deforms, which turns one solve into many and reintroduces every convergence question the hand methods had. The elegance above belongs to the linear case only.

The generalisation #

The habit worth taking from this is about where the difficulty in a problem actually sits.

The stiffness method is not clever about structures. It is clever about indexing — about the fact that if each member’s contribution is written in terms of the freedoms it touches, then equilibrium and compatibility are both enforced by addition into a shared array, with no reasoning at each step. The engineering insight was fully available in 1900; what was missing was a reason to organise it that way, and the reason arrived with machines that were good at arrays and bad at judgement.

That pattern is general. A great many methods in this subject are a rearrangement of something already known into a form that a particular instrument is good at: graphic statics rearranged equilibrium for a drawing board, influence lines rearranged a search for a pair of dividers, moment distribution rearranged compatibility for a slide rule, and matrix methods rearranged all of it for a machine.

The one thing that has not been rearranged is the free body. Every entry in every element matrix above came from cutting something out and insisting the sums cancel, and no amount of bookkeeping has removed the need to have done it.