曲率数值计算
曲率数值计算
提问于 2018-05-30 12:07:17
我想计算一下局部曲率，也就是每一点。我有一组数据点，在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)(固体绿色)的比较。
为什么两者之间有这么大的偏差？如何得到解析值的精确值。
帮助感激。
回答 1
Stack Overflow用户
回答已采纳
发布于 2018-05-30 13:35:13
我一直在玩你的值，我发现它们不够光滑，不足以计算曲率。事实上，即使是一阶导数也有缺陷。原因如下：
你可以看到蓝色你的数据看起来像一个抛物线，它的导数应该看起来像一条直线，但它不是。当你取二阶导数时，情况会变得更糟。在红色，这是一个光滑的抛物线计算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()我的建议是用抛物线插值你的数据，并计算出在这个插值中要处理的点数。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/50604298
复制相关文章
相似问题
腾讯云开发者
