A05: Numerical error, conditioning, roots, and approximation

A. Prerequisites and observable goals#

Level: Degree core, with numerical-analysis extensions clearly separated.

B. The central question: what did the computer actually approximate?#

A mathematical real number and a floating-point value are different objects. A finite floating-point format represents a finite set of values with nonuniform spacing. Arithmetic usually rounds to a representable result. In common IEEE binary64 round-to-nearest arithmetic, the unit roundoff is $2^{-53}$ for the usual relative-error model away from underflow and overflow. The model $\operatorname{fl}(a\circ b)=(a\circ b)(1+\delta)$ with $|\delta|\le u$ has hypotheses; it is not a rule for division by zero, nonfinite results, or every exceptional value. C09 develops the representation itself.

Numerical error includes rounding and algorithmic approximation. Modeling error is different: an exact solution to a frictionless equation may still predict a real system badly. Measurement uncertainty is another source. Reporting fifteen decimal places does not distinguish these causes or establish accuracy.

C. Conditioning belongs to the problem; stability belongs to the algorithm#

For a differentiable scalar function $y=f(x)$ with nonzero relevant quantities, a small relative input perturbation approximately produces a relative output perturbation amplified by

This local relative condition number comes from $\Delta y\approx f'(x)\Delta x$, divided by $y$ and compared with $\Delta x/x$. At $x=0$ or $f(x)=0$, choose an appropriate absolute or differently scaled measure instead of dividing by zero.

For $f(x)=1/(1-x)$ near $x=0.99$, $\kappa=|x/(1-x)|=99$. Small relative input uncertainty can be greatly amplified even by an excellent algorithm. A backward-stable algorithm produces the exact answer to a nearby problem, according to an appropriate perturbation measure. If the original problem is ill-conditioned, a nearby problem can have a very different solution. Therefore stability does not guarantee small forward error without suitable conditioning.

For linear systems, a small residual $r=b-A\widehat x$ is useful evidence but does not by itself guarantee $\widehat x$ is close to the true $x$. An almost singular $A$ can turn a small residual into a large solution error. A scaled residual such as $\|r\|/(\|A\|\|\widehat x\|+\|b\|)$ describes a different property from forward error. Be explicit about which one you measured.

D. Algebra can repair a numerically poor expression#

Consider $\sqrt{x+1}-\sqrt{x}$ for large positive $x$. The two square roots are close; each may already be rounded, and subtraction can remove their common leading digits while retaining rounding error. Rationalize the expression:

The identity is exact on the stated real domain, but the two floating-point evaluation paths can behave very differently. Near $x=10^{16}$, the direct subtraction can round to zero, while the rationalized expression yields about $5\times10^{-9}$. The issue is not that subtraction is always forbidden: subtracting exact nearby values can be well behaved. It is that errors already present in the operands may dominate the tiny desired difference.

Gaussian elimination has a related concern. A tiny pivot can create large multipliers and cancellation. Partial pivoting chooses a largest-magnitude available entry in the column, reducing those multipliers. It is a standard practical safeguard, not a universal proof that every problem is well-conditioned or every possible matrix has a small growth factor. The exact rational solver in A01 deliberately teaches a different issue: nonzero pivots and reversible algebra, not floating-point robustness.

E. Root finding with a guarantee you can state#

If $f$ is continuous on $[a,b]$ and endpoint signs differ, the Intermediate Value Theorem guarantees at least one root inside. Bisection evaluates the midpoint and retains a sign-changing half interval. After $k$ halvings, the interval width is $(b-a)/2^k$, so its midpoint is within $(b-a)/2^{k+1}$ of some root in the retained interval. This is an absolute location bound, not a proof of uniqueness. Continuity matters: a sign change across a pole does not establish a zero.

Newton’s method uses the local tangent approximation,

Near a suitable simple root, it can converge rapidly. It can also divide by a small derivative, leave the domain, cycle, or converge to an unwanted root. For $f(x)=x^3-2x+2$, starting at zero gives $x_1=1$ and then $x_2=0$, a cycle. A safeguarded method combines a bracket with faster steps only when they satisfy safety conditions. The lab implements bisection because its invariant and failure modes are especially clear.

A stopping rule based only on $|f(x)|$ can mislead when the function’s scale or derivative is unusual. A location bound, a scaled residual, an iteration budget, and a stagnation check answer complementary questions. The lab rejects nonfinite values and reports inability to reach its tolerance rather than silently returning an arbitrary last iterate.

F. Differentiation, integration, and interpolation are different approximations#

Forward differences approximate $f'(x)$ by $[f(x+h)-f(x)]/h$, with truncation error typically $O(h)$ under sufficient smoothness. Central differences use $[f(x+h)-f(x-h)]/(2h)$ and typically have $O(h^2)$ truncation error. Rounding effects can grow like $u/h$ after cancellation and division, so making $h$ smaller indefinitely is not an improvement. A rough central-difference balance suggests a step scale proportional to $u^{1/3}$ times relevant problem scales, not a universal hard-coded step for all functions.

Symbolic differentiation manipulates expressions according to mathematical rules. Automatic differentiation propagates local derivatives through the executed computational graph using the chain rule, avoiding a finite-difference step. It still uses finite precision and differentiates the chosen program, including its branch path and primitive definitions; discontinuities and unsupported operations require attention. Finite differences approximate derivatives through nearby function evaluations and are valuable as independent gradient checks.

Composite Simpson integration combines parabolic local approximations and requires an even positive number of subintervals. For sufficiently smooth $f$ on a finite interval, its classical error bound uses a bound on $|f^{(4)}|$ and is $O(h^4)$. A smoothness-based bound is not justified across a singularity. Module 23 derives the rule and bounds; the lab compares it with exact cubic integrals.

Interpolation constructs a function matching given data points. For distinct nodes, a polynomial of degree at most $n-1$ through $n$ scalar values is unique, but evaluating an expanded high-degree polynomial may be unstable and extrapolation may be poor. Piecewise low-degree interpolation and barycentric polynomial evaluation are useful continuation topics. “It passes through every measured point” does not mean it predicts unmeasured behavior correctly.

G. Independent practice#

F foundation C standard course work T transfer S optional stretch L runnable lab or project. Attempt the complete set before opening its solutions.

H. Complete solutions: A05#

Compare the first step where your work differs. A corrected answer is not yet an independent success.

Open every solution

"""A05.05: bounded numerical algorithms with explicit stopping conditions."""
import math
from collections.abc import Callable

def bisection(f: Callable[[float], float], a: float, b: float,
              xtol: float = 1e-10, max_steps: int = 200) -> float:
    if not all(math.isfinite(x) for x in (a, b, xtol)) or a >= b or xtol <= 0:
        raise ValueError("finite a < b and positive tolerance required")
    if type(max_steps) is not int or max_steps <= 0:
        raise ValueError("positive step budget required")
    fa, fb = f(a), f(b)
    if not math.isfinite(fa) or not math.isfinite(fb):
        raise ValueError("nonfinite endpoint value")
    if fa == 0: return a
    if fb == 0: return b
    if (fa > 0) == (fb > 0):
        raise ValueError("opposite endpoint signs required")
    for _ in range(max_steps):
        mid = a / 2 + b / 2
        fm = f(mid)
        if not math.isfinite(fm):
            raise ValueError("nonfinite function value")
        if fm == 0 or b - a <= 2 * xtol:
            return mid
        if mid == a or mid == b:
            raise ArithmeticError("floating-point stagnation")
        if (fm > 0) == (fa > 0):
            a, fa = mid, fm
        else:
            b, fb = mid, fm
    raise RuntimeError("iteration budget exhausted")

def central_difference(f: Callable, x: float, h: float) -> float:
    if not math.isfinite(h) or h <= 0 or x + h == x or x - h == x:
        raise ValueError("positive representable step required")
    return (f(x + h) - f(x - h)) / (2 * h)

def simpson(f: Callable, a: float, b: float, n: int) -> float:
    if type(n) is not int or n <= 0 or n % 2:
        raise ValueError("positive even subinterval count required")
    if not math.isfinite(a) or not math.isfinite(b):
        raise ValueError("finite endpoints required")
    h = (b - a) / n
    terms = [f(a), f(b)]
    terms += [(4 if i % 2 else 2) * f(a + i * h) for i in range(1, n)]
    if any(not math.isfinite(x) for x in terms):
        raise ValueError("nonfinite integrand sample")
    return h * math.fsum(terms) / 3

if __name__ == "__main__":
    print(bisection(lambda x: x * x - 2, 0, 2))
    print(central_difference(math.sin, 0.3, 1e-5))
    print(simpson(lambda x: x**3, 0, 1, 10))

I. Retention and readiness#

Retrieve limit versus point value, the Intermediate Value Theorem’s continuity condition, and units. Advance when a reported numerical answer includes its problem, approximation method, stopping rule, and an independent check. Further numerical analysis should cover backward error, robust factorizations, adaptive methods, and stability regions in greater depth.