optimization - Python: Finding multiple roots of nonlinear equation

Question

Assume the following function:

f(x) = x * cos(x-4)

With x = [-2.5, 2.5] this function crosses 0 at f(0) = 0 and f(-0.71238898) = 0.

This was determined with the following code:

import math
from scipy.optimize import fsolve
def func(x):
    return x*math.cos(x-4)
x0 = fsolve(func, 0.0)
# returns [0.]
x0 = fsolve(func, -0.75)
# returns [-0.71238898]

What is the proper way to use fzero (or any other Python root finder) to find both roots in one call? Is there a different scipy function that does this?

fzero reference

Answer 1

Define your function so that it can take either a scalar or a numpy array as an argument:

>>> import numpy as np
>>> f = lambda x : x * np.cos(x-4)

Then pass a vector of arguments to fsolve.

>>> x = np.array([0.0, -0.75])
>>> fsolve(f,x)
array([ 0.        , -0.71238898])

Answer 2

I once wrote a module for this task. It's based on chapter 4.3 from the book Numerical Methods in Engineering with Python by Jaan Kiusalaas:

import math

def rootsearch(f,a,b,dx):
    x1 = a; f1 = f(a)
    x2 = a + dx; f2 = f(x2)
    while f1*f2 > 0.0:
        if x1 >= b:
            return None,None
        x1 = x2; f1 = f2
        x2 = x1 + dx; f2 = f(x2)
    return x1,x2

def bisect(f,x1,x2,switch=0,epsilon=1.0e-9):
    f1 = f(x1)
    if f1 == 0.0:
        return x1
    f2 = f(x2)
    if f2 == 0.0:
        return x2
    if f1*f2 > 0.0:
        print('Root is not bracketed')
        return None
    n = int(math.ceil(math.log(abs(x2 - x1)/epsilon)/math.log(2.0)))
    for i in range(n):
        x3 = 0.5*(x1 + x2); f3 = f(x3)
        if (switch == 1) and (abs(f3) >abs(f1)) and (abs(f3) > abs(f2)):
            return None
        if f3 == 0.0:
            return x3
        if f2*f3 < 0.0:
            x1 = x3
            f1 = f3
        else:
            x2 =x3
            f2 = f3
    return (x1 + x2)/2.0

def roots(f, a, b, eps=1e-6):
    print ('The roots on the interval [%f, %f] are:' % (a,b))
    while 1:
        x1,x2 = rootsearch(f,a,b,eps)
        if x1 != None:
            a = x2
            root = bisect(f,x1,x2,1)
            if root != None:
                pass
                print (round(root,-int(math.log(eps, 10))))
        else:
            print ('\nDone')
            break

f=lambda x:x*math.cos(x-4)
roots(f, -3, 3)

roots finds all roots of f in the interval [a, b].

Answer 3

In general (i.e. unless your function belongs to some specific class) you can't find all the global solutions - these methods usually do local optimization from given starting points.

However, you can switch math.cos() with numpy.cos() and that will vectorize your function so it can solve for many values at once, e.g. fsolve(func, np.arange(-10,10,0.5)).