Как сравнить переменные с плавающей запятой в наномасштабе

Вам действительно не нужно «грубо» отвечать, вы можете сделать это с помощью математики.

Рассмотрим простую синусоидальную функцию y = sin(2 pi * (x - x0) / wl + p0) с начальными условиями p0 = 115 градусов в x0 = 0 и длиной волны wl = 2

from matplotlib import pyplot as plt
import numpy as np


wl = 2
p0 = 115 * np.pi/180
x0 = 0

def func(x):
    return np.sin(2 * np.pi * (x - x0) / wl + p0)


x = np.linspace(x0, 10, 10000)
y = func(x)

Фаза меняется на 2*pi каждый раз, когда мы перемещаем wl единиц на x. Чтобы выяснить, насколько мы должны переместиться в x, чтобы изменить фазу с p0 на 2*pi, нам нужно решить уравнение:

 wl            ?
----  =  -------------
2*pi      (2*pi - p0)

который дает:

? = wl * (2*pi - p0) / (2*pi)

Фаза этой волны 2*pi впервые в x = x0 + wl * (2pi - p0) / 2pi. Чтобы получить последующие местоположения, где фаза равна 2*pi, нам просто нужно перейти вперед на целое число, кратное длине волны. Общая формула будет такой:

xphase_2pi = x0 + wl * (n + (2*pi - p0) / (2*pi))

Реализуя это, мы видим:

^{Note: I used n = np.arange(0, 5) to simultaneously calculate 5 points where the phase is 2*pi}

xphase_2pi = x0 + wl * (np.arange(0, 5) + (2*np.pi - p0) / (2 * np.pi))
# array([1.36111111, 3.36111111, 5.36111111, 7.36111111, 9.36111111])

yphase_2pi = func(xphase_2pi)
# array([-2.44929360e-16, -4.89858720e-16, -7.34788079e-16, -9.79717439e-16,
       -1.22464680e-15])

Как и ожидалось, y фактически равно нулю, когда фаза синусоиды равна 2*pi.

Мы также можем построить эти значения, чтобы убедиться:

fig = plt.figure()
plt.plot(x, y)
plt.plot(x0, func(x0), 'ok')
plt.plot(xphase_2pi, yphase_2pi, 'xr')
plt.tight_layout()
plt.grid(True)

График показывает, что все красные символы x находятся в правильном месте: в нуле, когда волна поднимается.

Если вы настаиваете на том, чтобы сделать это численно, вы можете использовать один из многих удивительных алгоритмов нахождения нуля , с которыми поставляется scipy . Вам действительно не нужно строить свой собственный. Например, вы можете использовать scipy.optimize.fsolve, чтобы решить уравнение, начиная с вашего начального местоположения. Если вы найдете ноль там, где наклон увеличивается, у вас есть ответ! Если вы найдете ноль там, где наклон падает, добавьте к нему wl/2, и вы получите ответ!

>>> x_zero = scipy.optimize.fsolve(func, x0)
# array([0.36111111])
>>> xphase_2pi = x_zero + np.arange(0, 5) * wl
>>> delta_x = wl / 10000
... slope = (func(x_zero + delta_x) - func(x_zero - delta_x)) / (delta_x * 2)
# array([-3.14159245])
>>> if slope[0] < 0:
...    xphase_2pi += wl / 2

Что приводит к тому же xphase_2pi:

array([1.36111111, 3.36111111, 5.36111111, 7.36111111, 9.36111111])

Если вы все еще хотите построить свой собственный, Метод Ньютона — это классный алгоритм поиска нуля, который вы можете реализовать. Для удобства мы можем определить func_slope, функцию, которая будет возвращать наклон волны при заданном x. Наклон рассчитывается как центральная разность по 1/5000 длины волны, как мы делали выше:

def func_slope(func, x):
    h = wl / 10000
    return (func(x + h) - func(x - h))/ (2 * h)

def newton(func, x0, max_iter=100):
    i = 0
    x = x0
    fx = func(x)
    fdx = func_slope(func, x)
    while not np.isclose(fx, 0) and i < max_iter:
        x = x - fx / fdx
        i += 1
        fx = func(x)
        fdx = func_slope(func, x)

    return x, fx, fdx, i

Теперь вы можете вызвать эту функцию:

>>> x_zero, err, slope, iters = newton(func, x0)
# (0.36111111111454486, -1.0787248531787957e-11, -3.1415924468797582, 5)

>>> xphase_2pi = x_zero + np.arange(0, 5) * wl + (slope < 0) * wl / 2
# array([1.36111111, 3.36111111, 5.36111111, 7.36111111, 9.36111111])