Stabilny numerycznie sposób obliczania kątów między wektorami

Podczas stosowania klasycznej formuły kąta między dwoma wektorami:

α = \arccos \frac{v_{1} \cdot v_{2}}{‖ v_{1} ‖ ‖ v_{2} ‖}

$\alpha = \arccos \frac{\mathbf{v_1} \cdot \mathbf{v_2}}{\|\mathbf{v_1}\| \|\mathbf{v_2}\|}$

stwierdzono, że dla bardzo małych / ostrych kątów występuje utrata precyzji, a wynik nie jest dokładny. Jak wyjaśniono w tej odpowiedzi Przepełnienie stosu , jednym rozwiązaniem jest użycie arcus tangens zamiast:

α = \arctan 2 (‖ v_{1} \times v_{2} ‖, v_{1} \cdot v_{2})

$\alpha = \arctan2 \left(\|\mathbf{v_1} \times \mathbf{v_2}\|, \mathbf{v_1} \cdot \mathbf{v_2} \right)$

I to rzeczywiście daje lepsze wyniki. Zastanawiam się jednak, czy dałoby to złe wyniki dla kątów bardzo zbliżonych do $\pi / 2$ . Czy to prawda? Jeśli tak, to czy istnieje wzór umożliwiający dokładne obliczenie kątów bez sprawdzania tolerancji wewnątrz ifgałęzi?

algorithms precision numerical-limitations

— astrojuanlu
źródło

Będzie to zależeć od implementacji funkcji odwrotnej stycznej do dwóch parametrów. Wolne, stabilne wersje warunkowo przełączają się między pracą z x / y i y / x, aby zachować precyzję, podczas gdy te szybkie trzymają rzeczy we właściwej ćwiartce, a zatem nie są bardziej precyzyjne niż wersja z jednym parametrem.

— origimbo

Powinieneś zdefiniować „utratę precyzji”: załóżmy, że poprawną odpowiedzią jest

α

$\alpha$ a zamiast tego otrzymujesz

α + Δ

$\alpha + \Delta$ . Czy potrzebujesz

Δ ≪ α

$\Delta \ll \alpha$ czy

Δ ≪ π

$\Delta \ll \pi$ jest wystarczający?

— Stefano M

W tym przypadku poprawną odpowiedzią było a dostałem , oba .

α

$\alpha$

α \cdot 10^{8}

$\alpha \cdot 10^8$

≪ 1

$\ll 1$

— astrojuanlu

Odpowiedzi:

( Testowałem już to podejście i pamiętam, że zadziałało poprawnie, ale nie przetestowałem go specjalnie na to pytanie ).

O ile wiem, zarówno i . $\|\mathbf{v}_1\times \mathbf{v}_2\|$ $\mathbf{v}_1\cdot \mathbf{v}_2$

Zacznij od przeformułowania problemu jako znalezienia kąta trójkąta o długości boków,oraz(wszystkie są dokładnie obliczone w arytmetyce zmiennoprzecinkowej). Jest dobrze znany wariant Wzór Herona powodu Kahan ( miscalculating Obszaru oraz kąty igłą jak Triangle ), co pozwala na obliczenie powierzchni i kąt (między i ) trójkąta określonego przez jego długości bocznych, i zrób to stabilnie numerycznie. Ponieważ redukcja do tego podproblemu jest również dokładna, takie podejście powinno działać w przypadku dowolnych danych wejściowych. $a=|\mathbf{v}_1|$ $b=|\mathbf{v}_2|$ $c=|\mathbf{v}_1-\mathbf{v}_2|$ $a$ $b$

Cytując z tego artykułu (patrz str. 3), zakładając, że , Wszystkie nawiasy tutaj są umieszczone ostrożnie i mają one znaczenie; jeśli okaże się, że bierzesz pierwiastek kwadratowy z liczby ujemnej, wejściowe długości boków nie są długościami boków trójkąta. $a\geq b$

μ = {\begin{cases} c - (a - b), & if b \geq c \geq 0, \\ b - (a - c), & if c > b \geq 0, \\ invalid triangle, & otherwise \end{cases}

$\mu = \begin{cases} c-(a-b),&\text{if }b\geq c\geq 0,\\ b-(a-c),&\text{if }c>b\geq 0,\\ \text{invalid triangle},&\text{otherwise} \end{cases}$

a n g l e = 2 \arctan (\sqrt{\frac{((a - b) + c) μ}{(a + (b + c)) ((a - c) + b)}})

$\mathrm{angle} = 2\arctan\left( \sqrt{\frac{((a-b)+c)\mu}{(a+(b+c))((a-c)+b)}} \right)$

W pracy Kahana wyjaśniono, jak to działa, w tym przykłady wartości, dla których inne formuły zawodzą. Twoja pierwsza formuła dla to na stronie 4. $\alpha$ $C''$

Głównym powodem, dla którego sugeruję formułę Herona Kahana, jest to, że tworzy ona bardzo przyjemną prymitywną - wiele potencjalnie trudnych pytań dotyczących geometrii płaskiej można zredukować do znalezienia pola / kąta dowolnego trójkąta, więc jeśli możesz zredukować swój problem do tego, istnieje niezła, stabilna formuła i nie ma potrzeby wymyślać czegoś na własną rękę.

Edytuj Po komentarzu Stefano, stworzyłem wykres błędu względnego dla , ( kod ). Dwie linie to błędy względne dla i , przebiegające wzdłuż osi poziomej. Wygląda na to, że działa. $v_1=(1,0)$ $v_2=(\cos\theta, \sin\theta)$ $\theta=\epsilon$ $\theta=\pi/2-\epsilon$ $\epsilon$

— Cyryl
źródło

Dzięki za link i odpowiedź! Niestety druga formuła, którą napisałem, nie pojawia się w artykule. Z drugiej strony ta metoda może być nieco złożona, ponieważ wymaga projekcji w 2D.

— astrojuanlu

@astrojuanlu Nie ma tutaj rzutowania na 2d: cokolwiek to są dwa wektory 3d, definiują one pojedynczy (płaski) trójkąt między nimi - wystarczy znać jego długości boków.

— Kirill,

Masz rację, mój komentarz nie ma sensu. Myślałem o współrzędnych zamiast długości. Dzięki jeszcze raz!

— astrojuanlu,

@astrojuanlu Jeszcze jedna rzecz, na którą chcę zwrócić uwagę: wydaje się, że istnieje formalny dowód na to, że formuła pola jest dokładna w artykule Jak obliczyć powierzchnię trójkąta: Formalna wizyta , Sylvie Boldo , przy użyciu Flocq.

— Kirill,

Doskonała odpowiedź, ale zaprzeczam, że zawsze można dokładnie obliczyć w arytmetyki zmiennoprzecinkowej. W rzeczywistości, jeśli wówczas dochodzi do katastrofalnych odwołań przy obliczaniu składników .

c

$c$

c < ϵ \cdot min (a, b)

$c < \epsilon \cdot \min(a,b)$

(v_{1} - v_{2})

$(\mathbf{v}_1-\mathbf{v}_2)$

— Stefano M

Skuteczną odpowiedź na to pytanie nie jest zaskakujące w innej notatce Velvela Kahana :

α = 2 \arctan (‖ \frac{v_{1}}{‖ v_{1} ‖} + \frac{v_{2}}{‖ v_{2} ‖} ‖, ‖ \frac{v_{1}}{‖ v_{1} ‖} - \frac{v_{2}}{‖ v_{2} ‖} ‖)

$\alpha=2\arctan\left(\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}+\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|,\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}-\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|\right)$

gdzie używam jako kąta utworzonego przez z osią poziomą. (W niektórych językach może być konieczne odwrócenie kolejności argumentów). $\arctan(x,y)$ $(x,y)$

(Dałem Mathematica demonstrację wzoru Kahan jest tutaj ).

— JM
źródło

Masz na myśli ?

\arctan 2

$\arctan2$

— astrojuanlu

Jestem przyzwyczajony do przedstawiania dwuargumentowego arcus tangens jako , tak. W języku takim jak FORTRAN byłby to odpowiednik .

\arctan (x, y)

$\arctan(x,y)$ ATAN2(Y, X)

— JM