[mitchell-tesch]

I was trying out solve(geometrical_non_linear=True) on a simple beam and it raised LinAlgError: Singular matrix:

ss = SystemElements(EI=1810, EA=460_000)
ids = ss.add_multiple_elements([[0, 0], [6, 0]], 20)
ss.add_support_hinged(1)
ss.add_support_roll(ss.id_last_node)
ss.q_load(q=-5, element_id=ids)
ss.point_load(ss.id_last_node, Fx=-150)
ss.solve(geometrical_non_linear=True)  # LinAlgError: Singular matrix

There turned out to be three bugs, one per commit.

1. The buckling eigen-solve. #138 replaced scipy.linalg.eigvals(k0, kg) with np.linalg.inv(kg) * k0. kg is always singular for beam elements (it has no axial terms), and * is an element-wise product rather than a matrix product. I've changed it to take the eigenvalues of np.linalg.solve(k0, kg) instead, which never inverts kg and still only needs numpy. One small behaviour change: a structure with nothing in compression now gets a buckling factor of inf. Before, the abs() gave it a finite number.

2. The sign of the geometric stiffness matrix. Once the solve ran, compression was reducing the moments and tension was increasing them. geometric_stiffness_matrix converts the matrix from Gavin's notes (linked in the code, y up) to anaStruct's y direction by multiplying by [1, -1, 1, 1, -1, 1]. That only flips the columns; it needs to flip the rows as well, so it's now an np.outer of that vector with itself. I also fixed a typo in row 4 (s1 * c2 should be s2 * c2).

3. Buckling factors below 1. det_linear_buckling builds k0 + kg with a second solve(). When the loads are already beyond the buckling load, that matrix is no longer positive definite, so validate() rejected it with a StabilityError ("Check your support conditions") and no factor was returned. Only the assembled matrix is needed, so it's now built with solve(naked=True). A member at 1.2 N_cr reports a factor of 0.833 (test_buckling_factor_below_one_is_reported, horizontal and vertical).

Results

For the beam above (6 m, EI = 1810 kNm², 5 kN/m, 20 elements), compared with the closed-form beam-column solution:

Axial force	Max moment (exact)	Max moment (this PR)	Buckling factor
0.3 N_cr compression	32.42	32.42	3.333 (exact 3.333)
0.6 N_cr compression	57.27	57.27	1.667 (exact 1.667)
0.3 N_cr tension	17.17	17.17	inf

On master all three raise. The existing buckling test for example 18 (a vertical column, tests/test_e2e.py) still passes with a factor of 600.

The new tests cover the matrix itself (tests/test_geometric_stiffness.py) and the closed-form cases above for horizontal and vertical members (tests/test_geometrical_non_linear.py). black, mypy, pylint and pytest all pass locally.

A member with no transverse load, or using discretize_kwargs, still fails because of a separate bug, which I've fixed in #491. I also noticed #482 covers this as part of a much larger change. This is just the minimal fix, in case it's easier to review on its own.

[mitchell-tesch]

I also have a small follow-up that speeds up validate() and the buckling solve using symmetric eigen-solvers (~2× faster second-order solves on larger models). I'll open it once this is merged.