7.4 Tests before use

Source: Ideas from Computational Physics (Landau, Páez, Bordeianu), §7.4. Rewritten in original wording. Continues 7.3b Matrices in Python labs. Labs: 7.4b Tests-before-use labs.

Trust, but verify

A library solve or a hand-written Gauss–Jordan inverse returns numbers. Before you publish a trajectory or an eigenvalue, ask three questions:

1. Does $AA^{-1}$ look like $I$?
2. Does $A\mathbf x$ reproduce $\mathbf b$?
3. Does $A\mathbf v\approx\lambda\mathbf v$ for a claimed eigenpair?

Residuals near machine epsilon ($\sim 10^{-16}$ in double precision) mean the arithmetic succeeded. Large residuals mean a bug, a singular matrix, or an ill-conditioned problem — not “close enough for physics.”

A concrete $3\times 3$ test matrix

Take

$$A=\begin{pmatrix}4&-2&1\\3&6&-4\\2&1&8\end{pmatrix}.$$

An analytic inverse (useful as a gold standard) is

$$A^{-1}=\frac{1}{263}\begin{pmatrix}52&17&2\\-32&30&19\\-9&-8&30\end{pmatrix}.$$

Gauss–Jordan on the augmented block $[A\mid I]$ reduces the left half to $I$ and leaves a numerical inverse on the right. Compare that inverse to the analytic form, then form $AA^{-1}-I$ and report $\max|(\cdot)_{ij}|$. For this mild matrix you should see residuals at a few $\times 10^{-16}$ or better.

Solving $A\mathbf x=\mathbf b$ for known right-hand sides (e.g. $\mathbf b=(12,-25,32)^{\mathsf T}$ with exact $\mathbf x=(1,-2,4)^{\mathsf T}$) is another check: print $\mathbf x$ to many decimals and confirm $A\mathbf x-\mathbf b$ is tiny. Prefer elimination on $[A\mid\mathbf b]$ over forming $A^{-1}\mathbf b$ when you only need one solve — a full inverse does extra work and can magnify round-off.

Eigenpairs and complex roots

For

$$A_3=\begin{pmatrix}-2&2&-3\\2&1&-6\\-1&-2&0\end{pmatrix},$$

one real eigenpair is $\lambda=5$ with

$$\mathbf x_1=\frac{1}{\sqrt{6}}(-1,-2,1)^{\mathsf T}.$$

Verify $\|A_3\mathbf x_1-5\mathbf x_1\|$ at machine noise. Algebraically, a $2\times 2$ rotation-like matrix can produce a complex conjugate pair $\alpha\pm i\beta$ with eigenvectors $[1,\mp i]^{\mathsf T}$ — check by multiplying; you do not need a numerical eigensolver to confirm the algebra. Degenerate eigenvalues (e.g. $\lambda=-3$ with a two-dimensional eigenspace) require two independent vectors that both satisfy $(A+\lambda I)\mathbf v=\mathbf 0$; one vector is not enough to span the subspace.

Ill-conditioned warm-up: Hilbert

The Hilbert matrix $H_{ij}=1/(i+j+1)$ is famously poorly conditioned as $N$ grows. For $N=5$, if $\mathbf b$ is the first column of $H$, the exact solution of $H\mathbf x=\mathbf b$ is the first unit vector $\mathbf e_1$. Your solver should recover that — measure $\max\|\mathbf x-\mathbf e_1\|$, not a flaky signed-zero printout of components.

Quiz

Tests before use

1 / 4

A numerical inverse is most convincingly checked by…

Pick an answer to see if you are right.

Next: 7.4b Tests-before-use labs.