Linear algebra with NumPy

Why this matters

A surprising amount of engineering and science ends up as a set of linear equations. Kirchhoff's laws for a circuit, force balance at the joints of a truss, mass balances in a chemical process, calibrating a spectrometer against several standards: each gives you n equations in n unknowns. By hand, three unknowns is tedious and ten is hopeless. With NumPy it is one line, and the harder part becomes writing the system down correctly and checking the answer.

Core idea: write it as A x = b

Take a two-loop circuit. A 12 V source drives loop 1, which shares a 4 Ω resistor with loop 2. Writing Kirchhoff's voltage law around each loop with mesh currents i1 and i2:

loop 1:  (2 + 4) i1 -       4 i2 = 12
loop 2:      -4 i1 + (4 + 6) i2 =  0

Collect the coefficients into a matrix A, the unknowns into a vector x, and the right-hand sides into b:

A = [[ 6, -4],      x = [i1,      b = [12,
     [-4, 10]]           i2]           0]

Each row of A is one equation, each column belongs to one unknown.

import numpy as np

A = np.array([[6.0, -4.0],
              [-4.0, 10.0]])
b = np.array([12.0, 0.0])

x = np.linalg.solve(A, b)
print(x)          # [2.7273 1.0909]  amps

Matrix products and the @ operator

A @ x is the matrix-vector product (np.dot does the same). Plain A * x is element-wise and broadcasts, which is almost never what you want here.

print(A @ x)      # [12.  0.]  reproduces b

Other tools you will use:

• A.T: transpose.
• np.linalg.det(A): determinant (zero means no unique solution).
• np.linalg.inv(A): the inverse. Avoid it for solving; solve is faster and more accurate.
• np.linalg.norm(v): the length of a vector.
• np.linalg.eig(A): eigenvalues and eigenvectors (vibration modes, principal stresses).

Always check: the residual

A solver always returns numbers. Whether they are trustworthy is your job. The residual r = A x - b should be tiny compared with b:

r = A @ x - b
print(np.linalg.norm(r))                 # about 1e-15
print(np.allclose(A @ x, b))             # True

If you typed a coefficient into the wrong row, the residual against your original equations will still be small (the solver solved what you gave it). So also sanity-check physics: currents should satisfy Kirchhoff's current law, forces should sum to zero, concentrations should be positive.

Conditioning: when small errors blow up

Some systems are ill-conditioned: tiny changes in b (measurement noise) cause huge changes in x. The condition number measures this:

print(np.linalg.cond(A))        # about 3.5: well conditioned

B = np.array([[1.0, 1.0],
              [1.0, 1.0001]])
print(np.linalg.cond(B))        # about 40,000: nearly singular

Rule of thumb: you can lose about log10(cond) significant digits. A condition number of 1e4 with data good to 3 digits means the answer may have none. Ill-conditioning usually signals two nearly identical equations, for example two calibration standards with almost the same concentration.

Overdetermined systems: least squares

With more equations than unknowns (more measurements than parameters) there is usually no exact solution. Least squares finds the x that minimises ||A x - b||.

Worked example: a calibration line

A sensor's output y (mV) is measured at known temperatures T (°C). Fit y = m T + c:

T = np.array([0.0, 10.0, 20.0, 30.0, 40.0])
y = np.array([0.02, 0.41, 0.79, 1.22, 1.60])

A = np.column_stack([T, np.ones_like(T)])   # columns: T and 1
coef, residuals, rank, sv = np.linalg.lstsq(A, y, rcond=None)
m, c = coef
print(round(m, 4), round(c, 4))             # 0.0397 0.014

The design matrix A has one column per parameter. The same trick fits any model that is linear in its parameters (polynomials, sums of known functions) and is the basis of curve fitting in Module 3.

SciPy adds more specialised tools (scipy.linalg.lu_factor to reuse a factorisation, sparse solvers in scipy.sparse.linalg for systems with millions of unknowns), but numpy.linalg covers most day-to-day work.

Common mistakes

• * instead of @. Element-wise multiplication silently gives the wrong numbers.
• Integer matrices. Build matrices with floats (6.0, not 6) so later edits do not truncate.
• Using inv(A) @ b. It works but is slower and less accurate than solve(A, b).
• Mixing up rows and columns. One row per equation, one column per unknown, in the same order every row.
• Ignoring LinAlgError: Singular matrix. It means the equations are not independent: a missing equation or a duplicated one.

Quick recap

• Write the problem as A x = b, then x = np.linalg.solve(A, b).
• Use @ for matrix products.
• Check np.allclose(A @ x, b) and the physics; check np.linalg.cond(A) when inputs are noisy.
• For more equations than unknowns, build a design matrix and use np.linalg.lstsq.

Try it

In the Do tab you will solve a three-loop circuit, verify it with the residual, and fit a calibration line with least squares.