curvatu的数值计算

Question

我想计算局部曲率，即每个点的曲率。我有一组数据点，在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）（纯绿色）的比较

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

感谢帮助。在

Answer 1

我一直在玩你的值，我发现它们不够平滑，无法计算曲率。事实上，即使是一阶导数也有缺陷。 原因如下：

你可以看到蓝色的你的数据看起来像一条抛物线，它的导数应该看起来像一条直线，但事实并非如此。当你取二阶导数时，情况会变得更糟。在红色中，这是一条平滑的抛物线，用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()

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