Solving equations

Many mathematical problems eventually reduce to solving an equation of the formf(x) = 0, where f is a function of a single variable. Here, we try to find a value of x for which the equation holds. The values of x for which the equation holds are sometimes called roots of the equation. There are numerous algorithms for finding solutions to equations of this form. In this recipe, we will use the Newton-Raphson and secant methods to solve an equation of the form f(x) = 0.

The Newton-Raphson method (Newton's method) and the secant method are good, standard root finding algorithms that can be applied in almost any situation. These areiterative methods that start with an approximation of the root and iteratively improve this approximation until it lies within a given tolerance.

To demonstrate these techniques, we will use the function from the Symbolic calculus using SymPy recipe defined by

which is defined for all real values of x and has exactly two roots, one atx = 0 and one atx= 2.

Getting ready

The SciPy package contains routines for solving equations (among many other things). The root finding routines can be found in the optimizemodule from the scipy package.

If your equation is not in the form f(x) = 0, then you will need to rearrange it so that this is the case. This is usually not too difficult, and simply requires moving any terms on the right-hand side over to the left-hand side. For example, if you wish to find the fixed points of a function, that is, wheng(x)= x, then we would apply the method to the related function given by f(x) =g(x)- x.

How to do it...

The optimize package provides routines for numerical root finding. The following instructions describe how to use the newton routine from this module:

1. The optimizemodule is not listed in the scipynamespace, so you must import it separately:

from scipy import optimize

1. Then we must define this function and its derivative in Python:

from math import exp

def f(x):
    return x*(x - 2)*exp(3 - x)

1. The derivative of this function was computed in the previous recipe:

def fp(x):
    return -(x**2 - 4*x + 2)*exp(3 - x)

1. For both the Newton-Raphson and secant methods, we use the newton routine from optimize. Both the secant method and the Newton-Raphson method require the function and the first argument and the first approximation,x0, as the second argument. To use the Newton-Raphson method, we must provide thederivative off, using the fprimekeyword argument:

optimize.newton(f, 1, fprime=fp) # Using the Newton-Raphson method
# 2.0

1. To use the secant method, only the function is needed, but we must provide the first two approximations for the root; the second is provided as the x1keyword argument:

optimize.newton(f, 1., x1=1.5) # Using x1 = 1.5 and the secant method
# 1.9999999999999862

How it works...

The Newton-Raphson method for a functionf(x) with derivativef'(x) and initial approximation x₀is defined iteratively using the formula

for each integeri ≥0. Geometrically, this formula arises by considering the direction in which the gradient is negative (so the function is decreasing) iff(x_i)>0 or positive (so the function is increasing) iff(x_i) <o.

The secant method is based on the Newton-Raphson method, but replaces the first derivative by the approximation

whenx_i-x_i-1 is sufficiently small, which occurs if the method is converging, then this is a good approximation. The price paid for not requiring the derivative of the functionf is that we require an additional initial guess to start the method. The formula for the method is given by

Generally speaking, if either method is given an initial guess (guesses for the secant method) that is sufficiently close to a root, then the method will converge to that root. The Newton-Raphson method can also fail if the derivative is zero at one of the iterations, in which case the formula is not well defined.

There's more...

The methods mentioned in this recipe are general purpose methods, but there are others that may be faster or more accurate in some circumstances. Broadly speaking, root finding algorithms fall into two categories: algorithms that use information about the function's gradient at each iterate (Newton-Raphson, secant, Halley) and algorithms that require bounds on the location of a root (bisection method, regula-falsi, Brent). The algorithms discussed so far are of the first kind, and while generally quite fast, they may fail to converge.

The second kind of algorithms are those for which a root is known to exist within a specified intervala ≤x ≤b. We can check whether a root lies within such an interval by checking thatf(a) andf(b) have different signs, that is, one of f(a) <0<f(b) orf(b) <0<f(a) is true. (Provided, of course, that the function iscontinuous, which tends to be the case in practice.) The most basic algorithm of this kind is the bisection algorithm, which repeatedly bisects the interval until a sufficiently good approximation to the root is found. The basic premise is to split the interval betweenaand bat the mid-point and select the interval in which the function changes sign. The algorithm repeats until the interval is very small. The following is a rudimentary implementation of this algorithm in Python:

from math import copysign

def bisect(f, a, b, tol=1e-5):
    """Bisection method for root finding"""
    fa, fb = f(a), f(b)
    assert not copysign(fa, fb) == fa, "Function must change signs"
    while (b - a) > tol:
        m = (b - a)/2 # mid point of the interval
        fm = f(m)
        if fm == 0:
            return m
        if copysign(fm, fa) == fm: # fa and fm have the same sign
            a = m
            fa = fm
        else: # fb and fm have the same sign
            b = m
        return a

This method is guaranteed to converge, since at each step the distanceb-a is halved. However, it is possible that the method will require more iterations than Newton-Raphson or the secant method. A version of the bisection method can also be found inoptimize. This version is implemented in C and is considerably more efficient that the version presented here, but the bisection method is not the fastest method in most cases.

Brent's method is an improvement on the bisection method, and is available in theoptimize module asbrentq. It uses a combination of bisection and interpolation to quickly find the root of an equation:

optimize.brentq(f, 1.0, 3.0)  # 1.9999999999998792

It is important to note that the techniques that involve bracketing (bisection, regula-falsi, Brent) cannot be used to find the root functions of a complex variable, whereas those techniques that do not use bracketing (Newton, secant, Halley) can.