Scipy寻根方法答案

【问题标题】：Scipy root-finding methodScipy寻根方法
【发布时间】：2016-05-06 22:15:53
【问题描述】：

我正在编写一个具有数学函数作为属性的类，例如 f。

f 是：

在实际段上定义 [-w;+w]
正且以实 H 为界
偶数（对于 [-w;+w] 中的所有 x，f(x)=f(-x)）和 f(w)=f(-w)=0
在 [-w;+w] 上可微，其导数在 [-w;0] 上为正且连续

我的班级看起来像：

from scipy.misc import derivative
from scipy.integrate import quad
from math import cosh, sqrt

class Function(object):

    w = 1.
    PRECISION = 1e-6    

    def f(self, x):
        '''This is an example but f could be 
           any math function matching requirments above.
        '''
        return 0.5+1.07432*(1-cosh(x/1.07432)) 

    def deriv_f(self, x):
        return derivative(self.f, x, self.PRECISION)

    def x_to_arc_length(self, x):
        def func(x): 
            return sqrt(1+self.deriv_f(x)**2)
        return quad(func, -self.w, x)[0]

    def arc_length_to_x(self, L):
        bound = [-self.w, self.w]
        while bound[1]-bound[0] > self.PRECISION:
            mid= sum(bound)/2
            bound[(self.x_to_arc_length(mid)-L > 0)] = mid
        return sum(bound)/2

我使用二等分法来反转弧长方法，但我正在考虑将其更改为 scipy.optimize 寻根方法之一以提高速度。我是 scipy 的新手，必须承认我的数学有点生锈...... Scipy 让我在 brentq、brenh、ridder、bisect 和 newton 之间进行选择。

谁能指出最适合这种情况的方法？或者也许有更好的图书馆？

【问题讨论】：

试试看？为准确性设置试验，为速度设置路径，并比较各种方法。小心过早的优化。
@Evert ：澄清一下，f 在我的问题代码中设置为悬链线，但可以设置为符合要求的任何函数。我可以尝试一些函数，但我希望有一个数学方法来解决这个问题。我修改了问题以明确 f 是一个例子。
我知道的很少python，但在数值分析中，通常建议使用布伦特方法来寻找标量函数的根。看起来scipy 教程与this suggestion 一起使用（在链接页面中搜索“根查找”）。牛顿法在特定情况下可能更快，但通常更容易出错。请记住，对于所有这些方法，除了 Newton，您需要在输入中提供 f 更改符号的区间的极值。
感谢 DeltaV，我相信这就是我的问题的答案。我在 Wolfram.mathworld 上找到了这个描述：mathworld.wolfram.com/BrentsMethod.html。这看起来比简单的二分法好多了！ brenth 是二次方brentq 的双曲线版本。看来他们是相当的（docs.scipy.org/doc/scipy-0.14.0/reference/generated/…）。值得对它进行基准测试还是我应该坚持brentq？
@DeltaIV：非常感谢。我会看看我是否有时间研究这个，因为这些方法现在并不重要，还有很多工作要做......你为什么不发布你的 cmets 作为答案，以便我可以接受它？

标签： python scipy

【解决方案1】：

我不是 Python 方面的专家，但我从数值分析中知道，在您列出的方法（布伦特、二分法、里德方法和牛顿-拉夫森）中，布伦特方法通常是通用实标量函数的首选 f 的单个实变量 x。如您所见here，如果 f 是连续的并且该方法应用于区间 [a,b] 且 f(a)f(b)N^2 次迭代，其中 N 是平分以达到给定的公差。

另一方面，Newton's method 在收敛时通常比 Brent 更快地收敛，但在某些情况下它根本不收敛。对于同一个函数，牛顿法可能收敛也可能不收敛，这也取决于起点和根之间的距离。因此，在通用代码中使用风险更大。

关于brentq 和brenth 之间的选择，it looks like 它们应该非常相似，第一个测试更严格。因此，您可以选择brentq，或者，如果您有时间，可以在它们之间做一些基准测试。

【讨论】：

【解决方案2】：

根据 DeltaIV 的回答，我根据上面的示例对不同选项进行了基准测试。

 2000    0.005    0.000   10.395    0.005 diff.py:42(arc_length_to_x_newton)
 2000    0.005    0.000   16.842    0.008 diff.py:36(arc_length_to_x_brenth)
 2000    0.005    0.000   17.141    0.009 diff.py:30(arc_length_to_x_brentq)
 2000    0.005    0.000   26.375    0.013 diff.py:48(arc_length_to_x_ridder)
 2000    0.005    0.000   72.249    0.036 diff.py:54(arc_length_to_x_bisect)

在这种情况下，牛顿法似乎是最快的，可能是因为该函数表现良好（即连续甚至导数以 1 为界）。

不收敛的风险（平稳点、循环或导数问题）不适用于上述函数，所以我最终选择了newton。

【讨论】：