Блог

Я хотел бы вычислить градиент сплайна и визуализировать его на полярном графике.

Я использую tck, u = splprep() для получения сплайна с декартовыми координатами и splev(u, tck, der=1) для вычисления его частных производных в направлении x и y соответственно. Затем я вычисляю конечные точки стрелок, которые должны визуализировать градиент, и преобразую их в полярные координаты.

На первый взгляд график выглядит хорошо, но если я сравню расчетное направление градиента с аналитическим решением, будут существенные различия, даже если я увеличу количество точек.

MWE

 from matplotlib import pyplot as plt
import numpy as np
from scipy.interpolate import splprep, splev

if __name__ == '__main__':
    N = 11  # number of samples
    # x = np.arange(0, N)  # [Update] This did not decrease the step size for increasing N
    x = np.linspace(0, 10, N)
    y = np.sin(x)

    tck, u = splprep([x, y], s=0)  # spline

    theta, r = np.arctan2(y, x), np.hypot(x, y)  # convert to polar

    gradient = splev(u, tck, der=1) # compute first derivative

    # normalize
    gradient = gradient / np.hypot(gradient[0], gradient[1])

    # compare numerical and analytical solution
    # direction = np.arctan2(gradient[1], gradient[0])  # [Update] this was wrong
    slope = gradient[1] / gradient[0]
    print(np.cos(x) - slope)  # cos(x) should be the analytical solution

    endpoints_x = x   gradient[0]
    endpoints_y = y   gradient[1]

    # convert cartesian endpoints to polar
    endpoints_theta, endpoints_r = np.arctan2(endpoints_y, endpoints_x),
                                   np.hypot(endpoints_x, endpoints_y)


    fig, ax = plt.subplots(1, 1, subplot_kw=dict(polar=True))

    plt.scatter(theta, r, marker='o')

    plt.plot(np.stack((theta, endpoints_theta)), np.stack((r, endpoints_r)), 'r')

    plt.show()
  

Скриншот

Обновить

Я обнаружил две ошибки. Во-первых, размер шага в x не уменьшался за счет увеличения, N как я использовал np.arange(0, N) . Во-вторых, я ожидал, что численное решение будет arctan2(gradient[1], gradient[0]) но это просто наклон gradient[1] / gradient[0] в декартовых координатах.

Теперь все работает так, как ожидалось.