Python中的根查找

Question

编辑：这里有个大问题，就是 scipy.optimize.brentq 要求搜索区间的边界必须是相反的符号。如果你把搜索区间随便切成几个部分，然后像我下面做的那样，或者像Dan在评论中提到的那样在每个部分上运行 brentq，你会遇到很多没用的ValueErrors。有没有什么聪明的方法可以在Python中处理这个问题呢？

原始帖子：我在Python中反复搜索函数的最大零点。目前我使用 scipy.optimize.brentq 来寻找根，如果我的初始边界不行，我就用一种粗暴的搜索方法：

#function to find the 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
                continue
    return root

这有两个大问题。

首先，我假设如果一个区间内有多个根，brentq会返回最大的那个。我做了一些简单的测试，发现它从来只返回最大的根，但我不知道这在所有情况下是否都成立。

第二个问题是，在我使用的脚本中，这个函数在某些情况下总是返回零，尽管我传给 bigRoot 的函数在0处是发散的。如果我把搜索的步长从0.1改成0.01，那么在这些情况下它会返回一个恒定的非零值。我意识到这取决于我传给 bigRoot 的函数，但我觉得问题可能出在我搜索的方式上。

我的问题是，在Python中寻找函数的最大根，有什么更聪明的方法吗？

---

谢谢Dan；根据要求，下面是更多信息。

我搜索的函数在我感兴趣的区域表现良好。下面是一个示例图（代码在帖子末尾）。

唯一的奇点在0处（图的顶部的峰值是有限的），通常有两个或三个根。最大的根通常不大于1，但从来不会像逃跑到无穷大那样。根之间的间隔在域的低端变小，但它们永远不会变得极小（我认为它们总是会大于10^-3）。

from numpy import exp as e
#this isn't the function I plotted
def V(r):
    return  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 is the definition of the function in the plot
def f(r,b,E):
    return 1 - b**2/r**2 - V(r)/E

#the plot is of f(r,0.1,0.06)

Answer 1

这是个好问题，但其实是个数学问题，而不是Python问题。

如果没有一个公式可以直接算出函数的根（也就是函数等于零的点），那么就没办法保证你找到的是那个函数的最大根，即使是在一个特定的有限区间内。举个例子，我可以构造一个函数，它在接近1的时候，来回摆动在±1之间，而且摆动得越来越快。

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

这样的话，任何试图在区间[0,1)内找到最大根的数值方法都会出问题，因为对于任何一个根，区间内总是存在更大的根。

所以你需要提供一些关于这些函数特性的背景信息，才能对这个普遍的问题有更多的了解。

更新：看起来这些函数的表现还不错。brentq的文档提到，在这个区间内找不到最大或最小根是没有保证的。你可以尝试把区间分开，然后递归地寻找更小或更大的其他根。

from scipy.optimize import brentq

# This function should recursively find ALL the roots in the interval
# and return them ordered from smallest to largest.

from 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 at %g in [%g,%g] interval" % (one_root, a, b)
    except ValueError:
        print "No root in [%g,%g] interval" % (a, b)
        return [] # No root in the 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 = []

    return lesser_roots + [one_root] + greater_roots

我在你的函数上试过这个方法，找到了最大根，大约在0.14的位置。

不过，brentq有点棘手：

print find_all_roots(sin, 0, 10, ())

Root at 0 in [0,10] interval
Root at 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]

比如sin函数应该在0、π、2π、3π这些地方有根。但这个方法只找到了前两个。我意识到问题就在于文档中提到的：f(a)和f(b)必须是相反的符号。看起来所有的scipy.optimize根查找函数都有这个要求，所以随便分割区间是行不通的。