Numerical calculation of curvature

Question

I would like to calculate local curvature i.e at each point. I have a set of data points that a equally spaced in x. Below is the code for generating curvature.

data=np.loadtxt('newsorted.txt') #data with uniform spacing
x=data[:,0]
y=data[:,1]

dx = np.gradient(data[:,0]) # first derivatives
dy = np.gradient(data[:,1])

d2x = np.gradient(dx) #second derivatives
d2y = np.gradient(dy)

cur = np.abs(d2y)/(1 + dy**2))**1.5 #curvature

Below is the image of curvature (magenta) and its comparison with analytical (equation: -0.02*(x-500)**2 + 250)(solid green)

Why does there is so much deviation between the two? How to get exact values that of analytical.

Help appreciated.

What do you mean the points are equally spaced in `x`? The equation you have is for the curvature of a plane described as `x(t), y(t)`. If `x` is your independent variable then `dx/dx == 1`. Do you mean equally spaced in `t` or whatever your independent variable is? — FHTMitchell, May 30 '18 at 12:16
Might be helpful: https://stackoverflow.com/a/28270382/826983 — Stefan Falk, May 30 '18 at 12:23
what does your data look like? A `plot(x,y)` would be helpful — Jonas, May 30 '18 at 12:35
@newstudent Without seeing data it is impossible to tell what's going on (I guess you see noise amplified by deriving twice). Where does the analytical thing come from? — MB-F, May 30 '18 at 12:36
@kazemakase analytical plot comes from a fit. same plot. so..need to get same curvature. — newstudent, May 30 '18 at 12:38
Hmm I think you're going to need to give us the data to help you. At least a sample. — FHTMitchell, May 30 '18 at 12:46
I have a question that might help : does increasing the number of points in your data improves the result ? Your purple line is very symetric and i feel like your data is undersampled. — Thibault D., May 30 '18 at 12:51
@pLOPeGG Actually my original data had non uniform spacing (also symmetric - its a parabola), hence reduced the points. Also linked in the post. — newstudent, May 30 '18 at 12:55
I made the same experiment with 10000 points, the curvature looks fine (like expected). I'm still searchin why 200 points aren't enough — Thibault D., May 30 '18 at 13:10
Thanks. Have you used the data ? Is the data good? or too disoriented? — newstudent, May 30 '18 at 13:13
@newstudent The curvature is different because the data is different. But this should be obvious, what do you want to know exactly? Do you think the differences in curvature should be smaller? Why? — Stop harming Monica, May 30 '18 at 19:45

score 3 · Accepted Answer · answered May 30 '18 at 13:35

I've been playing a bit with your values and i've found that they aren't smooth enough to compute curvature. In fact, even the first derivative is flawed. Here is why :

You can see in blue your data looks like a parabola, and it's derivative should look like a straight line but it does not. And it get worse when you take the second derivative. In red this is a smooth parabola computed with 10000 points (tried with 100 points, it works the same : perfect lines and curvature). I made a little script to 'enrich' your data, increasing artificially the number of points but it only get worse, here is my script if you want to try.

import numpy as np
import matplotlib.pyplot as plt

def enrich(x, y):
    x2 = []
    y2 = []
    for i in range(len(x)-1):
        x2 += [x[i], (x[i] + x[i+1]) / 2]
        y2 += [y[i], (y[i] + y[i + 1]) / 2]
    x2 += [x[-1]]
    y2 += [y[-1]]
    assert len(x2) == len(y2)
    return x2, y2

data = np.loadtxt('newsorted.txt')
x = data[:, 0]
y = data[:, 1]

for _ in range(0):
    x, y = enrich(x, y)

dx = np.gradient(x, x)  # first derivatives
dy = np.gradient(y, x)

d2x = np.gradient(dx, x)  # second derivatives
d2y = np.gradient(dy, x)

cur = np.abs(d2y) / (np.sqrt(1 + dy ** 2)) ** 1.5  # curvature


# My interpolation with a lot of points made quickly
x2 = np.linspace(400, 600, num=100)
y2 = -0.0225*(x2 - 500)**2 + 250

dy2 = np.gradient(y2, x2)

d2y2 = np.gradient(dy2, x2)

cur2 = np.abs(d2y2) / (np.sqrt(1 + dy2 ** 2)) ** 1.5  # curvature

plt.figure(1)

plt.subplot(221)
plt.plot(x, y, 'b', x2, y2, 'r')
plt.legend(['new sorted values', 'My interpolation values'])
plt.title('y=f(x)')
plt.subplot(222)
plt.plot(x, cur, 'b', x2, cur2, 'r')
plt.legend(['new sorted values', 'My interpolation values'])
plt.title('curvature')
plt.subplot(223)
plt.plot(x, dy, 'b', x2, dy2, 'r')
plt.legend(['new sorted values', 'My interpolation values'])
plt.title('dy/dx')
plt.subplot(224)
plt.plot(x, d2y, 'b', x2, d2y2, 'r')
plt.legend(['new sorted values', 'My interpolation values'])
plt.title('d2y/dx2')

plt.show()

My recommendation would be to interpolate your data with a parabola and compute as many points on this interpolation to work with.

Thanks a lot for the answer. I got the same for dy/dx. Actually can you tell me what does "not smooth" mean ? also..If I run your code I get completely diff plots for everything. — newstudent, May 30 '18 at 13:52
That's weird if you don't get the same plots, how different are they ? And about the "smoothness", it refers to the differentiability of your function. In your case, your function does not look like differentiable when we plot the differential.That's why I suggest you to interpolate you function with a real parabola. — Thibault D., May 30 '18 at 14:05
Your Y axis is different on 3 of the plots, did you paste my code into a new python file ? Or maybe it's python / matplotlib version that differs ? — Thibault D., May 30 '18 at 14:19
That might be the reason, I'm working with python 3.6, matplotlib 2.2.2 and numpy 1.14.3. Try first to use python 3, and if it not enough look for upgrading numpy and matplotlib too. It should solve the problem — Thibault D., May 30 '18 at 14:29

Numerical calculation of curvature

1 Answers1