python - 曲率数值计算

我想计算局部曲率，即每个点的曲率。我有一组数据点，在x中等距分布。下面是生成曲率的代码。

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

下面是曲率图像（品红色）及其与解析曲线的比较（方程式：-0.02*(x-500)**2 + 250）（纯绿色）

为什么两者之间有这么大的偏差？如何得到精确的分析值。
帮助理解。

最佳答案

我对你的值做了一些修改，发现它们不够光滑，无法计算曲率。事实上，即使是一阶导数也有缺陷。
原因如下：
你可以用蓝色看到你的数据看起来像一条抛物线，它的导数应该看起来像一条直线，但它不是。如果你用二阶导数，情况会更糟。红色表示这是一条平滑的抛物线，计算10000个点（尝试100个点，其工作原理相同：完美的直线和曲率）。
我做了一个小脚本来“丰富”你的数据，人为地增加点数，但它只会变得更糟，这是我的脚本，如果你想试试。

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

我的建议是用一条抛物线来插值数据，并在这个插值上计算尽可能多的点。