Численный расчет кривизны - PullRequest
Новая
       27

Численный расчет кривизны

0 голосов
/ 30 мая 2018

Я хотел бы рассчитать локальную кривизну, т.е. в каждой точке.У меня есть набор точек данных, которые одинаково разнесены по 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) (сплошной зеленый) curvature

Почему между ними так много различий?Как получить точные значения, аналитические.

Помощь приветствуется.

1 Ответ

0 голосов
/ 30 мая 2018

Я немного поиграл с вашими значениями и обнаружил, что они недостаточно гладкие, чтобы вычислять кривизну.На самом деле, даже первая производная имеет недостатки.Вот почему: 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()

Я бы порекомендовал интерполировать вашиданные с параболой и вычислите столько точек на этой интерполяции, чтобы работать с ними.

...