optimization - root finding in python -

optimization - root finding in python -

edit: big issue here scipy.optimize.brentq requires limits of search interval have opposite sign. if piece search interval arbitrary sections , run brentq on each section below , dan in comments, wind throwing lot of useless valueerrors. there slick way deal in python?

original post: i'm repeatedly searching functions largest 0 in python. right i'm using scipy.optimize.brentq root , using brutish search method if initial bounds don't work:

#function find largest root of f def bigroot(func, pars): try: root = brentq(func,0.001,4,pars) except valueerror: s = 0.1 while true: try: root = brentq(func,4-s,4,pars) break except valueerror: s += 0.1 go on homecoming root

there 2 big problems this.

first assume if there multiple roots in interval, brentq homecoming largest. did simple testing , i've never seen homecoming except largest root don't know if that's true in cases.

the sec problem in script i'm using function homecoming 0 in cases though function pass bigroot diverges @ 0. if alter step size of search 0.1 0.01 homecoming constant nonzero value in cases. realize details of depend on function i'm passing bigroot think problem might way i'm doing search.

the question is, what's smarter way largest root of function in python?

thanks dan; little more info below requested.

the functions i'm searching behaved in regions i'm interested in. illustration plotted below (code @ end of post).

the singular point @ 0 (the peak goes off top of plot finite) , there either 2 or 3 roots. largest root isn't greater 1 never run away infinity. intervals between roots smaller @ low end of domain they'll never become extremely little (i'd they'll larger 10^-3).

from numpy import exp e #this isn't function plotted def v(r): homecoming 27.2*( 23.2*e(-43.8*r) + 8.74e-9*e(-32.9*r)/r**6 - 5.98e-6*e(-0.116*r)/r**4 + 0.0529*( 23*e(-62.5*r) - 6.44*e(-32*r) )/r - 29.3*e(-59.5*r) ) #this definition of function in plot def f(r,b,e): homecoming 1 - b**2/r**2 - v(r)/e #the plot of f(r,0.1,0.06)

good question, it's math problem rather python problem.

in absence of analytic formula roots of function, there's no way guarantee you've found largest root of function, on given finite interval. example, can build function oscillates between ±1 faster , faster approaches 1.

f(x) = sin(1/(1-x))

this bork numerical method tries find largest root on interval [0,1), since root there larger ones in interval.

so you'll have give background characteristics of functions in question in order more insight general problem.

update: looks functions well-behaved. brentq docs suggest there no guarantee of finding largest/smallest root in interval. seek partitioning intervals , recursively searching smaller , larger other roots.

from scipy.optimize import brentq # function should recursively find roots in interval # , homecoming them ordered smallest largest. scipy.optimize import brentq def find_all_roots(f, a, b, pars=(), min_window=0.01): try: one_root = brentq(f, a, b, pars) print "root @ %g in [%g,%g] interval" % (one_root, a, b) except valueerror: print "no root in [%g,%g] interval" % (a, b) homecoming [] # no root in interval if one_root-min_window>a: lesser_roots = find_all_roots(f, a, one_root-min_window, pars) else: lesser_roots = [] if one_root+min_window<b: greater_roots = find_all_roots(f, one_root+min_window, b, pars) else: greater_roots = [] homecoming lesser_roots + [one_root] + greater_roots

i tried on function , finds largest root, @ ~0.14.

there's tricky brentq, though:

print find_all_roots(sin, 0, 10, ()) root @ 0 in [0,10] interval root @ 3.14159 in [0.01,10] interval no root in [0.01,3.13159] interval no root in [3.15159,10] interval [0.0, 3.141592653589793]

the sin function should have roots @ 0, π, 2π, 3π. approach finding first two. realized problem right there in docs: f(a) , f(b) must have opposite signs. appears all of scipy.optimize root-finding functions have same requirement, partitioning intervals arbitrarily won't work.

python optimization numpy