Solution by Linear Interpolation

The bi-section method is very simple, but generally quite inefficient, in part because it only makes use of the sign of the function f(x) at each evaluation, while ignoring its magnitude. Thus it ignores significant information which could be used accelerate the finding of the root. A method based on interpolation makes use of this information by approximating the function on the interval [x₁,x₂] by the chord joining the points (x₁,f(x₁)) and (x₂,f(x₂)), that is the straight line:

$$\frac{y - y_1}{x - x_1} = \frac{y_2 - y_1}{x_2 - x_1}.$$ (6.1)

Solving this linear equation for y = 0 yields the new interval endpoint within the interval [x₁,x₂]:

$$x_3 = x_1 - f(x_1) \frac{x_2 - x_1}{f(x_2) - f(x_1)}.$$ (6.2)

The choice between the two intervals [x₁,x₃] and [x₃,x₂] is decided by evaluating f(x₃) and discarding the interval whose endpoints have the same sign, as was done in the bi-section method. This iteration process is repeated, but it will converge more quickly than the bi-section method, since the information about the magnitude of f(x) pushes the x₃ value more quickly towards the actual root. The iteration is continued until the end points of the interval show no change to the required precision in subsequent iterations. One important point to note is that frequently the interval does not become small, as it did in the bisection method, since one end may remain fixed while the other end closes in on the zero from one side only. The example shown in figure 6.1 illustrates this effect.

This method is also known as the regula falsi method. This method is also a bracketing method, since it also simply refines the definition of the interval containing the zero. We shall see that a related method without the bracketing requirement exists, called the secant method, where the linear interpolation becomes linear extrapolation.

Figure 6.1: Successive steps in the regula falsi iteration to find the zero of f(x). Notice that in this example the right hand end of the test interval remains fixed.

#!/usr/bin/env python

def sign(x):                   # determines the sign of its argument
  if x == abs(x) : return 1    # argument was positive or zero
  else: return -1              # argument was negative

# Solve f = 0 on interval [x1,x2] by interpolation, with tolerances
def interpol_solve(f,x1,x2,ftol,xtol):
  f1 = f(x1)
  if abs(f1) <= ftol : return x1
  s1 = sign(f1)
  f2 = f(x2)
  if abs(f2) <= ftol : return x2
  s2 = sign(f2)
  if s1 == s2 :
    sys.stderr.write("Same sign at %g to %g - exit!\n" % (x1,x2))
sys.exit(1)
  while abs(x2 - x1) > xtol :
    x3 = x2 - f2*(x2 - x1)/(f2 - f1)
    f3 = f(x3)
    if abs(f3) <= ftol : break
    s3 = sign(f3)
    if s3 == s1 :
      (x1,f1) = (x3,f3)            # replace pair (x1,f1) by (x3,f3)
    else :
      (x2,f2) = (x3,f3)            # replace pair (x2,f2) by (x3,f3)
  return x3

def quad(x):    # a simple test function with known zeroes
  return (x - 5.0)*(x - 2.0) 

# a simple main to test the regula falsi solver */
root = interpol_solve(quad,1.0,3.0,0.000001,0.000001);
print "ROOT = %g" % root

Figure 6.1 illustrates the regula falsi method and makes clear one fault with the method. If the function f(x) is strongly curved on the starting interval, the approach to the root can be one-sided in the sense that one end of the interval remains fixed rather than the endpoint that is discarded oscillating from one end to the other. This can be partially suppressed by artificially decreasing the value of f(x) at the stagnant end of the interval to f(x)/2 at each iteration. This does not alter the sign of f(x) at that endpoint, but does tend to accelerate the convergence towards the root by inducing a more equal distribution of the interval end which is dropped.