Root finding

Engineering is full of equations that cannot be solved with algebra. The friction factor in a pipe appears on both sides of the Colebrook equation. The depth at which a float settles satisfies a cubic. The natural frequencies of a beam satisfy cos βL cosh βL = −1. The interest rate that makes a project break even satisfies a polynomial of high degree. Each is a root-finding problem: find x such that f(x) = 0.

ENCM-100 introduced bisection as a programming exercise. This chapter treats root finding properly: bracketing and open methods, how fast they converge, when they fail, and how professional library routines combine them.

Setting up the problem

Any equation can be rearranged as f(x) = 0. Before reaching for a method, sketch or plot f: it shows how many roots there are, roughly where, and whether anything odd (a discontinuity, a near-double root) is lurking. A good starting guess, from a plot, a simplified model or engineering judgment, is half the battle.

Bracketing methods

Bisection

If f is continuous and f(a) and f(b) have opposite signs, the intermediate value theorem (MATH-100) guarantees a root in [a, b]. Bisection evaluates the midpoint m = (a + b)/2, keeps the half where the sign changes, and repeats. Each step halves the bracket, so after n steps the error is at most \((b - a)/2^{n+1}\), and the number of steps for a tolerance ε is

Bisection is slow (one binary digit per step, about 3.3 steps per decimal digit) but bulletproof: it cannot diverge. The false position (regula falsi) method replaces the midpoint by the point where the straight line through (a, f(a)) and (b, f(b)) crosses zero. It is often faster, but can stall with one end of the bracket stuck.

Open methods

Open methods start from one or two guesses that need not bracket the root. They are usually much faster, but they can diverge.

Fixed-point iteration

Rearrange f(x) = 0 as x = g(x) and iterate \(x_{k+1} = g(x_k)\). Near the root r, the error is multiplied by about \(g'(r)\) at each step, so the iteration converges if \(|g'(r)| < 1\), and faster the smaller \(|g'(r)|\) is. The same equation can be rearranged in ways that converge or diverge.

Newton's method

Newton's method (MATH-100) follows the tangent line at the current guess to where it crosses zero:

Near a simple root, Newton converges quadratically: the error is roughly squared at each step, \(e_{k+1} \approx C e_k^2\), so the number of correct digits roughly doubles each iteration.

Newton's method for √a, applied to \(x^2 - a = 0\), gives \(x_{k+1} = \frac{1}{2}(x_k + a/x_k)\), the method of the ancient Babylonians. From x₀ = 1, the errors in √2 are 0.086, 0.0025, 2.1 × 10⁻⁶, 1.6 × 10⁻¹² and then machine precision.

When Newton fails

Newton's method needs a derivative and a reasonable starting point. It can fail in several ways:

• A zero or tiny derivative sends the next guess far away (a horizontal tangent never meets the axis).
• From a poor start it can diverge or cycle. For \(f(x) = \tan^{-1}x\), starting from x₀ = 1.3 converges to 0, but from 1.5 the guesses go −1.69, 2.32, −5.11, 32.3, −1575: each overshoot is worse.
• At a multiple root (where f' is also zero, as in \((x - 1)^2\)), convergence slows to linear.

The secant method

When f' is unavailable or expensive, replace it with the slope through the last two points:

The secant method needs two starting values but only one new function evaluation per step. Its order of convergence is the golden ratio, about 1.618: slower than Newton per step, but often faster per unit of work. On the floating ball from 0.02 and 0.05, its errors are 2.2 × 10⁻³, 3.7 × 10⁻⁵, 2.3 × 10⁻⁷, 2.1 × 10⁻¹¹ and then machine precision.

Comparing convergence

If the error satisfies \(|e_{k+1}| \approx C|e_k|^p\), p is the order of convergence. Bisection is linear (p = 1, C = 1/2), fixed-point iteration is linear with C = |g'(r)|, the secant method has p ≈ 1.618 and Newton p = 2.

Professional solvers combine the strengths. Brent's method (scipy.optimize.brentq, MATLAB fzero) keeps a bracket, so it cannot fail, but takes secant or inverse-quadratic steps when they behave well and falls back to bisection when they do not. If you can bracket the root, this is the method to use.

Systems of nonlinear equations

Newton's method extends to n equations in n unknowns, \(\mathbf{F}(\mathbf{x}) = \mathbf{0}\). The derivative becomes the Jacobian matrix J of partial derivatives (MATH-209), and each step solves a linear system (Lesson 3):

This is how circuit simulators solve nonlinear circuits, how structural programs handle large deflections and how process simulators converge flowsheets. Convergence is again quadratic near the solution, and again a good starting guess matters.

from scipy.optimize import brentq, newton, fsolve
import numpy as np

f = lambda x: x**3 - 0.165*x**2 + 3.993e-4
print("ball depth (Brent):", brentq(f, 0, 0.11))
print("ball depth (Newton):", newton(f, 0.05, fprime=lambda x: 3*x**2 - 0.33*x))

def colebrook(u, rel=0.0002, Re=1e5):
    return u + 2*np.log10(rel/3.7 + 2.51*u/Re)
u = brentq(colebrook, 1, 20)
print(f"friction factor f = {1/u**2:.5f}")

F = lambda v: [v[0]**2 + v[1]**2 - 4, np.exp(v[0]) + v[1] - 1]
print("intersection:", fsolve(F, [1, -1]))

Output: ball depth 0.0623775815… m by both methods; friction factor f = 0.01901; intersection [1.00416874 −1.72963729].