如何使用numpy(和scipy)查找函数的全零？

Question

假设我f(x)在a和之间定义了一个函数b.此函数可以有很多零,但也有很多渐近线.我需要检索此函数的所有零.最好的方法是什么？

实际上,我的策略如下:

1. 我在给定数量的点上评估我的函数
2. 我检测到是否有符号变化
3. 我发现正在改变符号的点之间的零点

4. 我验证找到的零是否真的为零,或者这是否为渐近线

   U = numpy.linspace(a, b, 100) # evaluate function at 100 different points
   c = f(U)
   s = numpy.sign(c)
   for i in range(100-1):
       if s[i] + s[i+1] == 0: # oposite signs
           u = scipy.optimize.brentq(f, U[i], U[i+1])
           z = f(u)
           if numpy.isnan(z) or abs(z) > 1e-3:
               continue
           print('found zero at {}'.format(u))
   

这个算法似乎有效,除了我看到两个潜在的问题:

1. 它不会检测到没有穿过x轴的零(例如,在函数中f(x) = x**2)但是,我认为它不会出现在我正在评估的函数中.
2. 如果离散点太远,它们之间可能会有一个零,并且算法可能无法找到它们.

您是否有更好的策略(仍然有效)来查找函数的所有零？

---

我不认为这个问题很重要,但对于那些好奇的人,我正在处理光纤中波传播的特征方程.该函数看起来像(之前V和ell之前定义的,ell是一个正整数):

def f(u):
    w = numpy.sqrt(V**2 - u**2)

    jl = scipy.special.jn(ell, u)
    jl1 = scipy.special.jnjn(ell-1, u)
    kl = scipy.special.jnkn(ell, w)
    kl1 = scipy.special.jnkn(ell-1, w)

    return jl / (u*jl1) + kl / (w*kl1)

Answer 1

我看到的主要问题是你是否真的可以找到所有的根——正如评论中已经提到的，这并不总是可能的。如果您确定您的功能不是完全病态的（sin(1/x)已经提到过），那么下一个是您对丢失一个或多个根的容忍度。换句话说，这取决于你准备走多远以确保你没有错过任何一个——据我所知，没有通用的方法可以为你隔离所有的根源，所以你必须自己做。您所展示的已经是合理的第一步。几点评论：

• 布伦特的方法在这里确实是一个不错的选择。
• 首先，处理分歧。由于在您的函数中，分母中有贝塞尔，因此您可以首先求解它们的根- 更好地在例如阿布拉莫维奇和 Stegun（ Mathworld 链接）中查找它们。这比使用您正在使用的临时网格更好。
• 您可以做什么，一旦找到两个根或分歧，x_1并且x_2，在区间 中再次运行搜索[x_1+epsilon, x_2-epsilon]。继续下去，直到找不到更多的根（布伦特方法保证收敛到一个根，前提是有一个根）。
• 如果您无法枚举所有分歧，您可能需要更加小心地验证候选者确实是分歧：给定的不要x只是检查它f(x)是否很大，还要检查它，例如|f(x-epsilon/2)| > |f(x-epsilon)|对于epsilon (1e-8, 1e -9、1e-10，类似的东西）。
• 如果您想确保没有直接接触零的根，请查找函数的极值，并针对每个极值 ，x_e检查 的值f(x_e)。