Root-Finding Algorithms¶

Here we discuss some of the numerical methods for finding the roots of the equations of the form f(x) = 0, where 'f' is a function. We can find the roots in closed form for only some limited class of equations. Most of the time it is difficult to find roots of transedental equations. Sometimes, we can use graphical techniques to have rough idea of roots for such equations.

There are a variety of techniques to estimate the roots of the equation depending on the nature of the function involved. Broadly we classify the different methods into two categories:

Root Bracketing Methods
Iterative Methods

Bracketing method is used when the root 'c' of f(x) = 0 is known to lie in the interval (a,b). A possible guess of interval (a,b) for a given equation can be made by checking the signs of f(a) and f(b). If signs are opposite then there is a possibility of root of the equation.

Bisection Method is one of the simplest numerical techniques based on root bracketing methods. This method works if function 'f' is continuous and f(a)*f(b) < 0. We find the sign of f(x) at the mid of the interval (a,b) i.e., c = (a+b)/2 and check the sign of f(c). If it matches with the sign of f(a) the sign of then the interval (a,b) is replaced with (c,b) else with (a,b). This process is repeated and it can be stopped by one of the following mechanisms:

the absolute value of the difference between two consecutive estimates of mid-points of the interval satisfies some tolerance limit.
the absolute value of the function 'f' at the mid-point of the interval satisfies some tolerance limit.

Here is a sample python code to estimate the one of the roots of $$ 2x^4-7x^3-15x^2+34x+40 = 0 $$ using Bisection method. Clearly, $ f(x) = 2x^4-7x^3-15x^2+34x+40 $ and (a,b) can be taken as (0,3).

from __future__ import  division
import numpy as np

def f(x):
	return 2*x**4-7*x**3-15*x**2+34*x+40

a = input ("Enter the lower bound of the interval in which root may lie:")
b = input ("Enter the upper bound of the interval in which root may lie:")

c = (a+b)/2

if(f(a)*f(b) > 0) :
	print a, "  and ",b," are not appropriate for bisection method."
else: 
	while ((b-a) > 0.001):
		if(f(a)*f(c) < 0):
			b=c
			c=(a+b)/2
		else :
			a=c
			c=(a+b)/2
	print c," is the approximate root."

Enter the lower bound of the interval in which root may lie:0
Enter the upper bound of the interval in which root may lie:3
2.50012207031  is the approximate root.

Let us plot the above function $ f(x) = 2x^4-7x^3-15x^2+34x+40 $ to see the location of roots.

%matplotlib inline
import matplotlib.pyplot as plt

x=np.linspace(-3,5,100)
plt.xlabel('x')
plt.ylabel(r'$f(x)$')
plt.plot(x,2*x**4-7*x**3-15*x**2+34*x+40)
plt.plot(x,x*0)
plt.axis([-3,5,-30,60])
plt.show()

There is an another technique known as Regula-Falsi Method under Root Bracketing Method which has better convergence than Bisection Method. Here to choose the point 'c', we also consider the role of function f(x) along with chosen interval (a,b),

$$c = b-\frac{(b-a)*f(b)}{f(b)-f(a)} = \frac{a*f(b)-b*f(a)}{f(b)-f(a)}$$

Now let us discuss some Iterative Methods for numerical approximation of roots. Here, we start with a guess of the root '$x_0$' and using some recursion relation we get better and better approximations of the root. The method works if the approximate roots generated recursively converge.

Newton's Method or Newton-Raphson Method is one of such methods. Here to find the next approximation of root, stating with $x_0$, we calculate the x-intercept of tangent line to the curve f(x) at the point corresponding to $x_0$. $$x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}$$ In general the recursion relation can be written as, $$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$$ This recursion can be stopped by

fixing apriori the total number of iterations
setting a tolerance limit for $|x_{n+1}-x_n|$

This method works fine if f(x) is differentiable and convergence is better if the initial guess $x_0$ is close enough to actual unknown root.

Here is a sample python code to estimate the root of $x- exp(-x^2) = 0 $, using Newton's method. Here, $ f(x) = x- exp(-x^2) $, $ f'(x) = 1+2x* exp(-x^2) $ and the guess $x_0$ can be taken as 1.0

from __future__ import  division
import numpy as np

def f(x):
	return x-np.exp(-x**2)

def f_(x):
	return 1+2*x*np.exp(-x**2)     # or one can estimate derivative with 
                                   #     (f(x+h)-f(x-h))/(2*h)
 
x1 = input ("Enter the guess value of root:")

x2 = x1-(f(x1)/f_(x1))

while (abs(x2-x1)> 0.001):
	x1=x2
	x2 = x1-(f(x1)/f_(x1))
	
print x2, "is the approximate root."

Enter the guess value of root:1.0
0.652918640437 is the approximate root.

The plot of $f(x)=x-exp(-x^2)$ looks like

%matplotlib nbagg
import matplotlib.pyplot as plt
import numpy as np

x=np.linspace(-2,5,1000)
plt.xlabel('x')
plt.ylabel(r'$f(x)$')
plt.plot(x,x-np.exp(-x**2))
plt.plot(x,x*0)
plt.show()

There are many other iterative methods of root approximation like Secant Method, Muller's Method, which may be more efficient than Newton's Method. In Secant method, we use method of finite difference to evaluate the derivative (which can be helpful if $f(x)$ has complicated form).

Vinay Vaibhav

Root Finding Algorithms

Root-Finding Algorithms¶