Root-Finding Methods Compared: Convergence Analysis and the Colebrook-White Friction Factor

Every numerical-methods textbook proves that Newton-Raphson converges quadratically near a simple root, while bisection grinds along at its guaranteed-but-slow linear rate. This project tests those theorems against an equation that practicing engineers actually need solved: the Colebrook-White equation for the Darcy friction factor f in turbulent pipe flow,

1/√f = -2 log₁₀(ε/3.7D + 2.51/(Re√f)),

which is implicit in f — the unknown appears on both sides — so a fixed number of arithmetic operations cannot extract it. The study implements bisection, Newton-Raphson (with the analytically derived derivative), the secant method and the fixed-point iteration of Swamee and Jain's explicit approximation family in Python, then subjects each to a declared experimental protocol: 1,000 randomly drawn (Re, ε/D) pairs spanning Re = 4×10³ to 10⁸ and relative roughness 10⁻⁶ to 5×10⁻², root accuracy 10⁻¹⁰, and a cap of 100 iterations. For every run the project records iterations-to-convergence, final error, and failure cases, estimates the observed convergence order p from three consecutive errors, and explains every divergence. The deliverable is a reproducible notebook, error-and-order plots, and a written analysis connecting the observed behaviour to the underlying assumptions of each theorem.

The Colebrook-White equation defines the Darcy friction factor f for turbulent pipe flow, but because f appears inside its own logarithm, no closed-form rearrangement exists, so every evaluation of head loss in a real pipeline is secretly a root-finding problem. Undergraduate numerical-methods courses teach bisection, Newton-Raphson and secant almost exclusively on well-behaved polynomials, where convergence is quick and failure is rare; students rarely see how the same methods behave on an industrially important equation whose derivative is nearly flat at the root and whose parameters span five orders of magnitude. The result is a gap between theorem and practice: a student who can quote the quadratic convergence of Newton's method often cannot say how many iterations it will need on this equation, why a poor initial guess diverges, or which method a working engineer should trust in a spreadsheet. This project closes that gap for one equation with a measured, reproducible comparison, quantifying iterations, observed convergence order and failure rates across the physically meaningful range of Reynolds number and relative roughness.

Abstract. The Colebrook-White equation is the industry-standard correlation for turbulent friction in fully rough and transition pipes, adopted by codes and handbooks worldwide, yet it cannot be rearranged for the friction factor. This project compares four iterative root-finders — bisection, Newton-Raphson, secant and a fixed-point scheme — under a pre-declared protocol, quantifies their convergence orders empirically, and packages the analysis as a reproducible Python study.

Introduction. Iterative methods sit at the heart of applied mathematics, and the B.Sc syllabus meets them as theorems: the bisection convergence guarantee from the intermediate value theorem, Newton's quadratic order for simple roots with a good initial guess, the secant method's superlinear order (1+√5)/2. Textbook examples are polynomials where everything behaves. A famous implicit equation from fluid mechanics is a more honest opponent: extremely flat derivatives near the root, a logarithmic structure that punishes poor starting points, and physical parameter ranges that stretch over five orders of magnitude.

Existing work and the gap. Engineers sidestep the implicit equation with explicit approximations (Swamee-Jain, Haaland, and successors), and textbooks tabulate solutions. What standard treatment omits is a measured comparison: how many iterations each classic method really needs across the physical parameter space, where each fails, and whether the theoretical convergence orders survive contact with this equation. That measurement, with a stated protocol and reproducible code, is this project's contribution at the undergraduate level.

Proposed study. (1) Implement the four methods with careful attention to bracketing, derivative derivation and stopping rules. (2) Draw 1,000 random (Re, ε/D) points from the physical ranges by Latin-hypercube sampling. (3) Run each method to |Δf| < 10⁻¹⁰ or 100 iterations; record everything. (4) Estimate empirical convergence orders from error triples; compare against theory. (5) Analyse failures: Newton divergence from bad seeds, bisection bracket invalidity, fixed-point non-contraction. (6) Close with the engineering view — a worked example computing head loss for a campus water-supply line.

Feasibility. The mathematics needs only core calculus and numerical-methods course material; the coding is standard Python with NumPy and Matplotlib; no laboratory equipment is required. A personal laptop completes the full 4×1,000-run experiment in under a minute, and every figure and table regenerates from one notebook, making the study reproducible for the viva demonstration.

In scope. The single-equation root-finding problem for the Colebrook-White correlation: four iterative methods (bisection, Newton-Raphson with analytic derivative, secant, and one fixed-point form), a randomised parameter study over Re = 4×10³–10⁸ and ε/D = 10⁻⁶–5×10⁻², empirical convergence-order estimation, failure analysis, and one worked engineering example converting the computed friction factor into head loss for a water line. Implementation is Python with NumPy and Matplotlib; all experiments run from one notebook.

Out of scope. Systems of nonlinear equations, other pipe-flow correlations beyond comparison references (Haaland is cited for context only), non-turbulent (laminar) flow except as a boundary note, minor-loss and network calculations, and parallel or GPU acceleration. The study measures the mathematics, not fluid-dynamics instrumentation.

The study follows a fixed experimental protocol declared before any run.

Phase 1 — formulation (weeks 1–2). Write the residual g(f) = 1/√f + 2 log₁₀(ε/3.7D + 2.51/(Re√f)); derive g′(f) analytically and verify it against a central-difference approximation on ten random points. Establish the physical bracket f ∈ [0.008, 0.10] from Moody-chart bounds.

Phase 2 — implementation (weeks 3–4). Implement the four methods with one interface: each takes (Re, eps_over_D), returns f, iterations and error history. Unit checks: a known textbook case (Re = 10⁵, ε/D = 0.0001 → f ≈ 0.018) must reproduce to six decimals.

Phase 3 — experiment (weeks 5–7). Draw 1,000 Latin-hypercube samples in the two parameters. Run every method on every point with tolerance |Δf| < 10⁻¹⁰ and a 100-iteration cap. Store the full error history per run.

Phase 4 — analysis (weeks 8–10). Compute iterations-to-convergence distributions; estimate observed order p from three consecutive errors via the standard log-ratio formula; plot error curves for representative runs; map every failure onto the violated assumption (bracket, seed, contraction).

Phase 5 — write-up (weeks 11–12). Draft the report chapters, regenerate every figure from the notebook, rehearse the viva demonstration, and freeze the final results tables.

The study is a data pipeline in five stages, and its architecture is the flow of numbers from parameter sample to conclusion.

flowchart TD
  A["Latin-hypercube sampler<br/>1000 (Re, eps/D) points"] --> B["Method runners<br/>bisection / Newton / secant / fixed-point"]
  B --> C["Run log<br/>iterations, error history, status"]
  C --> D["Analysis layer<br/>order estimates, failure map, plots"]
  D --> E["Report artefacts<br/>tables, figures, worked example"]
  F["Reference solutions<br/>mpmath 50-digit run"] --> B

A 50-digit-precision reference solution (via the mpmath library with the same bisection on a verified bracket) anchors the error measurements, so every method's error is measured against an exact yardstick rather than against another iterate. The method runners are pure functions with no shared state, which is what makes the whole 4,000-run experiment reproducible from one notebook cell. The analysis layer never modifies run logs; it only reads them, so figures can be restyled without rerunning experiments.

Three natural extensions follow. First, the same protocol can rank the modern explicit approximations (Haaland, Romeo, and successors) for absolute error against the 50-digit reference across the same parameter space — an accuracy-per-arithmetic-operation study with genuine engineering value. Second, the single-equation framework extends to systems: the two-equation pipeline network problem introduces Jacobians and raises the question of whether Newton's advantage grows or shrinks in higher dimensions. Third, the run log invites statistical treatment — bootstrap confidence intervals on median iteration counts — connecting numerical analysis to the inference tools of the data-science electives.