Three of the most widely-used numerical algorithms for finding roots of 1-D functions are the bisection method, the modified Regula Falsi method, which I refer to as Regula Falsi, and Newton's method (also known as the Newton-Raphson method).
The bisection method is efficient for functions that are rough (irregular), while Regula Falsi and Newton's method are efficient for functions that are smooth. Many functions that are encountered in engineering and applied physics are rough at large scale but smooth at small scale. This suggests that an adaptive algorithm that combines bisection with Regula Falsi might perform better than either one of these alone.
It is sometimes necessary to repeatedly solve for roots of functions that have similar general behavior, i.e., one may have a family of functions such that the scale at which the transition from rough to smooth occurs is essentially the same for all members of the family. In such cases, a root-finding algorithm that provides tunable parameters can offer significant run-time performance gains.
I compare the bisection and secant algorithms against a novel root-finding algorithm that is both adaptive and tunable. While no root-finding algorithm is universally best, for a number of simple example functions, my algorithm performs better than either bisection or secant. Python implementations of all of these algorithms are provided.
Practical problems in engineering, science, finance, and other domains often involve the finding of roots, i.e., finding the value or values of \(x\)—the input to a function \(f\) of a single variable—such that the output of the function is zero. A problem in which the desired output is a constant value other than zero, or in which the outputs of two single-variable functions are to be made equal to one another, can be trivially reformulated as a problem in which one seeks the root of a single function. (In a more general form of the problem, the function may have multiple inputs and multiple outputs, and one must find a combination of inputs that simultaneously drives all outputs to zero. I don't consider such problems here).
In some applications of root-finding, the function is expressed in a simple form and directly evaluated. Consider the following geometry problem: "A right triangle has sides of lengths \(a\), \(b\), and \(c\), where \(c= 1\) is the length of the hypotenuse and \(b= 3 a\). Find the larger of the two non-right angles, using neither the Pythagorean Theorem nor inverse trigonometric functions". This problem can be solved by finding the root of the following function of \(\theta\), the unknown angle in degrees:
\(sin((\pi/180) \cdot \theta)- 3 \cdot cos((\pi/180) \cdot \theta)= 0\).
(The answer is approximately 71.56505).
In other applications of root-finding, the function has no closed-form expression. Examples include functions that are evaluated by performing numerical integration or via simulation. In such cases, the evaluation of \(f\) may be expensive, i.e., require a large amount of calculation.
Many numerical algorithms have been developed for finding the roots of a function in one variable. Some of these, such as Newton's method, require calculating the derivative (slope) of the function, either from an analytical formula or via numerical approximation. Such algorithms may fail if the derivative is undefined at one or more points. Here, I consider only algorithms that do not depend on derivatives, i.e., \(f\) must be continuous, but need not have a derivative anywhere.
Inputs to the typical root-finding algorithm include:
The algorithm iterates until either the tolerances are satisfied or the iteration limit has been reached.
Desirable properties of a root-finding algorithm include the following:
The complexity of the root-finding algorithm is typically not very important. When \(f\) is a simple function that can be evaluated with minimal computation, almost any root-finding algorithm that works (converges to the desired root) is acceptable. In the more interesting situation where the evaluation of \(f\) is expensive—this is the case that I consider here—speed of convergence may become critical. As we shall see, properties 2. and 3. are not in general simultaneously achievable, i.e., one may be forced to sacrifice one for the sake of the other.
Inputs to the bisection method (and most other root-finding algorithms) include two starting values of \(x\)—\(a\) and \(b\)—which are sometimes called initial guesses. The bisection algorithm requires that \(f(a)\) and \(f(b)\) have opposite signs. (This can be viewed as a drawback of this algorithm). Under the assumption that \(f\) is continuous, \(a\) and \(b\) bracket a root, i.e., the interval from \(min(a, b)\) to \(max(a, b)\) is guaranteed to contain a root. More precisely, this interval is guaranteed to contain an odd number of roots.
In the \(n\)th iteration, the interval between the current \(a\) and \(b\) is bisected, i.e., the midpoint \(x_n= (a+b)/2\) is calculated. \(f_n= f(x_n)\) is evaluated, the signs of \(f_n f(a)\) and \(f_n f(b)\) are compared to determine which subinterval contains the root, and the values of \(a\) and \(b\) are then updated accordingly. (If the root is exactly on the boundary, i.e., at \(f_n\), it does matter which subinterval we choose). The following figure illustrates the operation of the algorithm.
The algorithm iterates until tolerances are satisfied. Most root-finding algorithms count iterations and terminate if the count exceeds a specified limit. The above figure does not include such a test.
My Python implementation of this algorithm (see download link below) is
called find_root_bisection
.
The secant method, like the bisection method, takes initial guesses \(a\) and \(b\) as inputs. In each iteration, \(f(a)\) and \(f(b)\) are evaluated and a straight line is fitted through the points \((a, f(a))\) and \((b, f(b)\):
\(y= f(a) + (x-a) \cdot \frac{f(b)-f(a)}{b-a}\)
We can now set \(y= 0\) and solve for the value \(c\) where this line intersects the x-axis:
\(c= a - f(a) \cdot \frac{b-a}{f(b)-f(a)}\).
\(a\) and \(b\) are then updated as follows:
\(a \leftarrow b\)
\(b \leftarrow c\).
When the secant method is given a pair \((a, b)\) that brackets a root sufficiently closely—it is difficult to specify how close is close enough—the algorithm will not only converge to that root, but converge more rapidly than bisection search. (I will say more about this later). But, while the bisection method guarantees convergence, the secant method does not.
My Python implementation of this algorithm (see download link below) is
called find_root_secant
.
The Regula Falsi method (Ref. 2) is similar to the secant method, but differs as follows.
If the pair \((a, b)\) brackets a root, the updated pair is chosen according to the following rules:
If the pair \((a, b)\) does not bracket a root, the updated pair is \((b, c)\).
The benefit of this modification is that if at any iteration the algorithm bounds a root, it will continue to bound that same root, i.e., convergence is guaranteed. As we shall see, the price of reliable convergence is that the fast convergence of the secant method is sacrificed.
My Python implementation of this algorithm (see download link below) is
called find_root_Regula_Falsi
.
My root-finding algorithm is a hybrid of the Regula Falsi and bisection methods. Each iteration of the algorithm uses either bisection or Regula Falsi. The following parameters control the behavior of the algorithm:
a
and b
: These are the starting values of \(x\).
contraction_factor
: This parameter controls the application of the bisection method. When a pair of consecutive x's brackets a root, and the ratio of the new interval's length to that of the old exceeds this factor, the algorithm applies the bisection method at the next iteration and then reverts to using Regula Falsi. Allowed values for the contraction_factor
are 0.5 to 1. The default value is 0.7071= sqrt(2), chosen so that the algorithm will take no more than twice as many steps as bisection in the worst case. This value appears to typically work well, but may be tuned as needed.
use_bisection
: In the initial use_bisection
steps—this argument equals zero by default— the algorithm uses the bisection method.
My goals in designing this algorithm were the following:
My Python implementation of this algorithm (see download link below) is called find_root
.
Consider the following function, which has a simple root at x= 0.7:
\(f(x)= 2(x-0.7) + 0.03(x-0.7)^3\)
from find_roots import *
def f(x):
return 2.*(x-0.7) + 0.03*(x-0.7)**3
print("\nTesting `find_root_bisection` ...")
x= find_root_bisection(f, 0.6, 6.0, xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
print("\nTesting `find_root_secant` ...")
x= find_root_secant(f, 0.6, 6.0, xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
print("\nTesting `find_root` with default value of `contraction_factor` ...")
x= find_root(f, 0.6, 6.0, xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
print("\nTesting `find_root` with `contraction_factor= 1.0` ...")
x= find_root(f, 0.6, 6.0, xtol=1e-6, ftol=1e-6, contraction_factor=1.0, verbose=True)
print('%6f' % x)
find_root_secant
achieves the most rapid convergence, followed by find_root
, with find_root_bisection
in last place. The performance of find_root
is slightly improved by using a non-default value of the contraction_factor
parameter.
Consider the following function, which has a root of multiplicity 4 at x= 0.7:
\(f(x)= (x-0.7)^4\)
The function touches but does not cross the x-axis.
from find_roots import *
def f(x):
return (x-0.7)**4
print('\nTesting `find_root_bisection` ...')
try:
x= find_root_bisection(f, 0.6, 6.0, xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
except Exception as ex:
print(ex)
print('\nTesting `find_root_secant` ...')
x= find_root_secant(f, 0.6, 6.0, xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
print('\nTesting `find_root` with default value of `contraction_factor` ...')
x= find_root(f, 0.6, 6.0, xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
find_root_secant
and find_root
converge in exactly the same number of steps, while find_root_bisection
fails because \(f(x_1)\) and \(f(x_2)\) are both positive.
find_root
appears to have no advantage over find_root_secant
.
Consider the following function, which is piecewise-linear with three sections:
\(f(x)= \left \{ \begin{array}{cc} -1 & x < -1 \\ 0 & -1 \le x \le 1 \\ 1 & x > 1 \end{array} \right .\)
from find_roots import *
def f(x):
return numpy.clip(x, -1, 1)
print('\nTesting `find_root_bisection` ...')
try:
x= find_root_bisection(f, 0.6, 6.0, xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
except Exception as ex:
print(ex)
print("\nTesting `find_root_secant` ...")
try:
x= find_root_secant(f, 0.6, 6.0, xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
except Exception as ex:
print(ex)
print("\nTesting `find_root` with default value of `contraction_factor` ...")
x= find_root(f, 0.6, 6.0, xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
find_root_bisection
and find_root_secant
both fail.
find_root
achieves rapid convergence.
Consider the following function: \(f(x)= x e^ \left . -|x| \right .\) \(f()\) is continuous and has a first derivative everywhere. The function has peaks on either side of the single root, which occurs at x= 0.
from find_roots import *
def f(x):
return x * exp(-abs(x))
print('\nTesting `find_root_bisection` ...')
x= find_root_bisection(f, -0.5, 10., xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
print("\nTesting `find_root_secant` ...")
try:
x= find_root_secant(f, -0.5, 10., xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
except Exception as ex:
print(ex)
print("\nTesting `find_root` with default value of use_bisection ...")
x= find_root(f, -0.5, 10., xtol=1e-6, ftol=1e-6, verbose=True)
print('%6f' % x)
print("\nTesting `find_root` with use_bisection=2 ...")
x= find_root(f, -0.5, 10., xtol=1e-6, ftol=1e-6, use_bisection=2, verbose=True)
print('%6f' % x)
find_root_secant
fails.
find_root
and find_root_bisection
both converge.
find_root
with default parameters requires half as many function calls as find_root_bisection
.
find_root
, the number of function calls is modestly reduced by performing two bisection search steps at the outset.
find_root
and
find_root_secant
are more flexible
and convenient to use than
find_root_bisection
because they
do not require that the initial values of x bracket a root.
find_root
is more reliable than
find_root_secant
.
find_root
and
find_root_bisection
are both
applicable, find_root
often
converges more rapidly.