Compute gradient of spline in polar coordinates

Viewed 188

I would like to compute the gradient of a spline and visualize it in a polar plot.

I use tck, u = splprep() to obtain the spline with the cartesian coordinates and splev(u, tck, der=1) to compute its partial derivatives in x and y direction, respectively. Then, I compute the endpoints of the arrows that are supposed to visualize the gradient and convert them to polar coordinates.

The plot looks good on the first glance but if I compare the estimated direction of the gradient to the analytical solution, there are significant differences even if I increase the number of points.

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

Screenshot

Gradient of a spline in polar coordinates

Update

I found two mistakes. First, the step size in x was not decreasing by increasing N as I used np.arange(0, N). Second, I expected the numerical solution to be arctan2(gradient[1], gradient[0]) but it is simply the slope gradient[1] / gradient[0] in cartesian coordinates.

Everything works as expected now.

0 Answers
Related