在大型准周期数组中找到所有正向零交叉点

Question

我需要在大致周期函数的一维数组中找到零交叉点。这将是轨道卫星向北穿过地球赤道的点。

我已经制定了一个简单的解决方案，基于查找一个值为零或负值而下一个值为正值的点，然后使用带有

scipy.optimize.brentq

的二次或三次插值器来查找附近的零点。

插值器不会超出三次，在我学习使用更好的插值器之前，我首先想检查 numpy 或 scipy 中是否已经存在一种快速方法来查找大数组中的所有零交叉（n = 1E+06 至 1E+09).

问题： 所以我问 numpy 或 scipy 中是否已经存在一种比我所做的方法更快的方法来查找大数组（n = 1E+06 到 1E+09）中的所有零交叉在这里？

该图显示了插值零点与函数实际值之间的误差，较小的线是三次插值，较大的是二次插值。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.optimize import brentq

def f(x):
    return np.sin(x + np.sin(x*e)/e) # roughly periodic function

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
e = np.exp(1)

x = np.arange(0, 10000, 0.1)
y = np.sin(x + np.sin(x*e)/e)
crossings = np.where((y[1:] > 0) * (y[:-1] <= 0))[0]

Qd = interp1d(x, y, kind='quadratic', assume_sorted=True)
Cu = interp1d(x, y, kind='cubic', assume_sorted=True)

x0sQd = [brentq(Qd, x[i-1], x[i+1]) for i in crossings[1:-1]]
x0sCu = [brentq(Cu, x[i-1], x[i+1]) for i in crossings[1:-1]]

y0sQd = [f(x0) for x0 in x0sQd]
y0sCu = [f(x0) for x0 in x0sCu]

if True:
    plt.figure()
    plt.plot(x0sQd, y0sQd)
    plt.plot(x0sCu, y0sCu)
    plt.show()

Answer 1

问题：所以我问 numpy 或 scipy 中是否已经存在一种比我在这里完成的方法更快的方法来查找大数组（n = 1E+06 到 1E+09）中的所有零交叉？

假设您的意思是“基于插值器近似函数的根”，是的。不要使用

interp1d

便利函数，而是使用插值器的面向对象接口（例如

CubicSpline

，并使用其

roots

方法来查找所有根。

运行代码后：

from scipy.interpolate import CubicSpline
spline = CubicSpline(x, y, extrapolate=False)

# find all the roots
x0sCu2 = spline.roots()

# extract the ones with positive slope
dspline = spline.derivative()
slope = dspline(x0sCu2)
x0sCu2 = x0sCu2[slope > 0]

# confirm that the results agree with `brent`
np.testing.assert_allclose(x0sCu2[1:-1], x0sCu)

这比你的速度要快一点。

%timeit [brentq(Cu, x[i-1], x[i+1]) for i in crossings[1:-1]]
# 866 ms ± 11.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit x0sCu2 = spline.roots(); dspline = spline.derivative(); slope = dspline(x0sCu2); x0sCu2 = x0sCu2[slope > 0]
# 264 ms ± 5.34 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

实际上有比内置

roots

方法更快的方法。您的想法是正确的，但您可以使用矢量化寻根器，例如

newton

。

from scipy.optimize import newton
x0sCu3 = newton(Cu, x[crossings[1:-1]])
np.testing.assert_allclose(x0sCu3, x0sCu)

# Ignores the time spent finding `crossings`, but that's
# negligible, and can be improved (e.g. we don't need `where`)
%timeit newton(Cu, x[crossings[1:-1]])
# 3.93 ms ± 563 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

它几乎肯定会收敛，因为你能够提供一个非常好的猜测。然而，

brentq

的好处是它基本上可以保证收敛，因为你提供了一个括号。有矢量化括号寻根器，但它是私有的，因此不要指望它在未来版本的 SciPy 中的同一位置可用。

from scipy.optimize._zeros_py import _chandrupatla
res = _chandrupatla(Cu, x[crossings[1:-1]-1], x[crossings[1:-1]+1])
np.testing.assert_allclose(res.x, x0sCu)

%timeit _chandrupatla(Cu, x[crossings[1:-1]-1], x[crossings[1:-1]+1])
# 7.54 ms ± 1.24 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)