Numerical Error and Algorithm Quality

In pure mathematics, solutions live in the serene continuum of real numbers $\mathbb{R}$. We manipulate infinite Taylor series, take limits as step sizes $\Delta x \to 0$, and write down exact roots with complete confidence. But the moment we ask physical hardware to calculate, we step into a world of unavoidable physical constraints: numbers must fit into finite bit registers (typically 32 or 64 bits), infinite processes must terminate in finite steps, and continuous functions must be represented on discrete grids.

This note establishes the foundational vocabulary and mental models for computational mathematics. Across subsequent topics, root finding, polynomial interpolation, and numerical calculus, we will repeatedly evaluate algorithms through four essential virtues:

1. Accuracy: How close is the computed approximation to mathematical truth?
2. Efficiency: What is the computational cost (work, memory, iterations) required to reach a target tolerance?
3. Stability: Does the algorithm prevent small representation perturbations from exploding into catastrophic errors?
4. Robustness: Does the algorithm behave predictably across diverse inputs and fail gracefully when assumptions break down?

From phenomena to computation

Before writing an algorithm or writing code, we must understand the transformation that carries an observed physical phenomenon into machine memory. This pipeline spans four distinct layers: the physical phenomenon, the mathematical model, the mathematical problem, and the digital computation.

A conservation principle might produce a partial differential equation; a balance of forces might yield a system of nonlinear algebraic equations. The modelling step decides what physics to retain and what to neglect. Crucially, solving these equations with infinite precision does not eliminate modelling error: an analytical solution to an incomplete model remains an incomplete description of reality.

Concrete mathematical problems include finding a root $f(x)=0$, integrating a differential equation, or solving a linear system $Ax=b$. A problem statement is incomplete without its domain and underlying hypotheses: existence, uniqueness, differentiability, and nonzero denominators matter just as much as the algebraic formula itself.

This definition separates the exact continuous map $F$ from what silicon hardware can execute. While real numbers $\mathbb{R}$ contain uncountably infinite information, a floating-point register holds only finitely many bits. A computation therefore produces an approximation $\hat{y}$ rather than the idealized $y$.

A formula describes a static mathematical relationship; a computation describes an active, ordered process. Two programs implementing the exact same mathematical formula can produce dramatically different answers because they sequence operations differently or suffer from different numerical cancellation paths.

The discipline asks far more than “what number did the computer output?” It investigates whether the approximation converges to the true solution, how fast the error decays, how much work is required, how representation errors propagate, and where the method can break down.

Method, algorithm, and implementation

To reason clearly about software failures and numerical bugs, we distinguish between three layers of realization: the abstract method, the executable algorithm, and the machine implementation.

A method can be based on iteration, discretisation, interpolation, linearisation, or recurrence. Its mathematical formulation identifies both the tuning parameter (such as step size $h$, polynomial degree $n$, or iteration count $k$) and the theoretical sense in which the error is expected to shrink.

A single numerical method often admits multiple competing algorithms. For instance, evaluating an $n$-th degree polynomial can be done by computing powers individually or via Horner’s nested multiplication. While mathematically equivalent in infinite precision, their computational complexities and roundoff propagation profiles differ substantially.

def numerical_algorithm(data, tolerance, max_steps):
    state = initialise(data)
    for step in range(max_steps):
        new_state = update(state, data)
        if error_indicator(new_state, state) <= tolerance:
            return new_state
        state = new_state
    raise RuntimeError("requested reliability was not established")

Two foundational principles are illustrated here: stopping criteria are an integral part of an algorithm rather than an afterthought, and hitting max_steps is not a numerical answer. It is proof that the requested accuracy was not certified.

“The algorithm works” is too vague a claim. A theorem about an exact method does not guarantee that its floating-point implementation will succeed, and an error-free execution run does not prove that the underlying problem is well-conditioned.

These four properties interact constantly. A fast algorithm that frequently diverges is useless; an accurate algorithm that takes exponential time is impractical; and a stable algorithm applied to an ill-conditioned problem will still produce large forward errors because the problem itself is hyper-sensitive.

Accuracy: what “close” means

Every numerical result can inherit error from the model, the approximation, the machine arithmetic, and earlier computational steps. Reducing one source does not automatically reduce the others.

Truncation extends beyond Taylor series to any method with a finite budget: a fixed number of Newton iterations, a quadrature rule with $N$ nodes, or a mesh of spacing $h$. The Taylor remainder $R_n(x)$ below is one exact description of truncation error. Its order tells how the discarded part changes as $n$ or $h$ is refined.

Computers store only finitely many significant digits, so each represented number and arithmetic result can be rounded. Two practical failure modes deserve particular attention: catastrophic cancellation, when nearly equal numbers are subtracted and reliable digits are lost, and overflow, when an intermediate magnitude lies outside the available floating-point range. The stable quadratic formula later in this note illustrates how to avoid catastrophic cancellation.

Absolute error has the same units as the quantity. Relative error is dimensionless and reports error compared with the scale of the answer. Neither is universally superior: relative error is undefined at $x_{\mathrm{true}}=0$ and can be misleading when zero is a meaningful reference point.

The residual is often computable even when the exact error is not. For a linear system, $r = b - A\hat{x}$ measures how well the computed vector satisfies the equations. A small residual implies a small forward error only when the problem is not too sensitive.

Taylor series and order

The notation makes the source of the approximation explicit: $P_n(x)$ is the quantity we compute, while $R_n(x)$ is what was discarded. The displacement from the expansion point is $x-x_0$. When only its size matters, write $h = |x-x_0|$. If $f^{(n+1)}$ stays bounded near $x_0$, then $R_n(x) = O(h^{n+1})$. The first omitted power of $h$ determines the local order of the approximation.

On logarithmic axes, an error law $E(h) \approx C h^p$ appears as a line with slope $p$. A straight line in a log-log plot supports an asymptotic error model only over the tested range. Coarse discretisation may not yet be asymptotic, while excessive refinement may expose round-off or modelling errors.

Reporting digits without an error scale is incomplete. A more meaningful result is “$x_{\mathrm{approx}}$ with estimated relative error below $10^{-6}$,” together with the assumptions behind that estimate.

Efficiency: the cost of obtaining accuracy

If evaluating $f$ requires a large simulation, the number of function evaluations may matter more than the number of scalar additions. On modern hardware, memory traffic and parallel communication can dominate arithmetic cost.

Useful distinctions include operation complexity, storage complexity, per-step cost versus the number of required updates, convergence rate, work-precision efficiency (total work to reach a target error), and scalability.

Suppose a one-dimensional discretisation uses $N \approx 1/h$ degrees of freedom, costs $O(N)$ work, and has error $E(h) = O(h^p)$. To achieve $E(h) < \varepsilon$, we need roughly

This is a work-precision law. Increasing the order $p$ can reduce the cost of high accuracy dramatically, provided the higher-order method has reasonable constants and its smoothness assumptions are satisfied.

def horner(coefficients, x):
    value = 0.0
    for coefficient in reversed(coefficients):
        value = coefficient + x * value
    return value

Efficiency is best understood as the computational cost required to achieve a target level of confidence, rather than raw wall-clock time alone. A fast calculation that frequently fails, requires repeated retries, or cannot certify its error bounds is rarely economical in practice.

Big-O as a shared language

Big-O is an eventual upper bound up to a constant. It is not an equation for an exact value, and it does not say that two functions have the same leading constant. The statement $f(t) = O(g(t))$ can remain true when $g$ is a loose upper bound.

For a tighter classification, $f = \Theta(g)$ means both an asymptotic upper and lower bound. The notation $f = o(g)$ means $f/g \to 0$, so $f$ is asymptotically smaller than $g$.

The grammar is shared, but the interpretation changes. In cost analysis, smaller growth is desirable. In an error law $O(h^p)$ with $h \to 0$, larger $p$ is usually desirable because the error decays faster.

If $r(h) = O(h^p)$ and $s(h) = O(h^q)$ as $h \to 0$, then

The lower power normally dominates a sum near zero. For example, $3h^2 + 7h^3 = O(h^2)$. Cancellation can improve the order, but it must be shown. It cannot be assumed from the two separate bounds.

The formal order of an iteration is a preview for later root-finding methods, not a Week 1 requirement. Let $e_k = x_k - x^*$. If

for $\mu > 0$, the iteration has order $p$. For $p = 1$, convergence is linear when $0 < \mu < 1$. The case $p = 2$ is quadratic convergence. This definition differs from the discretisation statement $E(h) = O(h^p)$, even though both use the word “order.”

Big-O deliberately hides constants and the point at which asymptotic behaviour begins. A method with cost $1000n$ can be slower than one with cost $n^2$ over the entire practical range. Likewise, an $O(h^4)$ method with a large error constant may initially be less accurate than an $O(h^2)$ method.

Stability: controlling perturbations

For a scalar map $y = F(x)$, a local relative condition number is $\kappa(x) = |x F'(x) / F(x)|$ when the expression is defined. Roughly,

The condition number describes the problem before an algorithm is chosen. No algorithm can recover information that the represented input does not contain.

The phrase “does not amplify error” is useful intuition, but it needs a reference scale. Even a stable algorithm can exhibit a large forward error on an ill-conditioned problem, because the exact problem itself amplifies input uncertainty.

Forward stability directly bounds the difference between computed and exact output. Backward stability interprets the computed output as the exact solution of a nearby problem. Mixed analysis combines input and output perturbation bounds when a purely forward or backward statement is inconvenient. Backward stability is often powerful because it separates responsibilities: the algorithm introduces only a small input perturbation, and the condition number predicts how that perturbation affects the output.

For $ax^2 + bx + c = 0$, the textbook formula

can subtract nearly equal numbers when $b > 0$ and $4ac$ is small compared with $b^2$. The small root may lose most of its significant digits. A stable strategy computes the large root without cancellation and obtains the other from $x_1 x_2 = c/a$:

from math import copysign, sqrt

def stable_quadratic_roots(a, b, c):
    discriminant = b * b - 4 * a * c
    q = -0.5 * (b + copysign(sqrt(discriminant), b))
    return q / a, c / q

The mathematical formula has not changed. The route through intermediate values has. This is the practical use of stability analysis: identify where information is lost, and reorganise the computation so that the machine solves a nearby problem accurately.

If one step transforms an error approximately as $e_{k+1} \approx G'(x_k) e_k$, repeated steps multiply the amplification factors. Stability therefore concerns a process, not only an isolated arithmetic operation. The unstable recurrence above is one instance of this mechanism.

Conditioning describes sensitivity of the problem itself; numerical stability describes additional error introduced by its algorithm.

Robustness: reliable behaviour across cases

Robustness is wider than stability. A stable update formula may still fail because an initial guess lies outside its basin of attraction, a derivative vanishes, a stopping test is misleading, or an intermediate value leaves the representable range.

Useful dimensions include domain checks, dependence on the starting guess, mixed-magnitude data, behaviour near difficult parameter regimes, termination tests that distinguish success from stagnation, and whether the method remains meaningful under available arithmetic.

A robust algorithm may combine methods. For root finding, a bracketed method offers dependable interval reduction, while a derivative-based step can be much faster near the root. A hybrid can accept the fast step when it stays inside the bracket and fall back to bisection otherwise. The safeguard adds work per step but improves the range of inputs for which the procedure can be trusted.

Useful safeguards include checking domains and finite values before an update, scaling variables to comparable magnitudes, combining absolute and relative stopping tolerances, monitoring residual reduction as well as step size, imposing iteration limits with explicit failure states, and switching methods when progress stalls.

A case study connecting all four properties

The derivation of finite difference formulas begins with Taylor’s theorem:

where $\xi$ lies between $x$ and $x+h$. Rearranging gives

Thus the forward difference is first-order accurate. If $f$ is smooth, the centred difference

has truncation error $O(h^2)$. It is tempting to make $h$ as small as possible. In floating-point arithmetic, however, the numerator subtracts nearby function values and the division by $h$ magnifies their representation errors. A simplified total-error model is

where $u$ is the unit round-off. Refinement reduces truncation error but can increase round-off error. The total error has a useful intermediate scale.

Balancing the two leading terms suggests an optimal scale of order $h = O(u^{1/3})$ for this model. That conclusion is more useful than the slogan “smaller $h$ is better”: it combines truncation analysis, finite-precision stability, and the practical choice of a parameter.

A reusable framework

For a new numerical problem, keep the layers distinct:

1. Specify the problem: inputs, desired output, domain, assumptions, and a meaningful scale for error.
2. Assess conditioning: how uncertainty in the data can affect the exact answer.
3. Choose a method: identify its approximation parameter and convergence assumptions.
4. Design the algorithm: representation, updates, stopping tests, safeguards, and failure states.
5. Analyse accuracy, efficiency, and stability as three separate questions.
6. Test robustness on difficult scales and adversarial initialisations rather than typical inputs alone.
7. Report the approximation together with residuals, error estimates, tolerances, and relevant assumptions.

Numerical analysis studies the path from a mathematical problem to a trustworthy computed result. Big-O notation links several parts of that path: it describes how resource cost grows, how approximation error decays, and how iterative errors converge. Asymptotic order is only one dimension of quality. Constants, conditioning, floating-point stability, safeguards, and the intended use of the answer determine whether a method is genuinely useful.

The next note applies this language to bisection and fixed-point iteration. This page belongs to the Computational Mathematics reading path.