Мне нужно найти пересечения нуля в одномерном массиве функции с периодичностью c. Это будут точки, где орбитальный спутник пересекает экватор Земли в направлении на север.
Я разработал простое решение, основанное на поиске точек, в которых одно значение равно нулю или отрицательно, а другое положительно, а затем с помощью квадратичный c или кубический c интерполятор с scipy.optimize.brentq, чтобы найти близлежащие нули.
Интерполятор не go выходит за кубы c, и прежде чем я научусь использовать лучший интерполятор, сначала я хотел бы проверить, существует ли уже быстрый метод в numpy или scipy, чтобы найти все пересечения нуля в большом массиве (n = 1E + 06 до 1E + 09).
Вопрос: Поэтому я спрашиваю, существует ли более быстрый метод в numpy или scipy, чтобы найти все пересечения нуля в большом массиве (n = 1E + 06 до 1E + 09), чем то, как я это сделал здесь?
График показывает ошибки между интерполированными нулями и действительное значение функции: чем меньше строка интерполяции в кубических c, тем больше квадратичные c.
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()