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

@strimp099: f(x) is totally different from func(x). So which one are you trying to solve? — Roland Smith, Oct 24 '12 at 18:15
This series of message in the Scipy mailing list might be relevant: http://mail.scipy.org/pipermail/scipy-user/2007-September/013870.html — voithos, Oct 24 '12 at 18:43

halex · Answer 1 · 2015-02-04T08:22:40.693

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].

thanks for the script. in python version 2.x math.ceil returns a float, so you need to convert n to int(n) — Nicor Lengert, Feb 04 '15 at 08:10

score 13 · Accepted Answer · answered Oct 24 '12 at 18:44

13

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])

answered Oct 24 '12 at 18:44

John Vinyard

12,997
3
30
43

How do you pass in the jacobian in such a case if we wanted to?. Just passing in `jac=jac` with an analytical jacobian gives a `TypeError: Shape mismatch`. – Kid Charlamagne Aug 02 '17 at 21:22

score 5 · Answer 3 · answered Oct 24 '12 at 18:45

5

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)).

answered Oct 24 '12 at 18:45

Bitwise

7,577
6
33
50

And if your functions belong to a specific class, there's this: http://openopt.org/interalg – pv. Oct 24 '12 at 21:13

Python: Finding multiple roots of nonlinear equation

3 Answers3

Linked