曲率的数值计算

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

我的建议是使用抛物线对数据进行插值，并在此插值上计算尽可能多的点。