Нахождение местоположения заданной дуги расстояния вдоль параметризованного кривой в Python

Question

У меня есть параметризованный 2D кривой: (х, у) = F (T)

функция Р произвольно, но дифференцируемы и, следовательно, смогу выведите длину дифференциальной дуги ds вдоль кривой в любой точке, используя стандартные формулы. Я также могу получить полную длину дуги S (t) от начала до любой точки на кривой, численно интегрируя формулу длины дифференциальной дуги. Я могу контролировать точность вычисления.

Я хочу найти точку (x, y), которая имеет полную длину дуги S = ​​D от начала кривой. Еще лучше, если реализация была в python. Я буду делать это много раз, и это часть вычислительного приложения, где мне нужен жесткий контроль точности и некоторая уверенность в конвергенции.

Я не знаю, подходит ли поиск корня, но мой вопрос является эквивалентом проблемы поиска корней для g (t) = S (t) - D, где функция g (t) не вычисляется точно так как S (t) нет. Оценка неточной функции беспорядочна не только с точностью, но и монотонностью g (t). С самого начала я старался делать плотную численную интеграцию, но навсегда. Я уверен, что сходится к моей требуемой толерантности, алгоритм поиска корней должен был бы лениво контролировать точность интеграции, когда он продолжался, требуя небрежной оценки с самого начала и повышения точности по мере сглаживания корневого алгоритма.

Есть ли такая вещь в наличии? Есть ли альтернативный умный способ сделать это?

Цените помощь

Answer 1

Вы можете разместить код, и сказать нам, что случилось с ним?

Вот моя версия, что вычислить т, где S (т) == D:

from scipy.integrate import quad 
from scipy.optimize import fsolve 
from math import cos, sin, sqrt, pi 

def circle_diff(t): 
    dx = -sin(t) 
    dy = cos(t) 
    return sqrt(dx*dx+dy*dy) 

def sin_diff(t): 
    dx = 1 
    dy = cos(t) 
    return sqrt(dx*dx+dy*dy) 

def curve_length(t0, S, length): 
    return quad(S, 0, t0)[0] - length 

def solve_t(curve_diff, length):  
    return fsolve(curve_length, 0.0, (curve_diff, length))[0] 

print solve_t(circle_diff, 2*pi) 
print solve_t(sin_diff, 7.640395578) 

Answer 2

OK, @HYRY, вот фрагмент кода в значительной степени основаны на вашей. Вы дали мне подсказку, которая мне нужна для успеха: используйте «квадрат» вместо «квадратурного». Поэтому я, по крайней мере, проголосую за ваш ответ, но я хочу добавить к этой истории.

Во-первых, ваш код работал быстро, но около пяти мест, за исключением точности, которой я был после. Я добавил квадтол и opttol к вашему примеру, чтобы попытаться проиллюстрировать взаимодействие квадратурной и корневой точности поиска. Я также добавил цикл, основанный на высоких допусках по умолчанию, чтобы выявить различия в скорости.

Пример греха гораздо более чувствителен, чем круг по точности. Я также добавил памеризированную кривую, длина дуги которой задана гипергеометрической функцией, а параметр «brentq» закомментирован, потому что fsolve терпит неудачу в этом примере и в любом четном brentq равен или лучше по этой задаче.

«квадратура» является медленной, но проявляет ожидаемое поведение: скорость поиска, точность и успех в корне с квадратурным допуском.

Напротив, «квадрат», кажется, игнорирует запрошенный допуск и дает более точный ответ всегда. Эта невыполненная точность будет раздражать или приглашать объяснения, за исключением того, что она также работает так быстро на примерах, что я не уверен, что мой вопрос интересен больше. Благодаря!

from scipy.integrate import quad, quadrature 
from scipy.optimize import fsolve, brentq 
from math import cos, sin, sqrt, pi, pow 

def circle_diff(t): 
    dx = -sin(t) 
    dy = cos(t) 
    return sqrt(dx*dx+dy*dy) 

def sin_diff(t): 
    dx = 1 
    dy = cos(t) 
    return sqrt(dx*dx+dy*dy) 

def hypergeom_diff(t): 
    """ y = t^5 x = t^3 """ 
    dx = 3*t*t 
    dy = 5*pow(t,4) 
    return sqrt(dx*dx+dy*dy) 

def curve_length(t0, S, length,quadtol): 
    integral = quad(S, 0, t0,epsabs=quadtol,epsrel=quadtol) 
    #integral = quadrature(S, 0, t0,tol=quadtol,rtol=quadtol, vec_func = False) 
    return integral[0] - length 

def solve_t(curve_diff, length,opttol=1.e-15,quadtol=1e-10): 
    return fsolve(curve_length, 0.0, (curve_diff, length,quadtol), xtol = opttol)[0] 
    #return brentq(curve_length, 0.0, 3.2*pi,(curve_diff, length, quadtol), rtol = opttol) 

for i in range(1000): 
    y = solve_t(circle_diff, 2*pi) 

print 2*pi 
print solve_t(circle_diff, 2*pi) 
print solve_t(sin_diff, 7.640395578) 
print solve_t(circle_diff, 2*pi,opttol=1e-5,quadtol=1e-3) 
print solve_t(sin_diff, 7.640395578,opttol = 1e-12,quadtol=1e-6) 
print "hypergeom" 
print solve_t(hypergeom_diff, 2.0,opttol = 1e-12,quadtol=1e-12) 
print solve_t(hypergeom_diff, 2.0,opttol = 1e-12,quadtol=1e-6)