python - 如何使用numpy（和scipy）找到函数的所有零？

Question

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

实际上，我的策略如下：

我在给定数量的点上评估我的功能
我检测是否有符号变化
我发现正在改变符号的点之间的零

我验证找到的零是否真的是零，或者这是否是渐近线

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))

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

它不会检测到不穿过 x 轴的零（例如，在类似的函数中f(x) = x**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)

score 3 · Accepted Answer

你为什么受限于numpy？Scipy 有一个包可以完全满足您的需求：

http://docs.scipy.org/doc/scipy/reference/optimize.nonlin.html

我学到的一课：数值编程很难，所以不要这样做:)

无论如何，如果您对自己构建算法一无所知，scipy我链接的文档页面（需要永远加载，顺便说一句）为您提供了一个开始使用的算法列表。我之前使用过的一种方法是将函数离散化到您的问题所必需的程度。（也就是说，调整 \delta x 使其远小于问题中的特征大小。）这使您可以查找函数的特征（例如符号的变化）。并且，您可以很容易地计算线段的导数（可能从幼儿园开始），因此您的离散函数具有明确定义的一阶导数。因为您已将 dx 调整为小于特征大小，所以可以保证不会错过对您的问题很重要的函数的任何特征。

如果您想知道“特征大小”的含义，请查找函数的某些参数，其单位为长度或 1/长度。也就是说，对于某个函数 f(x)，假设 x 有长度单位，而 f 没有单位。然后寻找乘以 x 的东西。例如，如果要离散化 cos(\pi x)，则乘以 x 的参数（如果 x 具有长度单位）必须具有 1/长度的单位。所以 cos(\pi x) 的特征大小是 1/\pi。如果你的离散化比这小得多，你就不会有任何问题。可以肯定的是，这个技巧并不总是有效，所以你可能需要做一些修补。

score 1 · Accepted Answer

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

布伦特的方法在这里确实是一个不错的选择。
首先，处理分歧。由于在您的函数中分母中有贝塞尔，您可以首先求解它们的根——最好在例如 Abramovitch 和 Stegun 中查找它们（Mathworld 链接）。这将比使用您正在使用的临时网格更好。
你可以做什么，一旦你找到了两个根或分歧，x_1然后x_2在区间内再次运行搜索[x_1+epsilon, x_2-epsilon]。继续直到找不到更多的根（如果有一个根，Brent 的方法可以保证收敛到一个根）。
如果你不能列举出所有的分歧，你可能需要更小心地验证一个候选人确实是一个分歧：给定x不要只检查它f(x)是否很大，检查它，例如|f(x-epsilon/2)| > |f(x-epsilon)|对于epsilon (1e-8, 1e)的几个值-9, 1e-10，类似的东西）。
如果要确保没有简单地接触零的根，请查找函数的极值，对于每个极值x_e，请检查的值f(x_e)。

score 1 · Accepted Answer

我发现使用scipy.optimize.fsolve实现自己的根查找器相对容易。

想法：通过不断地调用改变来从区间(start, stop)和步长中找到任何零。使用相对较小的步长来找到所有的根。stepfsolvex0
只能在一个维度上搜索零（其他维度必须是固定的）。如果您有其他需求，我建议使用sympy来计算解析解。
注意：它可能并不总能找到所有的零，但我看到它给出了相对较好的结果。我也把代码放到了gist中，如果需要我会更新。

import numpy as np
import scipy
from scipy.optimize import fsolve
from matplotlib import pyplot as plt

# Defined below
r = RootFinder(1, 20, 0.01)
args = (90, 5)
roots = r.find(f, *args)
print("Roots: ", roots)

# plot results
u = np.linspace(1, 20, num=600)
fig, ax = plt.subplots()
ax.plot(u, f(u, *args))
ax.scatter(roots, f(np.array(roots), *args), color="r", s=10)
ax.grid(color="grey", ls="--", lw=0.5)
plt.show()

示例输出：

Roots:  [ 2.84599497  8.82720551 12.38857782 15.74736542 19.02545276]

放大：

RootFinder 定义

import numpy as np
import scipy
from scipy.optimize import fsolve
from matplotlib import pyplot as plt


class RootFinder:
    def __init__(self, start, stop, step=0.01, root_dtype="float64", xtol=1e-9):

        self.start = start
        self.stop = stop
        self.step = step
        self.xtol = xtol
        self.roots = np.array([], dtype=root_dtype)

    def add_to_roots(self, x):

        if (x < self.start) or (x > self.stop):
            return  # outside range
        if any(abs(self.roots - x) < self.xtol):
            return  # root already found.

        self.roots = np.append(self.roots, x)

    def find(self, f, *args):
        current = self.start

        for x0 in np.arange(self.start, self.stop + self.step, self.step):
            if x0 < current:
                continue
            x = self.find_root(f, x0, *args)
            if x is None:  # no root found.
                continue
            current = x
            self.add_to_roots(x)

        return self.roots

    def find_root(self, f, x0, *args):

        x, _, ier, _ = fsolve(f, x0=x0, args=args, full_output=True, xtol=self.xtol)
        if ier == 1:
            return x[0]
        return None

测试功能

scipy.special.jnjn不再存在，但我为案例创建了类似的测试功能。

def f(u, V=90, ell=5):
    w = np.sqrt(V ** 2 - u ** 2)

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

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

score 0 · Accepted Answer

我也遇到过这个问题来解决像 f(z)=0 这样的方程，其中 f 是一个全纯函数。我想确保不会错过任何零，最后开发了一种基于论证原理的算法。

它有助于找到复杂域中零的确切数量。一旦知道零的数量，就更容易找到它们。然而，有两个问题必须考虑：

注意多重性：当求解 (z-1)^2 = 0 时，你会得到两个零，因为 z=1 计数两次
如果函数是亚纯函数（因此包含极点），则每个极点会减少零的数量并中断对它们进行计数的尝试。

python - 如何使用numpy（和scipy）找到函数的所有零？

4 回答 4

示例输出：

RootFinder 定义

测试功能

Related

Reference