How to implement the following formula for derivatives in python?

Question

I'm trying to implement the following formula in python for X and Y points

I have tried following approach

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

There is something incorrect as it is not giving me correct result. I think I'm making some mistake while computing derivatives in first two lines. How can I fix that?

Here are some of the points which are in a data frame:

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)

After getting the derivatives I am putting the derivatives x_prim, x_prim_prim, y_prim, y_prim_prim in another dataframe using the following code:

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

after having everything in the data frame I am calling function for each row of the data frame to get curvature at that point using following code:

# 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

You do not specify exactly what the structure of the parameter pts is. But it seems that it is a two-dimensional array where each row has two values x and y and the rows are the points in your curve. That itself is problematic, since the documentation is not quite clear on what exactly is returned in such a case.

But you clearly are not getting the derivatives of x or y. If you supply only one array to np.gradient then numpy assumes that the points are evenly spaced with a distance of one. But that is probably not the case. The meaning of x' in your formula is the derivative of x with respect to t, the parameter variable for the curve (which is separate from the parameters to the computer functions). But you never supply the values of t to numpy. The values of t must be the second parameter passed to the gradient function.

So to get your derivatives, split the x, y, and t values into separate one-dimensional arrays--lets call them x and y and t. Then get your first and second derivatives with

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)

Then continue from there. You no longer need the t values to calculate the curvatures, which is the point of the formula you are using. Note that gradient is not really designed to calculate the second derivatives, and it absolutely should not be used to calculate third or higher-order derivatives. More complex formulas are needed for those. Numpy's gradient uses "second order accurate central differences" which are pretty good for the first derivative, poor for the second derivative, and worthless for higher-order derivatives.

Answer 2

I think your problem is that x and y are arrays of double values.

The array x is the independent variable; I'd expect it to be sorted into ascending order. If I evaluate y[i], I expect to get the value of the curve at x[i].

When you call that numpy function you get an array of derivative values that are the same shape as the (x, y) arrays. If there are n pairs from (x, y), then

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].

The curvature k will also be an array with n points:

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

Think of x and y as both being functions of a parameter t. x' = dx/dt, etc. This means curvature k is also a function of that parameter t.

I like to have a well understood closed form solution available when I program a solution.

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

Now you have a nice formula for curvature as a function of x.

If you want to parameterize it, use

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

See if you can plot those and make your Python solution match it.