Численно устойчивый способ вычисления углов между векторами

14

При применении классической формулы для угла между двумя векторами:

α=arccosv1v2v1v2

обнаруживается, что при очень малых / острых углах наблюдается потеря точности, и результат не является точным. Как объясняется в ответе на переполнение стека , одно из решений - использовать вместо этого арктангенс:

α=arctan2(v1×v2,v1v2)

И это действительно дает лучшие результаты. Тем не менее, мне интересно, даст ли это плохие результаты для углов, очень близких к π/2 . Это так? Если да, то есть ли какая-нибудь формула для точного вычисления углов без проверки допуска внутри ifветви?

astrojuanlu
источник
1
Это будет зависеть от реализации двухпараметрической обратной касательной функции. Медленные, стабильные версии условно переключаются между работой с x / y и y / x для поддержания точности, в то время как быстрые версии просто вставляют вещи в правильный квадрант и, следовательно, не более точны, чем версия с одним параметром.
origimbo
Вы должны определить «потерю точности»: предположим, что правильный ответ α а вместо него вы получите α+Δ . Вам нужен Δα или достаточно Δπ ?
Стефано М
В этом случае правильный ответ был и я получил , оба . αα1081
astrojuanlu

Ответы:

18

( Я проверял этот подход раньше, и я помню, что он работал правильно, но я не проверял его специально для этого вопроса. )

Насколько я могу судить, оба и может пострадать от катастрофической отмены, если они почти параллельны / перпендикулярны - atan2 не может дать вам хорошую точность, если любой из входов отключен.v1×v2v1v2

Начните с переформулировки задачи как нахождения угла треугольника с длинами сторон,и(все они точно рассчитаны в арифметике с плавающей запятой). Существует хорошо известный вариант формулы Герона из-за Кахана ( Miscalculation Area and Angles of Needle-подобные треугольники ), который позволяет вычислять площадь и угол (между и ) треугольника, заданного длиной его стороны, и делать это численно стабильно. Поскольку приведение к этой подзадаче также является точным, этот подход должен работать для произвольных входных данных.a=|v1|b=|v2|c=|v1v2|ab

Цитирование из этой статьи (см. Стр.3), предполагая, что , Все круглые скобки здесь размещены аккуратно, и они имеют значение; если вы обнаружите, что берете квадратный корень из отрицательного числа, длина входной стороны не равна длине стороны треугольника.ab

μ={c(ab),if bc0,b(ac),if c>b0,invalid triangle,otherwise
angle=2arctan(((ab)+c)μ(a+(b+c))((ac)+b))

Существует объяснение того, как это работает, включая примеры значений, для которых другие формулы не работают, в статье Кахана. Ваша первая формула для - это на странице 4.αC

Основная причина, по которой я предлагаю формулу Керона Кахана, состоит в том, что она делает очень хороший примитив - множество потенциально хитрых вопросов плоской геометрии может быть сведено к нахождению площади / угла произвольного треугольника, так что если вы можете свести свою проблему к этому, хорошая стабильная формула для этого, и нет необходимости придумывать что-то самостоятельно.

Редактировать После комментария Стефано я сделал график относительной ошибки для , ( код ). Две линии - это относительные ошибки для и , идущих вдоль горизонтальной оси. Кажется, это работает. v1=(1,0)v2=(cosθ,sinθ)θ=ϵθ=π/2ϵϵвведите описание изображения здесь

Кирилл
источник
Спасибо за ссылку и ответ! К сожалению, вторая формула, которую я написал, не появляется в статье. С другой стороны, этот метод может быть немного сложным, так как требует проекции в 2D.
astrojuanlu
2
@astrojuanlu Здесь нет проекции на 2d: какими бы ни были эти два 3d-вектора, они определяют один (плоский) треугольник между ними - вам нужно знать только длину его сторон.
Кирилл
Вы правы, мой комментарий не имеет смысла. Я думал в координатах, а не в длине. Еще раз спасибо!
astrojuanlu
2
@astrojuanlu Еще одна вещь, которую я хочу отметить: кажется, что есть формальное доказательство того, что формула площади точна в разделе Как вычислить площадь треугольника: формальное пересмотр , Сильви Болдо , с использованием Flocq.
Кирилл
Отличный ответ, но я спорю, что вы всегда можете точно вычислить в арифметике плавающей запятой. Фактически, если то при вычислении компонентов происходят катастрофические сокращения . cc<ϵmin(a,b)(v1v2)
Стефано М
7

Эффективный ответ на этот вопрос, что неудивительно, в другой заметке Велвела Кахана :

α=2arctan(v1v1+v2v2,v1v1v2v2)

где я использую в качестве угла с горизонтальной осью. (Возможно, вам придется изменить порядок аргументов в некоторых языках.)arctan(x,y)(x,y)

(Я дал Mathematica демонстрацию формулы Кахана здесь .)

JM
источник
Вы имеете в виду ? arctan2
astrojuanlu
1
Я привык изображать арктангенс с двумя аргументами как , да. На таком языке, как FORTRAN, эквивалент будет . arctan(x,y)ATAN2(Y, X)
JM