Root-Finding and Iterative Methods

A startling amount of physics ends at the same wall: an equation $f(x)=0$ that no rearrangement will solve in closed form. Dispersion relations, the transcendental matching conditions of a square well, Kepler's equation for the orbit of a planet, and the Bohr–Sommerfeld quantization condition that opens the door to the WKB Approximation all deliver a curve and ask where it crosses zero — and algebra simply refuses to answer. Root-finding is the craft of answering anyway, by replacing the impossible demand for an exact inverse with a recursion $x_{n+1}=g(x_n)$ engineered so that its limit, if it exists, is the root. Two questions decide everything. Does the iteration converge, and how fast? Convergence is a statement about whether the map $g$ pulls nearby points together — the contraction-mapping principle of Banach, one sentence of analysis that governs every scheme here. Speed is the order of convergence $p$ , defined by $|e_{n+1}|\le C|e_n|^{p}$ with $e_n=x_n-x_*$ : bisection crawls linearly ( $p=1$ ), the secant method sprints superlinearly ( $p\approx1.618$ ), and Newton's method doubles its correct digits every step ( $p=2$ ). The engineer's answer is never one method but a marriage of two. Bracketing schemes trap a root inside a known sign change and can never lose it, so they always converge, slowly and safely; open schemes fly toward a root once near one but can overshoot, cycle, or diverge from a poor start. Production solvers — Brent's method, the workhorse behind every `scipy.optimize.brentq` — braid bisection's guarantee together with interpolation's speed. The same fixed-point idea, lifted to vectors and function spaces, becomes the Jacobian Newton step, Picard iteration, and the self-consistency loops of Hartree–Fock and the renormalization group.

Key takeaways

Root-finding solves $f(x)=0$ for nonlinear $f$ by generating a sequence $x_n$ that converges to a root $x_*$ ; every method is a fixed-point map $x_{n+1}=g(x_n)$ .
Bisection brackets a sign change and halves the interval each step — guaranteed but only linear, $|e_{n+1}|\approx\tfrac12|e_n|$ , so about $0.30$ decimal digit per step.
Newton–Raphson linearizes $f$ at $x_n$ , giving $x_{n+1}=x_n-f(x_n)/f'(x_n)$ with quadratic convergence at a simple root: correct digits double each iteration.
Fixed-point iteration converges iff the map contracts, $|g'(x_*)|<1$ , at the linear rate $|g'(x_*)|$ ; Newton is the special map engineered so $g'(x_*)=0$ .
The secant method replaces $f'$ by a finite difference and converges at order $\varphi=(1+\sqrt5)/2\approx1.618$ without any derivative.
Speed and safety are independent qualities — real codes hybridize a bracketer with an open method (Brent) to get both.
Key idea: every root-finder is a fixed-point map $x_{n+1}=g(x_n)$ , and its convergence and order are dictated entirely by the leading non-vanishing derivative of $g$ at the fixed point $x_*$ .

The core insight

Every root-finder is a fixed-point map $x_{n+1}=g(x_n)$ ; whether it converges, and how fast, is fixed entirely by how $g$ behaves near the fixed point — specifically by the first non-vanishing derivative of $g$ at $x_*$ .

Where this lives today

Root-finding is invisible infrastructure: it runs inside almost every quantitative calculation without ever being named. When a physicist calls `scipy.optimize.brentq` to find a bound-state energy, they are running Brent's 1973 hybrid — bisection braided with inverse-quadratic interpolation — which guarantees convergence while usually achieving superlinear speed. In electronic-structure codes the self-consistent field loop is nothing but high-dimensional fixed-point iteration: the map takes a trial density to the density it produces, and the notorious "charge sloshing" instability is exactly the failure of the contraction condition $|g'|<1$ , tamed by density mixing (damped iteration). Spacecraft navigation at JPL solves the same Kepler equation Newton's contemporaries did, only billions of times per mission. And at the frontier, machine-learning training is fixed-point iteration in a space of billions of parameters, where the analogue of $|g'|<1$ is precisely what decides whether a run converges or diverges.

The physical picture

Think of $f(x)=0$ as the question "where does this curve cross the axis?" Newton's method answers it with a ruler: stand at your current guess, draw the tangent line, and slide down it to the axis — that landing point is your next guess. Where the curve is nearly straight the tangent is an almost perfect stand-in, so the new guess is dramatically better than the old one. Bisection answers the same question with a single yes/no probe: it keeps a bracket that is known to straddle the root and asks only "is the root in the left half or the right half?", throwing away the wrong half each time. The first figure shows Newton's tangent construction on a concrete curve, with real numbers attached.

Figure

Newton's tangent construction on f(x) = x² − 2 near √2: the tangent at x₀ = 1.5 with slope f′ = 3 strikes the axis at x₁ = 1.41667, and one more step lands on 1.41422 — six correct digits against the true root √2 = 1.41421356

Fixed-point iteration is the same idea stripped to its skeleton. Rewrite $f(x)=0$ as $x=g(x)$ and feed each output back as the next input; geometrically you bounce between the curve $y=g(x)$ and the diagonal $y=x$ . If $g$ is shallower than the diagonal near the crossing ( $|g'|<1$ ) the bounces spiral inward and you converge; if it is steeper they fly outward and you diverge. The whole theory lives in that one slope. And the payoff of a well-chosen map is dramatic: quadratic convergence means the number of correct digits roughly doubles each step — $3\to6\to12$ in three iterations — whereas linear convergence adds a fixed number of bits per step, so bisection needs some fifty iterations to reach machine precision. The second figure makes that gulf concrete by racing the three methods on one root.

Figure

Error versus iteration for the root √2: bisection loses a fixed 0.30 digit per step (linear), the secant climbs superlinearly at order 1.618, and Newton doubles its correct digits each step, the errors falling 8.6×10⁻² then 2.5×10⁻³ then 2.1×10⁻⁶ then 1.6×10⁻¹²

Order of convergence is a strictly local statement, and that is the recurring trap. A method that converges quadratically once it is near the root can still be sent diverging by a bad starting guess; speed near the root and reliability far from it are independent virtues. That single fact — that you can be fast or safe but rarely both from one formula — is why every serious library ships a hybrid.

The core equations

The object of study is the root and the error that measures our distance from it. Fixing this notation once makes every convergence statement below a claim about a single sequence of numbers $e_n$ and how fast it collapses to zero.

f(x_*)=0,\qquad e_n\equiv x_n-x_* \qquad (1)

where: $f$ — continuous (often smooth) real function whose zero is sought; $x_*$ — the exact root; $x_n$ — the $n$ -th iterate; $e_n$ — the error of the $n$ -th iterate. Every method below is judged solely by how $|e_n|$ decays. The safest scheme needs almost nothing: only that $f$ is continuous and changes sign across a bracket. By the Intermediate Value Theorem a sign change guarantees a root inside, and halving the interval halves the guaranteed error without ever consulting the shape of $f$ .

\operatorname{sign} f(a_n)\ne\operatorname{sign} f(b_n)\ \Rightarrow\ c_n=\frac{a_n+b_n}{2} \qquad (2)

where: $[a_n,b_n]$ — the current bracket, guaranteed to contain a root; $c_n$ — its midpoint, which replaces whichever endpoint keeps the sign change straddling the root. The error obeys $|e_n|\le(b_0-a_0)/2^{\,n}$ deterministically. The fastest common scheme trades that guarantee for speed. Newton–Raphson models $f$ by its tangent line at the current iterate and jumps to where the tangent hits the axis, which is exact for a straight line and superb for a gently curving one.

x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)} \qquad (3)

where: $f'(x_n)$ — the derivative at the current iterate, which must be non-zero. Reading it: the correction $-f/f'$ is large where the function is far from zero or the slope is shallow, and small where you are close or the slope is steep. The three-line construction is this: linearize, $f(x)\approx f(x_n)+f'(x_n)(x-x_n)$ ; set the model to zero, $0=f(x_n)+f'(x_n)(x_{n+1}-x_n)$ ; and solve for $x_{n+1}$ . Every method on this page is secretly the same object — a map whose fixed point is the root. Writing $f(x)=0$ as $x=g(x)$ exposes the common skeleton and lets one theorem govern them all.

x_{n+1}=g(x_n),\qquad x_*=g(x_*) \qquad (4)

where: $g$ — the iteration map, chosen so that $f(x)=0\Leftrightarrow x=g(x)$ ; its fixed point is the sought root. Newton is the choice $g(x)=x-f(x)/f'(x)$ ; bisection is a set-valued version of the same idea. Speed is quantified by a single exponent that compares consecutive errors. This one definition sorts every method into a class and tells you how many iterations a target precision will cost.

|e_{n+1}|\le C\,|e_n|^{p} \qquad (5)

where: $p$ — the order ( $1$ linear, $2$ quadratic; $1<p<2$ superlinear); $C$ — the asymptotic error constant, which for $p=1$ must satisfy $C<1$ for convergence at all. Reading it: with $p=2$ the number of correct digits doubles per step; with $p=1$ it grows by a fixed amount $-\log_{10}C$ per step. When the derivative is expensive or unavailable, the secant method estimates it from the two most recent points, buying most of Newton's speed for one function evaluation per step and no calculus.

x_{n+1}=x_n-f(x_n)\,\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})} \qquad (6)

where: $x_{n-1},x_n$ — the two previous iterates, through which the secant line is drawn as a stand-in for the tangent. Its order is the golden ratio $\varphi=(1+\sqrt5)/2\approx1.618$ , derived below. The same fixed-point machinery lifts verbatim to systems of equations, where the scalar slope becomes the Jacobian matrix and division becomes a linear solve.

\mathbf{x}_{n+1}=\mathbf{x}_n-J(\mathbf{x}_n)^{-1}\mathbf{f}(\mathbf{x}_n) \qquad (7)

where: $\mathbf{f}$ — the vector of nonlinear equations $\mathbf{f}(\mathbf{x})=\mathbf{0}$ ; $J$ — the Jacobian, $J_{ij}=\partial f_i/\partial x_j$ . In practice one never inverts $J$ but solves the linear system $J\,\delta=-\mathbf{f}$ for the step $\delta$ . Underneath all convergence proofs sits one condition: the map must pull points together. This inequality is the entire content of the Banach fixed-point theorem and the master switch for whether any iteration works.

|g(x)-g(y)|\le L\,|x-y|,\quad L<1 \qquad (8)

where: $L$ — the Lipschitz (contraction) constant; $L<1$ guarantees a unique fixed point and convergence from any start in the interval. For a differentiable $g$ the local rate is $|g'(x_*)|$ , and $L\ge|g'(x_*)|$ .

Derivation AQuadratic convergence of Newton's method

The goal is to prove that Newton's method squares the error each step at a simple root, and to pin down the exact constant that multiplies $e_n^2$ . We assume a simple root — $f(x_*)=0$ with $f'(x_*)\ne0$ — and $f\in C^2$ near it, then expand $f$ to second order about the current iterate $x_n$ , insert the Newton step, and watch the first-order error terms cancel identically, leaving a clean $e_n^2$ law whose constant is the curvature-to-slope ratio at the root. This single cancellation is the reason Newton is the gold standard, and its failure is why multiple roots misbehave.

Step 1. Taylor-expand $f$ about $x_n$ and evaluate at the exact root, keeping the exact remainder: $0=f(x_*)=f(x_n)+f'(x_n)(x_*-x_n)+\tfrac12 f''(\xi_n)(x_*-x_n)^2$ , with $\xi_n$ between $x_n$ and $x_*$ , and note $x_*-x_n=-e_n$ . Why: evaluating the model at the point where $f$ genuinely vanishes is what lets the unknown root enter the algebra exactly rather than approximately.

Step 2. Divide through by $f'(x_n)$ and use $x_*-x_n=-e_n$ to obtain $0=\dfrac{f(x_n)}{f'(x_n)}-e_n+\dfrac{f''(\xi_n)}{2f'(x_n)}e_n^{2}$ . Why: dividing by the slope manufactures precisely the Newton correction $f/f'$ as a term, so the next step can substitute it away.

Step 3. Insert the Newton step. The next error is $e_{n+1}=x_{n+1}-x_*=\big(x_n-\tfrac{f(x_n)}{f'(x_n)}\big)-x_*$ , and $x_n-x_*=e_n$ , so $e_{n+1}=e_n-\dfrac{f(x_n)}{f'(x_n)}$ . Why: this rewrites the abstract "next error" purely in terms of quantities that Step 2 already relates, closing the loop.

Step 4. Combine Steps 2 and 3: the line from Step 2 says $\dfrac{f(x_n)}{f'(x_n)}=e_n-\dfrac{f''(\xi_n)}{2f'(x_n)}e_n^{2}$ , and substituting gives $e_{n+1}=\dfrac{f''(\xi_n)}{2f'(x_n)}\,e_n^{2}$ . Why: the $e_n$ terms cancel exactly — that cancellation, not any approximation, is the origin of quadratic order.

Step 5. Take the limit $x_n\to x_*$ , so $\xi_n\to x_*$ , and the boxed error law follows with order $p=2$ ; the constant diverges precisely when $f'(x_*)=0$ , i.e. at a multiple root. Why: the constant's blow-up is a warning label — quadratic convergence is a promise made only at simple roots.

\lim_{n\to\infty}\frac{e_{n+1}}{e_n^{2}}=\frac{f''(x_*)}{2f'(x_*)}

Verified: matches Burden & Faires Thm 2.9 and the Newton convergence result in Arfken–Weber; reproduced numerically for $f=x^2-2$ , where errors fall $8.6\times10^{-2},\,2.5\times10^{-3},\,2.1\times10^{-6},\,1.6\times10^{-12}$ .

Go deeper · Derivation B · 2 minConvergence of fixed-point iteration by the Banach contraction principle

Derivation BConvergence of fixed-point iteration by the Banach contraction principle

The goal is to show that a map which contracts distances has exactly one fixed point and that simple iteration reaches it from any starting point in the interval, at a linear rate set by the contraction constant. We suppose $g:[a,b]\to[a,b]$ is continuous and Lipschitz with constant $L<1$ , then establish existence from the Intermediate Value Theorem, uniqueness from the contraction inequality, and geometric decay of the error by iterating that same inequality — the complete Banach argument specialized to the real line. The final step shows how engineering $g'(x_*)=0$ upgrades this linear result to the quadratic one Newton enjoys.

Step 1. Existence: since $g$ maps $[a,b]$ into itself, the continuous function $h(x)=g(x)-x$ satisfies $h(a)=g(a)-a\ge0$ and $h(b)=g(b)-b\le0$ , so by the Intermediate Value Theorem $h$ has a zero — a fixed point $x_*$ exists. Why: mapping the interval into itself is what forces the graph of $g$ to cross the diagonal, which is geometrically what a fixed point is.

Step 2. Uniqueness: if $x_*$ and $y_*$ were both fixed points, then $|x_*-y_*|=|g(x_*)-g(y_*)|\le L|x_*-y_*|$ , and with $L<1$ this can only hold if $|x_*-y_*|=0$ . Why: a genuine contraction cannot keep two distinct points at their distance, so it can pin down only one.

Step 3. Error recursion: apply the contraction bound to the iterate and the root, $|e_{n+1}|=|g(x_n)-g(x_*)|\le L|x_n-x_*|=L|e_n|$ . Why: this is the workhorse inequality — one application per step, and it says the error shrinks by at least the factor $L$ every iteration.

Step 4. Iterate the bound to get $|e_n|\le L^{n}|e_0|\to0$ : convergence is linear with asymptotic rate $|g'(x_*)|$ , since for differentiable $g$ the sharpest local $L$ is $|g'(x_*)|$ . Why: a geometric sequence with ratio below one collapses to zero, and its ratio is exactly the local slope of the map.

Step 5. Higher order: if in addition $g'(x_*)=0$ , a Taylor expansion gives $e_{n+1}\approx\tfrac12 g''(x_*)\,e_n^{2}$ — the iteration is now quadratic, and Newton's map $g=x-f/f'$ is precisely the one built to have $g'(x_*)=0$ . Why: this unifies the page — Newton is not a separate miracle but the contraction whose first derivative was deliberately annihilated.

|e_{n+1}|=|g(x_n)-g(x_*)|\le L\,|e_n|\ \Rightarrow\ |e_n|\le L^{n}|e_0|\to0

Verified: the standard Banach fixed-point theorem, Stone–Goldbart on contraction mappings; the rate $|g'(x_*)|$ is confirmed numerically for $g(x)=\cos x$ , whose fixed point $x_*\approx0.73909$ has $g'(x_*)=-\sin x_*\approx-0.674$ .

Go deeper · Derivation C · 1 minThe golden-ratio order of the secant method

Derivation CThe golden-ratio order of the secant method

The goal is to explain why replacing Newton's exact derivative by a two-point finite difference costs so little — the order drops only from $2$ to the golden ratio $\varphi\approx1.618$ — and to derive that number from the secant error recurrence. We first establish that the secant error obeys a product law $e_{n+1}\approx C\,e_n e_{n-1}$ with the same Newton constant $C=f''(x_*)/2f'(x_*)$ , then impose the self-similar ansatz $|e_n|\sim A|e_{n-1}|^{p}$ that any power-law convergence must satisfy, and match exponents to land on $p=1+1/p$ , whose positive root is the golden ratio — which is why the secant is the efficient derivative-free cousin of Newton.

Step 1. Expanding $f$ to second order at $x_n$ and $x_{n-1}$ and substituting into the secant formula (6) gives the error recurrence $e_{n+1}\approx\dfrac{f''(x_*)}{2f'(x_*)}\,e_n\,e_{n-1}\equiv C\,e_n e_{n-1}$ . Why: the secant line's deviation from the tangent is second order, so the leading error is bilinear in the two most recent errors rather than quadratic in one.

Step 2. Posit power-law convergence, $|e_{n+1}|=A|e_n|^{p}$ , which also means $|e_n|=A|e_{n-1}|^{p}$ , hence $|e_{n-1}|=(|e_n|/A)^{1/p}$ . Why: if the method has a definite order $p$ at all, consecutive errors must be related by that same power — this is the only self-consistent asymptotic form.

Step 3. Substitute both into Step 1: $A|e_n|^{p}\approx |C|\,|e_n|\,(|e_n|/A)^{1/p}=|C|A^{-1/p}\,|e_n|^{\,1+1/p}$ . Matching the powers of $|e_n|$ forces $p=1+\dfrac1p$ . Why: for the relation to hold as $|e_n|\to0$ the exponents on both sides must agree — the constants only fix $A$ , not the order.

Step 4. Solve $p^2-p-1=0$ , whose positive root is $p=\dfrac{1+\sqrt5}{2}=1.61803\ldots$ — superlinear, achieved with a single new function evaluation and no derivative. Why: this is the quantitative payoff — the secant is nearly as fast as Newton yet needs neither $f'$ nor a second sample of $f$ per step, so when $f'$ is costly it is often the better buy.

p=1+\frac1p\ \Longrightarrow\ p=\frac{1+\sqrt5}{2}=\varphi\approx1.618

Verified: the recurrence $e_{n+1}=C\,e_n e_{n-1}$ and order $\varphi$ appear in Numerical Recipes §9.2 and Burden & Faires; reproduced for $f=x^2-2$ , where the secant reaches $2\times10^{-16}$ in seven steps versus Newton's four.

Limiting cases

How each scheme behaves is dictated by the local structure of $f$ and of the iteration map $g$ at the root:

Case	Result	Meaning
Simple root, good start	$\vert e_{n+1}\vert=\dfrac{\vert f''(x_)\vert}{2\vert f'(x_)\vert}\,e_n^{2}$	Newton converges quadratically; the error constant is the curvature-to-slope ratio at the root.
Multiple root, multiplicity $m$	$e_{n+1}\approx\Big(1-\dfrac1m\Big)e_n$	Newton degrades to linear at rate $1-1/m$ (so $1/2$ for a double root); modified Newton $x_{n+1}=x_n-m\,f/f'$ restores quadratic order.
Fixed point with $g'(x_*)=0$	$\vert e_{n+1}\vert\approx\tfrac12\vert g''(x_*)\vert\,e_n^{2}$	Any map whose first derivative vanishes is automatically second order — exactly what Newton engineers.
Contraction $0<\vert g'(x_*)\vert<1$	$\vert e_{n+1}\vert\approx\vert g'(x_*)\vert\,\vert e_n\vert$	Plain fixed-point iteration: linear, gaining $-\log_{10}\vert g'\vert$ digit per step.
Bisection	$\vert e_n\vert\le\dfrac{b_0-a_0}{2^{\,n}}$	Error halves deterministically whatever the shape of $f$ ; global but only linear at rate $1/2$ ( $\approx0.30$ digit per step).

Common trapNewton's quadratic convergence is a local promise, not a global guarantee

It is tempting to read "Newton converges quadratically" as "Newton always converges fast." Both halves can fail. The quadratic law holds only near a simple root: at a multiple root $f'(x_*)=0$ , the error constant $f''/2f'$ diverges and convergence collapses to merely linear, $e_{n+1}\approx(1-1/m)e_n$ . And even at a perfectly simple root, a poor starting guess can send Newton overshooting past the root, cycling between two points forever, or diverging to infinity when it wanders onto a region where $f'$ is small and the tangent nearly horizontal. The lesson every numerical library encodes: never trust an open method naked. Wrap it in a bracket, check that each step reduces $|f|$ , and fall back to bisection whenever the Newton step tries to leave the bracket.

Worked examples

1 · Newton for √2, digit by digit (confidence builder)

Solve $f(x)=x^2-2=0$ by Newton's method from $x_0=1.5$ , and watch the correct digits double.

Step 1. Form the iteration. With $f=x^2-2$ and $f'=2x$ , the Newton step is $x_{n+1}=x_n-\dfrac{x_n^2-2}{2x_n}=\dfrac12\Big(x_n+\dfrac{2}{x_n}\Big)$ — the ancient Babylonian square-root rule, which is Newton in disguise.

Step 2. Iterate from $x_0=1.5$ . One gets $x_1=\tfrac12(1.5+2/1.5)=\mathbf{1.4166667}$ , then $x_2=\mathbf{1.4142157}$ , then $x_3=\mathbf{1.41421356237}$ . Against $\sqrt2=1.41421356237$ the errors are $e_0=8.6\times10^{-2}$ , $e_1=2.5\times10^{-3}$ , $e_2=2.1\times10^{-6}$ , $e_3=1.6\times10^{-12}$ .

Step 3. Read the pattern. The exponents of the error roughly double each step ( $-2,-3,-6,-12$ ), the signature of quadratic order. The predicted constant is $\dfrac{f''(x_*)}{2f'(x_*)}=\dfrac{2}{2\cdot2\sqrt2}=\dfrac{1}{2\sqrt2}=0.3536$ ; the check $e_2/e_1^2=2.1\times10^{-6}/(2.5\times10^{-3})^2\approx0.35$ confirms it to two figures. Three iterations take a two-digit guess to full double precision.

2 · Kepler's equation for a planet's position

A body orbits with eccentricity $\epsilon=0.2$ . At a mean anomaly $M=0.5$ rad, find the eccentric anomaly $E$ from Kepler's equation $M=E-\epsilon\sin E$ , which has no closed-form inverse.

Step 1. Cast it as a root problem. Define $f(E)=E-\epsilon\sin E-M$ ; then $f'(E)=1-\epsilon\cos E$ , which is bounded between $0.8$ and $1.2$ here, so the tangent never goes flat and Newton is safe. The iteration is $E_{n+1}=E_n-\dfrac{E_n-\epsilon\sin E_n-M}{1-\epsilon\cos E_n}$ .

Step 2. Start from $E_0=M=0.5$ (an excellent seed for modest eccentricity). The residual $f(0.5)=0.5-0.2\sin0.5-0.5=-0.0959$ , giving $E_1=\mathbf{0.6162972}$ . The next residual is $6.9\times10^{-4}$ , giving $E_2=\mathbf{0.6154682}$ ; then $E_3=\mathbf{0.61546817}$ with residual $4\times10^{-8}$ .

Step 3. Confirm and interpret. One more step pins $E=\mathbf{0.615468169}$ rad to machine precision — four iterations total. The residual fell $9.6\times10^{-2}\to6.9\times10^{-4}\to4.0\times10^{-8}\to10^{-16}$ , quadratic as promised. This is the workhorse of orbital mechanics: every ephemeris and spacecraft-tracking routine solves Kepler's equation this way. E = 0.615468 rad, reached in four Newton steps.

Physics of the portrait

Root-finding carries one name above all others — the method that bears it was the first to make correct digits explode, and its inventor built it from the calculus he had just created.

Origins · 1 minIsaac Newton and how this was discovered

Isaac Newton

1643–1727 · the calculus of successive approximation

Isaac Newton described his root-finding scheme around 1669 in De analysi and again in Method of Fluxions, but not in the form students learn today: he applied it to a specific cubic, $x^3-2x-5=0$ , generating corrections term by term without ever writing down a derivative or a general tangent-line formula. Joseph Raphson simplified and generalized it in 1690 into the iterative shape $x_{n+1}=x_n-f(x_n)/f'(x_n)$ that now carries both their names, and Thomas Simpson in 1740 was the first to cast it explicitly in the language of fluxions — derivatives — and to extend it to systems of two equations, making it the general method we use. The convergence theory came much later: the sharp statement that a simple root gives quadratic order, with the error constant $f''/2f'$ , is a nineteenth- and twentieth-century refinement, and the contraction-mapping principle that unifies all of it was proved by Stefan Banach in his 1922 thesis. The threads meet in the modern picture — a root-finder is a fixed-point map, and its order is read off from the derivatives of that map at the root — a synthesis that turns Newton's clever cubic trick into one line of a universal theory that now runs inside every scientific library on Earth.

History

1669

— Newton solves a cubic by successive correction in De analysi

1690

— Raphson recasts the scheme in its modern iterative form

1740

— Simpson expresses it via fluxions and extends it to systems

1879

— the secant and regula falsi methods are analyzed as derivative-free variants

1922

— Banach proves the contraction-mapping fixed-point theorem

1973

— Brent publishes his hybrid bracketing algorithm, still the library default

Going further

Texts. Stone & Goldbart, Mathematics for Physics (CUP, 2009), on contraction mappings · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed. · Press, Teukolsky, Vetterling & Flannery, Numerical Recipes, 3rd ed. (CUP, 2007), Ch. 9 · Burden & Faires, Numerical Analysis, 9th ed., Ch. 2 · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, 3rd ed.

Provenance. The Newton error law $e_{n+1}\to\big(f''/2f'\big)e_n^2$ and the contraction theorem follow standard textbook derivations (Burden & Faires 2.9; Stone–Goldbart). The secant order $p=\varphi=1.6180\ldots$ follows from $p=1+1/p$ ; the bisection bound $2^{-n}$ and the $20$ -iteration count $\lceil6/\log_{10}2\rceil$ are exact. All worked numbers were recomputed here: Newton on $x^2-2$ (errors $8.6\times10^{-2},2.5\times10^{-3},2.1\times10^{-6},1.6\times10^{-12}$ ), the secant reaching $2\times10^{-16}$ in seven steps, Kepler $\epsilon=0.2,M=0.5\Rightarrow E=0.615468169$ , and the Dottie number $\cos x=x\Rightarrow x_*=0.73909$ with rate $0.674$ . A verifier should re-run the Taylor cancellation, the golden-ratio recurrence, and the multiplicity- $m$ rate $1-1/m$ against Arfken–Weber and Riley–Hobson–Bence.

Prerequisites. Taylor's Theorem · Continuity and the Intermediate Value Theorem · Order Symbols and Gauge Functions. Leads to. Bohr–Sommerfeld Quantization · Connection Formulae and Turning Points · WKB Approximation.

Ask the tutor: Give me any equation you cannot solve by hand — a dispersion relation, Kepler's equation, a quantization condition — and I will walk you through choosing a method, bracketing the root, and watching the digits converge.

When should I reach for bisection instead of Newton?

Whenever reliability matters more than speed, or you have a bracket but no good starting guess. Bisection needs only a sign change and continuity, and it cannot fail: the error is guaranteed to halve every step regardless of how badly behaved $f$ is between the endpoints. Newton is far faster near a simple root but needs the derivative, a decent initial guess, and a well-behaved $f'$ ; give it a bad start and it can diverge. The honest answer in practice is "neither alone" — use a hybrid like Brent's method, which keeps a bracket at all times (so it inherits bisection's guarantee) but takes fast interpolation steps whenever they stay inside the bracket. That is why library routines default to Brent rather than to raw Newton.

Why does Newton's method sometimes diverge or cycle?

Because the tangent-line model is only trustworthy near the root. Far away, three things go wrong. If $f'(x_n)$ is small, the tangent is nearly horizontal and its intercept flies off to a distant point — a single flat spot can throw you across the whole domain. If the curvature is strong, the step can overshoot the root and land on the far side, worse than before. And for some functions the iteration settles into a stable cycle, bouncing between two points forever without converging — the classic example is $f(x)=x^3-2x+2$ started at $x_0=0$ , which oscillates between $0$ and $1$ . The cure is always the same: bracket the root, verify each step actually decreases $|f|$ , and reject any Newton step that tries to leave the bracket, substituting a bisection step instead.

What does "order of convergence" actually mean operationally?

It is the exponent $p$ in $|e_{n+1}|\approx C|e_n|^{p}$ , and the cleanest way to feel it is in digits. Linear convergence ( $p=1$ ) adds a fixed number of correct digits per step — bisection buys about $0.30$ of a decimal digit each halving, so roughly ten steps per three digits. Superlinear ( $1<p<2$ ) accelerates: the secant's $p=1.618$ means each step's digit-gain itself grows. Quadratic ( $p=2$ ) doubles the number of correct digits every step: three correct digits become six, then twelve, then twenty-four. This is why Newton reaches double precision in a handful of iterations while bisection needs fifty from the same bracket. The order is a property of the method and the local behavior of $f$ , not of the particular equation.

Is the secant method just a worse Newton?

Not worse — different trade-offs. The secant method sacrifices a little order (1.618 versus 2) but eliminates the derivative entirely, using only values of $f$ you already computed. When $f'$ is expensive, unavailable, or noisy — think of an $f$ that is itself the output of a long simulation — the secant is often the better choice, because the right measure of cost is work per correct digit, not order alone. Counting function evaluations, Newton needs $f$ and $f'$ each step while the secant needs only one new $f$ , so their efficiency indices ( $p^{1/(\text{evals})}$ ) are $\sqrt2\approx1.41$ for Newton versus $1.618$ for the secant. By that honest measure the secant can actually win. Its cousin, the derivative-free Brent method, is what most libraries ship.

How do I find all the roots of a function, not just one?

No single starting guess can guarantee you every root — open methods converge to whichever root their basin of attraction leads to, and those basins can be fractally tangled. The reliable strategy is to bracket first: scan the domain on a grid fine enough to catch every sign change of $f$ , and each sign change brackets at least one root that bisection or Brent will then isolate. For the transcendental equations of physics — say $\tan(ka)=\kappa/k$ for a finite square well, with infinitely many roots marching off to large $k$ — this is exactly right, because the roots interleave the poles of $\tan$ predictably, giving a natural bracket around each bound state. Newton alone would jump unpredictably between branches; grid-and-bracket extracts the energies one at a time, in order.

What is Brent's method, and why do libraries prefer it?

Brent's method (1973) is the practical resolution of the speed-versus-safety dilemma. It always maintains a bracket $[a,b]$ that straddles the root, so like bisection it can never diverge and always converges. But instead of blindly halving, it fits an inverse quadratic through the three most recent points and jumps to that curve's root — a superlinear step — accepting it only if the result stays inside the bracket and is making adequate progress; otherwise it falls back to a guaranteed bisection step. The result inherits bisection's ironclad convergence and, on smooth functions, achieves nearly the secant's speed. This is precisely the hybrid philosophy the whole subject points toward, and it is why `brentq` is the default one-dimensional root-finder in SciPy, GSL, and Numerical Recipes.

How does root-finding connect to WKB and quantization?

Directly — quantization conditions are root-finding problems. In the WKB approximation the Bohr–Sommerfeld condition $\oint p\,\mathrm dq=(n+\tfrac12)h$ becomes, once the action integral is evaluated, a transcendental equation $S(E)=(n+\tfrac12)\pi\hbar$ for the allowed energy $E$ . Because the classical turning points depend on $E$ , the action $S(E)$ is a nonlinear function with no closed-form inverse, so each energy level is literally a root, extracted by Newton on $S(E)-(n+\tfrac12)\pi\hbar=0$ or, more safely, by bisection between successive levels. The same is true of exact bound-state conditions like $\tan(ka)=\kappa/k$ . This is why root-finding sits directly upstream of the WKB chapter: every energy eigenvalue you will compute there is the output of an iteration of exactly the kind derived on this page.