SciPy - 10 (Root finding)

A number of root finding tools are available in scipy's optimize module. If you try to find a x such that f(x) = 0; then this is called root finding. Note that problems like g(x) = h(x) fall in this category as you can rewrite them as f(x) = g(x) - h(x) = 0.

Root finding using the bisection method

Let's compute the roots of f(x) = x³-2x² . This function has a (double) root at x = 0 and another root which is located between x = 1:5 (where f(1:5) = -1:125) and x = 3 (where f(3) = 9). It is pretty straightforward to see that this other root is located at x = 2. The following program determines this root numerically:

from scipy . optimize import bisect
def f(x):
      """ returns f(x)=x^3 -2x^2. Has roots at
      x=0 ( double root ) and x=2 """
      return x ** 3 - 2 * x ** 2
# main program starts here
x = bisect (f, 1.5 , 3, xtol =1e -6)
print "The root x is approximately x =%14.12g ,\n"\"the error is less than 1e -6." % (x)
print "The exact error is %g." % (2 - x)

The output of the program is as follows:

The root x is approximately x= 2.00000023842 ,
the error is less than 1e -6.
The exact error is -2.38419e -07.

------------------

(program exited with code: 0)

Press any key to continue . . .

The bisect() method takes three compulsory arguments: (i) the function f(x), (ii) a lower limit a (for which we have chosen 1.5 in our example) and (ii) an upper limit b (for which we have chosen 3). The optional parameter xtol determines the maximum error of the method.

One of the requirements of the bisection method is that the interval [a; b] has to be chosen such that the function is either positive at a and negative at b, or that the function is negative at a and postive at b. In other words: a and b have to enclose a root.

Root finding using the fsolve funcion

This algorithm is implemented in fsolve() function for root finding of (multidimensional) functions. It needs only one starting point close to the suspected location of the root. See the following example:

from scipy . optimize import fsolve
def f(x):
    return x ** 3 - 2 * x ** 2
x = fsolve (f, 3) # one root is at x =2.0
print "The root x is approximately x =%21.19 g" % x
print "The exact error is %g." % (2 - x)

The output of the program is as follows:

The root x is approximately x= 2.000000000000006661
The exact error is -6.66134e -15.

------------------

(program exited with code: 0)

Press any key to continue . . .

The return value of fsolve() is a numpy array of length n for a root nding problem with n variables. In the example above, we have n = 1.

This is the end of today's post. Will catch up again with a new topic! Till then keep practicing and learning Python as Python is easy to learn!