在python中查找根

Question

编辑：这里的一个大问题是scipy.optimize.brentq 要求搜索间隔的限制具有相反的符号。如果您将搜索间隔分成任意部分并在每个部分上运行brentq，就像我在下面所做的那样以及 Dan 在 cmets 中所做的那样，您最终会抛出很多无用的 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 中寻找函数的最大根的更聪明的方法是什么？

---

谢谢丹；根据要求，下面提供了更多信息。

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

唯一的奇异点在 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π。但是这种方法只找到前两个。我意识到问题就在那里in the docs：f(a) 和 f(b) 必须有相反的符号。看来scipy.optimize 根的所有 - 查找函数具有相同的要求，因此任意划分区间是行不通的。