曲率的数值计算

Question

曲率的数值计算

6

我希望计算每个点的局部曲率。我有一组数据点，这些点在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）的比较（实线绿色）。

为什么两者之间会有如此大的偏差？如何获得准确的解析值。

希望得到帮助。

- newstudent

你说的“x轴上的点等距离”是什么意思？你所提供的方程式描述的是一个平面的曲率，该平面由x(t), y(t)描述。如果x是你的自变量，那么dx/dx==1。你是指在t或其他自变量上等距离分布吗？ - FHTMitchell

3

没有数据的话就无法确定发生了什么（我猜你看到的是由于二次求导而被放大的噪声）。分析性的东西是从哪里来的？ - MB-F

@Jonas已更新。 - newstudent

@kazemakase 分析图来自拟合。同样的图形，所以需要获得相同的曲率。 - newstudent

我用10000个点做了同样的实验，曲率看起来很好（像预期的那样）。我仍在寻找为什么200个点不够的原因。 - Thibault D.

显示剩余8条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Thibault D. · Accepted Answer

我尝试调整您的数值并发现它们不够平滑以计算曲率。事实上，即使是一阶导数也存在缺陷。原因如下：

Few plots of given data and my interpolation

您可以看到，蓝色的数据看起来像一个抛物线，它的导数应该像一条直线，但实际上并非如此。当您进行二阶导数时，情况会变得更糟。红色的图形是使用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()

我的建议是使用抛物线对数据进行插值，然后计算出尽可能多的插值点以供使用。