14 Numerical Experiments and Their Prerequisites

14.1 Learning objectives

After completing this unit, you will be able to:

choose numerical methods whose mathematical prerequisites you understand;
record invariants, tolerances, iteration limits, and stopping criteria;
distinguish problem-formulation error, discretization error, rounding error, and implementation error;
check results through substitution, residuals, bounds, or comparison methods;
run convergence experiments without treating them as general proofs; and
reject library usage that hides assumptions about the domain.

This unit offers several routes. The root-finding route requires A30 and the concept of a continuous function. The quadrature route requires B30, the linear systems route requires B40, and the differential equations route requires B70. Defer any route whose prerequisites you have not yet met.

The B80 core boundary

Only the sections on methods and theorems, A30 bisection, the SciPy comparison, and sources of error, together with Exercises 1 and 2 and SciPy exercise s01, form this unit’s B80 core. The B30, B40, and B70 routes remain available as extensions, but must not count as evidence of B80 completion before their respective prerequisites have been met.

14.2 Methods come with theorems

A software command does not guarantee that a problem satisfies a method’s assumptions. Bisection, for example, maintains an interval on which the function has opposite signs at the two endpoints. Guaranteeing a root in that interval requires continuity. A program can check the endpoint signs; it cannot infer continuity from a few samples. This unit’s implementation also rejects nonfinite endpoints, midpoints, or function values.

Every numerical experiment states:

the problem and domain;
the mathematical assumptions;
the method and representation;
the parameters and stopping criteria;
the invariants being checked;
the error measure or residual; and
the limits of the conclusions.

14.3 A30 route: bisection

To find a root of $f(x)=0$ , start with $[a,b]$ such that $f(a)f(b)<0$ . Each step selects the midpoint and keeps the half-interval that still has a sign change.

If the initial width is $w_0=b-a$ , after $k$ steps the width is $w_k=w_0/2^k$ . This relationship provides an interval-based error bound, rather than relying only on a small value of $f(x)$ .

from source.code.unit11_numerical import bisection

result = bisection(lambda x: x*x - 2.0, 1.0, 2.0, width_tolerance=1e-12)
(result.midpoint, result.width, result.iterations)

Listing 14.1

(1.414213562372879, 9.094947017729282e-13, 40)

That value is an approximation. The exact statement $x=\sqrt{2}$ comes from the definition of the positive root of $x^2=2$ , not from the printed digits.

14.4 A30 route: an executable SciPy comparison

The original implementation makes each bisection step readable. A second implementation, scipy.optimize.root_scalar, provides an independently developed comparison. Both take a specified function, initial interval, method, and iteration limit, but their tolerance parameters must not be treated as equivalent. original_width_tolerance bounds the final interval width in the original implementation; scipy_xtol and scipy_rtol specify absolute and relative closeness criteria for the root returned by SciPy. Their numerical values may coincide in this example, but their meanings remain distinct and are recorded separately.

from source.code.unit11_numerical import (
    SCIPY_VERSION,
    compare_bisection_with_scipy,
    square_minus_two,
)

comparison_u11 = compare_bisection_with_scipy(
    square_minus_two,
    1.0,
    2.0,
    original_width_tolerance=1e-12,
    scipy_xtol=1e-12,
)

assert comparison_u11["scipy_version"] == SCIPY_VERSION
assert comparison_u11["verification"]["strict_sign_change"]
assert comparison_u11["scipy"]["converged"]
assert comparison_u11["verification"]["roots_inside_initial_bracket"]
assert comparison_u11["verification"]["roots_agree"]

(
    comparison_u11["scipy_version"],
    comparison_u11["original"]["midpoint"],
    comparison_u11["scipy"]["root"],
    comparison_u11["verification"]["absolute_root_difference"],
)

Listing 14.2

('1.15.2', 1.414213562372879, 1.4142135623724243, 4.547473508864641e-13)

The record stores the version of SciPy actually imported, not just the library name. strict_sign_change checks $f(a)f(b)<0$ before the call. converged and flag record SciPy’s report. Both roots must remain inside the initial interval, and their difference must not exceed the recorded agreement tolerance.

Agreement between two implementations improves the chance of detecting local coding errors. It does not prove that the function is continuous or that the interval contains exactly one root, nor does it turn approximate digits into an exact value. Continuity and the existence guarantee rest on the assumptions and the Intermediate Value Theorem; the version and convergence checks provide computational evidence.

14.5 B30 route: quadrature and refinement

B30 prerequisite gate - deferred extension

The quadrature section and Exercise 3 require B30. Both can be used for enrichment once that prerequisite has been met, but neither counts as B80 core content or evidence of B80 completion.

The composite trapezoidal rule replaces the curve on each subinterval with a line segment. A refinement experiment uses $n,2n,4n,\ldots$ subintervals and records changes in the result. Values that appear stable support a conjecture of convergence, but do not replace error analysis.

For a linear function, the trapezoidal rule is exact in exact arithmetic. This case provides a stronger implementation check than an arbitrarily chosen reference number.

14.6 B40 route: linear systems and residuals

B40 prerequisite gate - deferred extension

The linear systems section and Exercise 4 require B40. Neither counts as B80 core content or evidence of B80 completion.

If a program proposes $\widehat{x}$ as a solution of $Ax=b$ , compute the residual $r=b-A\widehat{x}$ . A small residual means that the equation is nearly satisfied on the scale being used. A small residual does not always mean a small error in the solution; the relationship depends on the matrix’s conditioning.

For a two-by-two system with rational coefficients, an exact solution using Fraction can provide an independent comparison for the floating-point route.

14.7 B70 route: Euler steps

B70 prerequisite gate - deferred extension

The differential equations section and Exercise 5 require B70. Neither counts as B80 core content or evidence of B80 completion.

For the problem $y'=g(t,y)$ , an Euler step uses

$y_{k+1}=y_k+h g(t_k,y_k).$

Reducing $h$ and comparing with a known solution can reveal an error pattern in the example. This tests the implementation and the behavior of that case; it does not prove an order of convergence for every function $g$ .

14.8 Four sources of disagreement

When a result differs from expectations, distinguish:

problem-formulation error: the model or data does not represent the question;
discretization error: the method replaces a continuous problem with a finite one;
rounding error: floating-point representation changes the operations; and
implementation error: the code does not carry out the intended method.

Changing a tolerance addresses only some of these sources. A report must therefore not label a single number “numerical error” without defining it.

14.9 Experiment artifacts

Run the local examples whose prerequisites you have met:

python source/code/unit11_numerical.py --output output/unit11-results.json

The output records the final bisection interval, the quadrature refinement table, system solutions and residuals, and Euler-step errors. The function identities, domains or initial data, tolerances, refinement grids, step counts, and reference values used by the checks are stored alongside the results.

The record also includes the SciPy comparison, the SciPy version, interval checks, convergence, agreement between results, and the limits of the evidence. The quadrature, linear systems, and Euler tables are retained as extension regression checks; their presence in the JSON does not make them part of the B80 core. The curriculum_routes field makes this separation explicit.

14.10 Exercises

14.10.1 Exercise 1 - bisection invariants

Name two invariants that should be checked at each bisection step, and explain why the midpoint value alone is insufficient.

Hint

Consider the sign change and the shrinking interval.

Answer and solution

The interval endpoints must continue to bracket a sign change, and the width must halve at each step. The midpoint value alone does not establish that a root remains bracketed or that the implementation selected the correct half.

14.10.2 Exercise 2 - an iteration bound

What is the minimum number of steps guaranteeing an interval width of at most $10^{-6}$ if the initial width is 1?

Hint

Find $k$ such that $2^{-k}\leq 10^{-6}$ .

Answer and solution

Since $2^{19}=524288<10^6$ and $2^{20}=1048576>10^6$ , 20 steps are needed. After 20 steps, the width is $2^{-20}<10^{-6}$ .

14.10.3 SciPy exercise - two implementations for a root of a cubic

Use compare_bisection_with_scipy to find a root of $f(x)=x^3-x-2$ on $[1,2]$ with a width tolerance of $10^{-10}$ . Verify the sign change, SciPy version, convergence report, inclusion of both results in the interval, and agreement between the roots. Explain one matter that still requires a mathematical argument.

Hint

Define the Python function first. The comparison result has fields named endpoint_values, scipy_version, scipy, original, and verification.

Self-check

from source.code.unit11_numerical import SCIPY_VERSION, compare_bisection_with_scipy

def cubic_function_u11(x):
    return x**3 - x - 2

cubic_result_u11 = compare_bisection_with_scipy(
    cubic_function_u11,
    1.0,
    2.0,
    original_width_tolerance=1e-10,
    scipy_xtol=1e-10,
)

assert cubic_result_u11["endpoint_values"][0] < 0
assert cubic_result_u11["endpoint_values"][1] > 0
assert cubic_result_u11["scipy_version"] == SCIPY_VERSION
assert cubic_result_u11["scipy"]["converged"]
assert cubic_result_u11["verification"]["roots_inside_initial_bracket"]
assert cubic_result_u11["verification"]["roots_agree"]
assert (
    cubic_result_u11["verification"]["absolute_root_difference"]
    <= cubic_result_u11["verification"]["agreement_tolerance"]
)

Answer and solution

Solution 14.1.

from source.code.unit11_numerical import SCIPY_VERSION, compare_bisection_with_scipy

def solution_cubic_function_u11(x):
    return x**3 - x - 2

cubic_solution_u11 = compare_bisection_with_scipy(
    solution_cubic_function_u11,
    1.0,
    2.0,
    original_width_tolerance=1e-10,
    scipy_xtol=1e-10,
)

cubic_verification_u11 = cubic_solution_u11["verification"]
assert cubic_solution_u11["method"] == "bisect"
assert cubic_solution_u11["scipy_version"] == SCIPY_VERSION
assert cubic_solution_u11["scipy"]["flag"] == "converged"
assert cubic_verification_u11["strict_sign_change"]
assert cubic_verification_u11["roots_inside_initial_bracket"]
assert cubic_verification_u11["roots_agree"]

Since $f(1)=-2$ and $f(2)=4$ , the endpoint values have opposite signs. Assuming that $f$ is continuous, the Intermediate Value Theorem guarantees at least one root in the interval. The two implementations agree within the recorded tolerance, and SciPy reports convergence. However, the two programs prove neither continuity nor uniqueness of the root; both require separate mathematical arguments.

14.10.4 Exercise 3 - a quadrature check

Curriculum status: deferred B30 extension; not B80 core.

Apply the trapezoidal rule to $f(x)=3x+2$ on $[0,4]$ with one subinterval. Compare with the exact integral.

Hint

The trapezoid’s area is its width times the average endpoint height.

Answer and solution

The endpoint values are 2 and 14, so the trapezoid’s area is $4(2+14)/2=32$ . The exact integral is $[3x^2/2+2x]_0^4=24+8=32$ . This exact agreement is expected because the graph of a linear function on the interval is itself a line segment.

14.10.5 Exercise 4 - a residual is not a solution error

Curriculum status: deferred B40 extension; not B80 core.

Explain why a small residual does not automatically guarantee that an approximate solution is close to the exact solution.

Hint

Consider a system that is highly sensitive to small changes in the right-hand side.

Answer and solution

The residual measures how well $\widehat{x}$ satisfies the given equation. For an ill-conditioned matrix, a small change in the right-hand side can correspond to a large change in the solution. A residual must therefore be interpreted alongside scale and conditioning information, not directly as the error in the solution.

14.10.6 Exercise 5 - an Euler experiment

Curriculum status: deferred B70 extension; not B80 core.

For $y'=y$ , $y(0)=1$ , compare the errors at $t=1$ for $h=1/10$ and $h=1/20$ . What conclusions are justified?

Hint

Use the reference solution $e^t$ , and do not turn two cases into a theorem.

Answer and solution

Euler’s method gives $(1+h)^{1/h}$ . For $h=0.1$ the value is approximately 2.5937; for $h=0.05$ it is approximately 2.6533; whereas $e\approx2.7183$ . The error decreases in these two cases. This supports the expected behavior and checks the implementation, but two step sizes do not prove convergence or an error order in general.

14.11 Summary

Numerical methods come with assumptions and theorems.
The original implementation and scipy.optimize.root_scalar(method="bisect") provide an A30 cross-check that must record the version, interval, tolerances, convergence, and agreement between results.
Stopping criteria must be specified before inspecting the results.
Residuals, interval bounds, and exact cases provide different checks.
Problem-formulation, discretization, rounding, and implementation errors must not be conflated.
Refinement experiments provide computational evidence; proving convergence requires a mathematical argument.
B30 quadrature, B40 linear systems, and B70 Euler steps remain extensions and must not count as B80 core content.