Как реализовать следующую формулу для производных в Python?

Question

Я пытаюсь реализовать следующую формулу в Python для точек X и Y

Я попробовал следующий подход

def f(c):
    """This function computes the curvature of the leaf."""
    tt = c
    n = (tt[0]*tt[3] - tt[1]*tt[2])
    d = (tt[0]**2 + tt[1]**2)
    k = n/d
    R = 1/k # Radius of Curvature
    return R

Есть что-то неправильное, поскольку оно не дает мне правильного результата. Я думаю, что я делаю некоторую ошибку, вычисляя производные в первых двух строках. Как я могу это исправить?

Вот некоторые из точек, которые находятся в кадре данных:

pts = pd.DataFrame({'x': x, 'y': y})


           x            y
    0.089631    97.710199
    0.089831    97.904541
    0.090030    98.099313
    0.090229    98.294513
    0.090428    98.490142
    0.090627    98.686200
    0.090827    98.882687
    0.091026    99.079602
    0.091225    99.276947
    0.091424    99.474720
    0.091623    99.672922
    0.091822    99.871553
    0.092022    100.070613
    0.092221    100.270102
    0.092420    100.470020
    0.092619    100.670366
    0.092818    100.871142
    0.093017    101.072346
    0.093217    101.273979
    0.093416    101.476041
    0.093615    101.678532
    0.093814    101.881451
    0.094013    102.084800
    0.094213    102.288577


pts_x = np.gradient(x_c, t)  # first derivatives
pts_y = np.gradient(y_c, t)
pts_xx = np.gradient(pts_x, t)  # second derivatives
pts_yy = np.gradient(pts_y, t)

После получения производных я помещаю производные x_prim, x_prim_prim, y_prim, y_prim_prim в другой кадр данных, используя следующий код:

d = pd.DataFrame({'x_prim': pts_x, 'y_prim': pts_y, 'x_prim_prim': pts_xx, 'y_prim_prim':pts_yy})

после того, как все находится в кадре данных, я вызываю функцию для каждой строки кадра данных, чтобы получить кривизну в этой точке, используя следующий код:

# Getting the curvature at each point
for i in range(len(d)):
    temp = d.iloc[i]
    c_temp = f(temp)
    curv.append(c_temp)

Answer 1

Вы не указываете точно, какова структура параметра pts. Но кажется, что это двумерный массив, в котором каждая строка имеет два значения x и y, а строки - это точки на вашей кривой. Это само по себе проблематично, поскольку в документации не совсем ясно, что именно возвращается в таком случае.

Но вы явно не получаете производные от x или y. Если вы предоставляете только один массив для np.gradient, то numpy предполагает, что точки равномерно расположены на расстоянии один. Но это, вероятно, не тот случай. Значение x' в вашей формуле является производной от x по отношению к t, переменной параметра для кривой (которая отделена от параметров компьютерных функций). Но вы никогда не предоставляете значения t для numpy. Значения t должны быть вторым параметром, передаваемым в функцию gradient.

Чтобы получить производные, разделите значения x, y и t на отдельные одномерные массивы - давайте назовем их x и y и t. Затем получите свои первые и вторые производные с

pts_x = np.gradient(x, t)  # first derivatives
pts_y = np.gradient(y, t)
pts_xx = np.gradient(pts_x, t)  # second derivatives
pts_yy = np.gradient(pts_y, t)

Тогда продолжайте оттуда. Вам больше не нужны значения t для расчета кривизны, которая является точкой формулы, которую вы используете. Обратите внимание, что gradient на самом деле не предназначен для вычисления вторых производных, и его абсолютно не следует использовать для вычисления производных третьего или более высокого порядка. Для них нужны более сложные формулы. Numpy's gradient использует «точные центральные различия второго порядка», которые довольно хороши для первой производной, плохи для второй производной и бесполезны для производных высшего порядка.

Answer 2

Я думаю, что ваша проблема в том, что x и y являются массивами двойных значений.

Массив x является независимой переменной; Я ожидаю, что это будет отсортировано в порядке возрастания. Если я оценю y [i], я ожидаю получить значение кривой в точке x [i].

Когда вы вызываете эту функцию numpy, вы получаете массив производных значений, которые имеют ту же форму, что и массивы (x, y). Если есть n пар из (x, y), то

y'[i] gives the value of the first derivative of y w.r.t. x at x[i]; 

y''[i] gives the value of the second derivative of y w.r.t. x at x[i].

Кривизна k также будет массивом с n точками:

k[i] = abs(x'[i]*y''[i] -y'[i]*x''[i])/(x'[i]**2 + y'[i]**2)**1.5

Думайте о x и y как о функциях параметра t. x '= dx / dt и т. д. Это означает, что кривизна k также является функцией этого параметра t.

Мне нравится, когда у меня есть понятное решение в закрытой форме, когда я программирую решение.

y(x) = sin(x) for 0 <= x <= pi
y'(x) = cos(x)
y''(x) = -sin(x)
k = sin(x)/(1+(cos(x))**2)**1.5

Теперь у вас есть хорошая формула для кривизны как функции от x.

Если вы хотите параметризовать его, используйте

x(t) = pi*t for 0 <= t <= 1
x'(t) = pi
x''(t) = 0

Посмотрите, сможете ли вы построить их и сделать так, чтобы ваше решение Python соответствовало им.