Errors and root finding

Computers approximate real numbers, so every algorithm must control error.

1. Floating-point numbers and rounding (Def. 1.1); relative error and cancellation (Example 1.2).
2. Bisection (Prop. 1.3).
3. Fixed-point iteration and its convergence (Thm 1.4); order of convergence (Def. 1.5).
4. Newton's method converges quadratically (Thm 1.6); the secant method (Prop. 1.7).

A binary floating-point number has the form $\pm(1.d_1d_2\dots d_t)_2\times2^e$ with $t$ mantissa bits and exponent $e$ in a range. Rounding to nearest gives $\operatorname{fl}(x) = x(1 + \delta)$ with $|\delta|\le u = 2^{-t-1}$, the unit roundoff ($\varepsilon_{\text{mach}} = 2^{-t}$; IEEE double precision has $t = 52$, $u\approx1.1\times10^{-16}$). Each arithmetic operation satisfies $\operatorname{fl}(a\circ b) = (a\circ b)(1 + \delta)$, $|\delta|\le u$.

Subtracting nearly equal numbers loses relative accuracy: if $a\approx b$ carry relative errors $u$, then $a - b$ has relative error up to $\frac{u(|a| + |b|)}{|a - b|}$. Examples and remedies:

• $\sqrt{x + 1} - \sqrt x$ for large $x$: use $\frac{1}{\sqrt{x + 1} + \sqrt x}$.
• Roots of $x^2 - bx + c$ with $b^2\gg|c|$: compute $x_1 = \frac{b + \operatorname{sign}(b)\sqrt{b^2 - 4c}}{2}$, then $x_2 = \frac c{x_1}$.
• $1 - \cos x$ for small $x$: use $2\sin^2\frac x2$.

If $f$ is continuous on $[a, b]$ with $f(a)f(b) < 0$, bisection (halve the interval, keep the half with a sign change) produces intervals $[a_n, b_n]$ of length $\frac{b - a}{2^n}$ containing a root. To guarantee error below $\varepsilon$ requires $n\ge\log_2\frac{b - a}\varepsilon$ steps.

Let $g$ map $[a, b]$ into itself with $|g'(x)|\le k < 1$ there. Then $g$ has a unique fixed point $x^*\in[a, b]$, and $x_{n+1} = g(x_n)$ converges to it from any $x_0\in[a, b]$, with $|x_n - x^*|\le k^n|x_0 - x^*|$. If $g'(x^*)\ne0$, then $\frac{e_{n+1}}{e_n}\to g'(x^*)$.

A sequence $x_n\to x^*$ converges with order $p\ge1$ if $\lim\frac{|e_{n+1}|}{|e_n|^p} = C$ with $0 < C < \infty$ (and $C < 1$ if $p = 1$). $p = 1$ is linear, $p = 2$ quadratic. Fixed-point iteration with $g'(x^*) = \dots = g^{(p-1)}(x^*) = 0\ne g^{(p)}(x^*)$ has order $p$.

Newton's method $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$ is fixed-point iteration with $g(x) = x - \frac{f}{f'}$. If $f$ is $C^2$ near a simple root $x^*$ ($f'(x^*)\ne0$), then for $x_0$ close enough to $x^*$ it converges quadratically:

At a root of multiplicity $m > 1$ it converges only linearly, with ratio $1 - \frac1m$.

The secant method replaces $f'(x_n)$ by the slope through the last two iterates: $x_{n+1} = x_n - f(x_n)\frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}$. Near a simple root, $e_{n+1}\approx\frac{f''(x^*)}{2f'(x^*)}e_ne_{n-1}$, giving order $\phi = \frac{1 + \sqrt5}2\approx1.618$. It needs one new function value per step and no derivative.

• Bracket first (sign change) for safety; then switch to Newton or secant for speed.
• Estimate the order from successive errors: $p\approx\frac{\log|e_{n+1}/e_n|}{\log|e_n/e_{n-1}|}$.
• Stopping: $|x_{n+1} - x_n|$ small and $|f(x_n)|$ small; beware stalls at roundoff level.
• Multiple roots slow Newton down to linear.