Systems of linear equations

The formula

a11·x1 + a12·x2 + ... + a1n·xn = b1
a21·x1 + a22·x2 + ... + a2n·xn = b2
  ...
an1·x1 + an2·x2 + ... + ann·xn = bn

n statements about n unknowns, each one a straight line, plane, or hyperplane. The solution is where all of them intersect.

A·x = b — the same thing in matrix shorthand. A holds the physics you know (stiffnesses, resistances, geometry), b holds what you applied or measured, x is what you're after.

det(A) ≠ 0 — the square system has exactly one solution. Determinant zero means the equations are not independent: either no solution or infinitely many.

x1 = (b1·a22 − b2·a12) / (a11·a22 − a12·a21) — Cramer's rule for the 2×2 case. Fine on paper at n = 2 or 3, ruinously expensive beyond that.

flops ≈ (2/3)·n³ — the cost of Gaussian elimination. Double the unknowns, pay eight times the arithmetic. A 1,000-unknown dense system runs about 667 million operations ^[1].

κ(A) = ‖A‖ · ‖A⁻¹‖ — the condition number. Expect to lose roughly log10(κ) digits of accuracy to the problem itself, before your algorithm loses any ^[7].

Where you meet it

Circuit bench. Kirchhoff's current law at every node of a network hands you one linear equation per node. Mesh analysis does the same per loop. SPICE builds and solves A·x = b at every timestep of every simulation you've ever run — that's the whole engine.
Test stand calibration. A six-axis force balance doesn't read out clean channels; every gage picks up a little of every load. The calibration matrix that untangles cross-coupling comes from applying known loads and solving the resulting linear system. Same story for any multi-channel sensor cal with interaction terms.
Structures review board. The FEA slide showing a million-element model is showing you K·u = F — stiffness matrix times displacements equals applied loads. The solver's job is one enormous sparse linear system, and the runtime the analyst is apologizing for is the n³-flavored cost of solving it.
Truss statics. Sum of forces at each joint, two equations per pin in 2D. A statically determinate truss is exactly the case where the system is square and the determinant is nonzero.

How it works

Every practical method is some form of elimination: use one equation to knock an unknown out of the others, repeat until you're down to one equation in one unknown, then substitute back up the chain. Modern solvers organize this as an LU factorization — split A into a lower and an upper triangular factor once, then solve L·y = b and U·x = y by cheap substitution. The payoff shows up when you have many right-hand sides against one matrix: factor once at (2/3)·n³, then each new b costs only about 2·n².

A worked case, mesh currents on a two-loop circuit:

 5·I1 − 2·I2 = 10
−2·I1 + 3·I2 = 0

Determinant is 15 − 4 = 11, so I1 = 30/11 ≈ 2.727 A and I2 = 20/11 ≈ 1.818 A. Plug them back in; both equations close.

The gotcha that bites working engineers is conditioning. Take this pair:

x1 +        x2 = 2
x1 + 1.0001·x2 = 2.0001

Solution: x = (1, 1). Now perturb the last right-hand-side entry from 2.0001 to 2.0002 — one part in twenty thousand, well inside the noise of any real measurement — and the solution jumps to x = (0, 2). Nothing is wrong with the algebra. The two lines are nearly parallel, the intersection slides a long way for a tiny nudge, and the condition number κ ≈ 4.0·10⁴ says so up front: expect around 4 to 5 digits gone.

Three mistakes to retire:

Judging health by the determinant. Scale a matrix by 10 and its determinant changes by 10ⁿ while the geometry of the problem hasn't moved. Turing made this point in 1948 with matrices sharing one determinant and wildly different condition numbers ^[5]^[6]. Check κ, not det.
Trusting a small residual. ‖A·x − b‖ can be tiny while x is far from the true answer — that is precisely what ill-conditioning means. A calibration fit that "closes to 0.01%" can still carry garbage coefficients.
Computing A⁻¹ explicitly. Inverting and multiplying costs roughly three times the arithmetic of a factorization and is generally less accurate. Every serious library's solve routine factors; none of them invert.

Two more habits worth keeping. Pivoting — swapping rows to divide by the largest available element — is not optional; textbook elimination without it can fail on perfectly solvable systems. And in double precision you start with about 16 digits, so a κ of 10¹² (a 10×10 Hilbert-type fit matrix gets there easily — κ ≈ 1.6·10¹³ at n = 10) leaves you 3 or 4. That is why polynomial calibration fits go sideways at high order: the equations are fine, the matrix is nearly singular.

History

The method is older than almost anything else engineers use daily. The Nine Chapters on the Mathematical Art, compiled in Han-dynasty China with material dating to roughly 150 BC, devotes its eighth chapter — fangcheng, "rectangular arrays" — to exactly this: coefficients laid out on a counting board and eliminated column by column, in what a modern reader recognizes immediately as Gaussian elimination, some two thousand years before Gauss ^[1]^[2].

In Europe the method surfaced in Newton's notes of 1669–70, published in 1707 as Arithmetica Universalis, and from there it spread into the algebra textbooks as an unremarkable standard technique — Gauss himself later called it "common" ^[1]^[4]. Determinant-based solutions arrived in parallel: Maclaurin's posthumous 1748 Treatise of Algebra proved the 2×2 and 3×3 cases, and Gabriel Cramer published the general n×n rule in 1750 ^[2]^[3].

Gauss's name got attached almost by accident. Recovering the orbit of the asteroid Pallas from observations between 1803 and 1809, he ran least-squares fits that reduced to six equations in six unknowns, and around 1810 devised a compact notation for the symmetric elimination involved. Professional hand computers adopted that notation for a century of surveying and astronomy work, and by the 1950s the "common" method had become "Gaussian elimination" ^[1]^[2]^[4].

The modern chapter opens with Alan Turing. His 1948 paper "Rounding-off errors in matrix processes" recast elimination as the factorization A = LU, introduced the condition number, and showed that the determinant tells you nothing about conditioning — the framework every numerical library still lives in ^[5]^[6]. James Wilkinson's backward error analysis of 1961 finished the job, explaining why elimination with pivoting is trustworthy on real hardware ^[5]^[6].

Related tools

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →