Skip to content
Numerical Calculus

Root finding: bisection, Newton, secant and fixed point

Bisection, Newton's method with proof of quadratic convergence, the secant method, fixed-point iteration, Newton for systems and SPICE.

4 min read 0 comments

Finding where a function is zero is one of the most frequent problems: the operating point of a diode in a circuit (where the diode’s exponential equation meets the load line), the frequency at which a filter attenuates by 3 dB, the launch angle to hit a target. For almost every nonlinear equation there is no formula, so one proceeds by iteration: start from an estimate and improve it until it is accurate enough.

The bisection method

If f is continuous on [a, b] and f(a) and f(b) have opposite signs, by the intermediate value theorem of Analysis 1 there is at least one root. Bisection halves the interval at each step:

  1. compute the midpoint c = (a + b)/2;
  2. if f(a) and f(c) have opposite signs the root is in [a, c], otherwise in [c, b];
  3. repeat until the interval is smaller than the tolerance.

After n steps the interval is (b - a)/2^n wide. Going from an interval of 1 to a tolerance of 10^{-6} takes 20 steps (2^{20} \approx 10^6). The method is slow but it always converges, and the error bound is guaranteed. It gains one decimal digit every 3.3 steps or so.

Newton’s method

The idea is to replace the function with its tangent line and take as the new estimate the point where the tangent crosses the axis:

    \[ x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)} \]

The first iterations of Newton's method with tangents approaching the root
Starting from x_0, the tangent to the graph at x_0 crosses the axis at x_1; the tangent at x_1 gives x_2, and so on. The estimates approach the root α very quickly.

Example: the square root. With f(x) = x^2 - a the formula becomes x_{k+1} = \frac12\left(x_k + \frac{a}{x_k}\right), the Babylonian method. For \sqrt2 starting from x_0 = 3 (as in the figure): 1.8333; 1.4621; 1.41500; 1.4142138; 1.41421356237311. The correct decimal digits go 0, 1, 3, 6, 13: they roughly double at every step.

Proof of quadratic convergence

Let \alpha be the root and e_k = x_k - \alpha the error. Expanding f with Taylor around x_k, evaluated at \alpha:

    \[ 0 = f(\alpha) = f(x_k) - f'(x_k)\,e_k + \tfrac12 f''(\xi)\,e_k^2 \]

Dividing by f'(x_k) and using the definition of the iteration, x_k - \frac{f(x_k)}{f'(x_k)} - \alpha = e_{k+1}:

    \[ e_{k+1} = \frac{f''(\xi)}{2f'(x_k)}\,e_k^2 \]

The new error is proportional to the square of the old one: if e_k = 10^{-4}, e_{k+1} is of the order of 10^{-8}. This is quadratic convergence, far faster than bisection, which is linear.

Newton’s limits

  • It needs the derivative, which is not always available.
  • Convergence is guaranteed only when starting close enough to the root. Far away, the tangents can send the estimate anywhere or make it oscillate: with f(x) = \arctan x and x_0 = 2 the estimates diverge.
  • If f'(\alpha) = 0 (a double root) convergence becomes only linear.

In practice the two methods are combined: bisection until the interval is small, then Newton to finish in a few steps. This is what functions such as MATLAB’s fzero or SciPy’s brentq do.

The secant method

When the derivative is not available it is approximated by the difference quotient of the last two estimates:

    \[ x_{k+1} = x_k - f(x_k)\,\frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})} \]

Convergence is superlinear, of order \frac{1 + \sqrt5}{2} \approx 1.618 (the golden ratio), and each step needs only one evaluation of f: it is often more efficient than Newton.

Fixed-point iteration

Every equation f(x) = 0 can be rewritten as x = g(x), and a fixed point sought with x_{k+1} = g(x_k). It converges if |g'(x)| \lt 1 near the solution: the error is multiplied by about |g'(\alpha)| at each step.

Example. x = \cos x: starting from any value and repeatedly pressing the cosine key of a calculator (in radians) you reach 0.739085…, because |g'| = |\sin x| \lt 1 near the solution. Newton itself is a fixed-point iteration with g(x) = x - f/f', which has g'(\alpha) = 0: that is why it converges faster.

Nonlinear systems

For a system \mathbf F(\mathbf x) = \mathbf 0 of n equations in n unknowns, Newton generalises with the Jacobian matrix of Analysis 2:

    \[ J(\mathbf x_k)\,\boldsymbol\Delta_k = -\mathbf F(\mathbf x_k),\qquad \mathbf x_{k+1} = \mathbf x_k + \boldsymbol\Delta_k \]

at each step a linear system is solved. It is the heart of SPICE: for a circuit with diodes and transistors, the simulator linearises the components around the current estimate and solves the linear system, repeating until the voltages stop changing. When a simulation “fails to converge”, it is Newton failing to find the solution from the initial guess.

Stopping criteria

Stop when the difference between two successive estimates is below a tolerance, |x_{k+1} - x_k| \lt \text{tol}, or when the residual |f(x_k)| is small enough, and in any case after a maximum number of iterations, so as not to get stuck if the method does not converge. A small residual alone is not enough: if the function is very flat near the root, |f| can be tiny even far from the solution.

Related articles

Subscribe
Notify of
guest
0 Comments
Oldest
Newest Most Voted
0
Would love your thoughts, please comment.x
()
x