Prev
Unit 2 · Nonlinear Equations
Bracketing Methods
60m100 XP
Next
🔍Nonlinear Equations

Bracketing Methods

The intermediate value theorem, bisection, false position, the Illinois fix and Brent's method

London, 1939
Water flowing through a pipe loses pressure to friction. How much depends on one dimensionless number, the Darcy friction factor ff, and every pump, pipeline, cooling loop and fuel feed line is sized from it. In 1939 Cyril Colebrook, working with data from smooth and sand-roughened pipes, published the law that engineers still use for turbulent flow: 1f=−2log⁡10 ⁣(ε/D3.7+2.51Re f)\dfrac{1}{\sqrt f} = -2\log_{10}\!\left(\dfrac{\varepsilon/D}{3.7} + \dfrac{2.51}{Re\,\sqrt f}\right) Look at where ff is. It sits on the left, and again on the right inside the logarithm. No algebra gets it out. There is no formula for ff. Five years later Lewis Moody plotted the solutions as a chart — friction factor against Reynolds number, one curve per roughness — and for decades engineers solved Colebrook's equation by reading a graph with a ruler. The Moody diagram is still printed in every fluid mechanics textbook. A computer cannot read a chart. It needs a procedure that turns an equation it cannot rearrange into a number — and ideally a guarantee that the procedure will work.
Every equation can be written as g(x)=0g(x) = 0: move everything to one side. A value x∗x^* with g(x∗)=0g(x^*) = 0 is a root. For Colebrook,
g(f)=1f+2log⁡10 ⁣(ε/D3.7+2.51Ref)g(f) = \dfrac{1}{\sqrt f} + 2\log_{10}\!\left(\dfrac{\varepsilon/D}{3.7} + \dfrac{2.51}{Re\sqrt f}\right)
and the friction factor is its root. The whole of this unit is about finding roots of functions we can evaluate but not invert.
💡The intermediate value theorem
If gg is continuous on [a,b][a, b] and g(a)g(a) and g(b)g(b) have opposite signs, then gg has at least one root in (a,b)(a, b). An interval with a sign change is called a bracket. A bracket is a certificate: as long as you keep one, a root is trapped inside it, whatever else happens. The methods in this lesson never let go of their bracket — that is where their guarantee comes from.
Bisection is the simplest way to shrink a bracket. Evaluate gg at the midpoint m=(a+b)/2m = (a+b)/2. If g(a)g(a) and g(m)g(m) differ in sign, the root is in [a,m][a, m]; otherwise it is in [m,b][m, b]. Keep that half and repeat.
Each step halves the bracket, whatever gg looks like. After nn steps the width is (b−a)/2n(b-a)/2^n, so reaching a tolerance τ\tau takes
n  ≥  log⁡2 ⁣b−aτn \;\ge\; \log_2\!\frac{b-a}{\tau}

Bisection's step count is known before you start. Each step buys exactly one binary digit, so about 3.3 steps per decimal digit.

Derivation
Counting bisection steps
Each step keeps exactly half of the bracket, whatever gg is. Starting from width b−ab - a, after nn steps the width is (b−a)/2n(b - a)/2^n. We stop once it is below the tolerance τ\tau:
b−a2n≤τ⟺2n≥b−aτ⟺n≥log⁡2b−aτ.\frac{b - a}{2^n} \le \tau \quad\Longleftrightarrow\quad 2^n \ge \frac{b - a}{\tau} \quad\Longleftrightarrow\quad n \ge \log_2\frac{b - a}{\tau}.
One more decimal digit means dividing τ\tau by 10, which adds log⁡210=3.32\log_2 10 = 3.32 steps. For Colebrook, log⁡2(0.092/10−10)=29.8\log_2(0.092/10^{-10}) = 29.8; the lab's loop stops on the half-width, one halving earlier, after 29 steps.
📐Colebrook by bisection
Take Re=105Re = 10^5 and ε/D=10−4\varepsilon/D = 10^{-4}. Friction factors for turbulent flow lie between about 0.008 and 0.1, and g(0.008)>0g(0.008) > 0, g(0.1)<0g(0.1) < 0: a bracket. For a tolerance of 10−1010^{-10} the formula predicts log⁡2(0.092/10−10)≈29.8\log_2(0.092/10^{-10}) \approx 29.8 halvings; the loop that stops when the half-width drops below 10−1010^{-10} takes 29. The result is f=0.0185138660f = 0.0185138660. For comparison, Haaland's explicit formula of 1983, designed precisely so that nobody has to iterate, gives 0.018265 — 1.3% low. Good enough to size a pump; not good enough to be the reference that a pipe-network solver converges to.
When do you stop? Three tests are common, and each can fool you on its own.
  • Bracket width ∣b−a∣<τ|b - a| < \tau — the honest test for bisection, but use a relative width τ(∣a∣+∣b∣)\tau(|a|+|b|) when the root may be large, or you will ask for more digits than a double holds.
  • Residual ∣g(x)∣<τ|g(x)| < \tau — meaningless without units: gg can be scaled by any constant. A flat function has a tiny residual far from its root; a steep one has a large residual right next to it.
  • Iteration cap — always have one. A bug or a function that is not continuous should end in an error message, not an infinite loop.
when you can measure what you are speaking about, and express it in numbers, you know something about it
– William Thomson, Lord Kelvin, lecture on electrical units of measurement (1883)
False position (regula falsi) uses more information. Instead of the midpoint, it takes the point where the straight line through (a,g(a))(a, g(a)) and (b,g(b))(b, g(b)) crosses zero:
x=a g(b)−b g(a)g(b)−g(a)x = \dfrac{a\,g(b) - b\,g(a)}{g(b) - g(a)}
Near a root a smooth function is almost a straight line, so this should be much closer than the midpoint. Often it is. But on a function that is convex across the bracket, one end of the bracket never moves: every new point lands on the same side. The bracket does not shrink to zero, and convergence slows to a crawl. On x10−1x^{10} - 1 over [0,1.3][0, 1.3], false position takes over 100 steps to do what bisection does in 34.
Derivation
Where the false-position point comes from
The straight line through (a,g(a))(a, g(a)) and (b,g(b))(b, g(b)) has slope (g(b)−g(a))/(b−a)\big(g(b) - g(a)\big)/(b - a), so
y(x)=g(a)+g(b)−g(a)b−a (x−a).y(x) = g(a) + \frac{g(b) - g(a)}{b - a}\,(x - a).
Set y=0y = 0 and solve for xx:
x=a−g(a) b−ag(b)−g(a)=a(g(b)−g(a))−g(a)(b−a)g(b)−g(a)=a g(b)−b g(a)g(b)−g(a).x = a - g(a)\,\frac{b - a}{g(b) - g(a)} = \frac{a\big(g(b) - g(a)\big) - g(a)(b - a)}{g(b) - g(a)} = \frac{a\,g(b) - b\,g(a)}{g(b) - g(a)}.
Because g(a)g(a) and g(b)g(b) have opposite signs, the denominator is never zero and xx always lies strictly between aa and bb — the bracket is never lost. The secant method of the next lesson uses exactly this formula, but with the two most recent points instead of the two ends of a bracket.
💡The Illinois fix
When the same end of the bracket survives twice in a row, halve the function value stored at that end. The next secant line then lands on the other side of the root, and the stuck end finally moves. The change is two lines of code. On x10−1x^{10} - 1 it cuts the work from over 100 evaluations to 14, and the method keeps its bracket — so it keeps the guarantee. Dowell and Jarratt analysed it in 1971 and showed its order of convergence is about 1.44.
Brent's method (Richard Brent, 1973, building on Theodorus Dekker's 1969 algorithm) is the industrial answer. It keeps a bracket at all times, but at each step it tries a fast step — a secant step or inverse quadratic interpolation through the last three points — and accepts it only if it lands inside the bracket and shrinks it fast enough. Otherwise it falls back to bisection.
The result converges superlinearly on smooth functions and is never much slower than bisection on bad ones. It is scipy.optimize.brentq, MATLAB's fzero, and the root finder inside most engineering codes. When you have a bracket, it is the default.
⚠What a bracket cannot do
A sign change finds an odd number of roots; it says nothing about which. Double roots — x2x^2 at 0 — have no sign change at all and are invisible to bracketing. And a function with a pole, like tan⁡x\tan x across π/2\pi/2, changes sign without a root: bisection will happily converge to the pole. Plot first, or check that ∣g∣|g| is small at the answer.
💡Where you will use this
Hydraulic Engineering: Colebrook inside every pipe-network solver. Geotechnical Engineering: the critical depth of a slip surface and the factor of safety at which a slope is just stable. Rocket Propulsion: the exit Mach number of a nozzle from its area ratio, where you must choose the subsonic or supersonic branch — a bracket chooses it for you. Orbital Mechanics: the time of flight in Lambert's problem as a monotonic function of one parameter.
📚Further reading
  • Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 5, 'Nonlinear Equations' — bisection, false position and safeguarded methods.
  • Richard Brent, Algorithms for Minimization without Derivatives (Prentice-Hall, 1973; Dover reprint 2002), chapter 4 — Brent's own account of his zero finder, with the ALGOL code it was published in.
  • Cyril Colebrook, 'Turbulent flow in pipes, with particular reference to the transition region between the smooth and rough pipe laws', Journal of the Institution of Civil Engineers 11, 133–156 (1939).
↓But there is a problem
Bisection is slow in a very precise sense: 29 steps for ten digits, 52 steps for all sixteen, whatever the function. It never uses the most obvious piece of information — which way the curve is heading. In 1669 Isaac Newton wrote down a different idea. Replace the curve by its tangent line, find where the line hits zero, and repeat. His first example was a cubic. From a guess of 2 it gave the root to 10 digits in three steps.