Rozwiązanie równania kwartalnego

Czy istnieje otwarta implementacja C dla rozwiązania równań kwartalnych:

a x ⁴ + b x ³ + c x ² + d x + e = 0

$ax⁴+bx³+cx²+dx+e=0$

Mam na myśli wdrożenie rozwiązania Ferrari. Na Wikipedii czytam, że rozwiązanie jest stabilne obliczeniowo tylko dla niektórych możliwych kombinacji znaków współczynników. Ale może mam szczęście ... Mam pragmatyczne rozwiązanie, rozwiązując analitycznie za pomocą komputerowego systemu algebry i eksportując do C. Ale jeśli istnieje przetestowana implementacja, wolałbym to wykorzystać. Szukam szybkiej metody i wolę nie używać ogólnej wyszukiwarki rootów.

Potrzebuję tylko prawdziwych rozwiązań.

polynomials nonlinear-equations roots

— wysoki
źródło

Czy potrzebujesz wszystkich (prawdziwych) rozwiązań jednocześnie? Jak mówi GertVdE poniżej, jeśli masz problemy ze stabilnością w rozwiązaniu o zamkniętej formie, nie ma naprawdę dobrego powodu, aby nie używać algorytmu znajdowania root.

— Godric Seer

Zabawne jest dla mnie, że oznaczono to algebrą nieliniową, ponieważ można po prostu obliczyć wartości własne macierzy towarzyszącej, która jest już w formie Hessenberga, a zastosowanie wymiarów QR byłoby dość proste.

— Victor Liu

Rzuć okiem na solwery sześcienne / kwartalne opublikowane w ACM TOMS (Algorytm 954) . Kod, który trafia do tego dziennika, jest zazwyczaj bardzo wysokiej jakości. Sam papier znajduje się za zaporą, ale kod można pobrać z tego linku .

— GoHokies

... (późniejsza edycja) kod ACM jest napisany w FORTRAN 90, ale moje pierwsze wrażenie jest takie, że można go wywołać z C bez większego wysiłku.

— GoHokies

@ GoHokies Myślę, że powinieneś zamienić swój komentarz w odpowiedź, ponieważ uważam, że to dobra odpowiedź na to pytanie. Zwłaszcza, że połączony papier pozwala uniknąć zwykłych niestabilności numerycznych, a to absolutnie nie jest trywialna rzecz.

— Kirill

Odpowiedzi:

Zdecydowanie odradzam stosowanie rozwiązań w formie zamkniętej, ponieważ są one zwykle bardzo niestabilne numerycznie. Musisz zachować szczególną ostrożność w sposobie i kolejności ocen dyskryminatora i innych parametrów.

Klasycznym przykładem jest równanie kwadratowe . Obliczenie pierwiastków jako spowoduje kłopoty dla wielomianów, w których od tego czasu można anulować w licznik ułamka. Musisz obliczyć . $ax^2+bx+c=0$

x_{1, 2} = \frac{- b \pm \sqrt{b^{2} - 4 a c}}{2 a}

$x_{1,2} = \frac{-b \pm \sqrt{b^2-4ac}}{2a}$

b ≫ 4 a c

$b\gg 4ac$

x_{1} = \frac{- (b + s i g n (b) \sqrt{b^{2} - 4 a c})}{2 a}; x_{2} = \frac{c}{a} \frac{1}{x_{1}}

$x_1 = \frac{-(b+sign(b)\sqrt{b^2-4ac})}{2a}; x_2 = \frac{c}{a}\frac{1}{x_1}$

Higham w swoim arcydziele „Dokładność i stabilność algorytmów numerycznych” (wydanie 2, SIAM) wykorzystuje metodę bezpośredniego wyszukiwania, aby znaleźć współczynniki wielomianu sześciennego, dla których klasyczne analityczne rozwiązanie sześcienne daje bardzo niedokładne wyniki. Przykład, który podaje, to . W przypadku tego wielomianu korzenie są dobrze rozdzielone, a zatem problem nie jest źle uwarunkowany. Jeśli jednak wyliczy pierwiastki za pomocą metody analitycznej i oceni wielomian w tych korzeniach, otrzyma pozostałość , stosując metodę stabilnego standardu (metoda macierzy towarzyszącej) , pozostałość jest rzędu $[a,b,c] = [1.732, 1, 1.2704]$ $\mathcal{O}(10^{-2})$ $\mathcal{O}(10^{-15})$ . Proponuje niewielką modyfikację algorytmu, ale nawet wtedy może znaleźć zestaw współczynników prowadzących do reszt co zdecydowanie nie jest dobre. Zobacz str. 480–481 wyżej wymienionej książki. $\mathcal{O}(10^{-11})$

W twoim przypadku zastosowałbym metodę Bairstow . Wykorzystuje iteracyjną kombinację iteracji Newtona na formach kwadratowych (a następnie rozwiązane są pierwiastki kwadratowe) i deflacji. Można go łatwo wdrożyć, a w Internecie dostępne są nawet niektóre wdrożenia.

— GertVdE
źródło

Czy mógłbyś wyjaśnić, co masz na myśli przez „zdecydowanie odradzam stosowanie rozwiązań w formie zamkniętej, ponieważ są one zazwyczaj bardzo niestabilne numerycznie”. Czy dotyczy to wyłącznie wielomianów czwartego stopnia, czy jest to ogólna zasada?

— NoChance,

@EmmadKareem Zaktualizowałem swoją odpowiedź powyżej.

— GertVdE

Zobacz te:

Rozwiązywanie Quartics i Cubics dla Graphics , pierwotnie opublikowany w Graphics Gems V . Oryginalny kod jest tutaj . Zobacz także to i to .
Uniwersalna metoda rozwiązywania równań kwartalnych .

— lhf
źródło

Używając tego kodu na wielomianu ze współczynnikami podanymi w mojej odpowiedzi, znajduję, co następuje: , który ma względny błąd porównaniu do rzeczywistego katalogu głównego (obliczane przy użyciu polecenia root Octave, które korzysta z metody macierzy towarzyszących). Ma resztę podczas gdy root metody macierzy towarzyszącej ma resztę . Od Ciebie zależy, czy jest to wystarczająco dobre (w przypadku grafiki komputerowej może być, w przypadku innych aplikacji nie będzie)

x_{1} = - 1.602644912244132 e + 00

$x_1 = -1.602644912244132e+00$

O (10^{- 8})

$\mathcal{O}(10^{-8})$

O (10^{- 7})

$\mathcal{O}(10^{-7})$

O (10^{- 15})

$\mathcal{O}(10^{-15})$

— GertVdE

Receptury numeryczne w c zapewniają zamkniętą ekspresję formy dla prawdziwych pierwiastków kwadratowych i sześciennych, które przypuszczalnie mają przyzwoitą precyzję. Ponieważ algebraiczne rozwiązanie kwartyku polega na rozwiązaniu sześciennej, a następnie rozwiązaniu dwóch kwadratów, być może kwantyk o zamkniętej formie z dobrą precyzją nie jest wykluczony.

— Nemocopperfield
źródło

Właśnie dostałem pierwiastek sześciennego przykładu cytowanego w granicach 2e-16 (odrobinę ponad precyzję moich pływaków) przy użyciu receptur numerycznych w formułach sześciennych c (press et al). Jest więc powód do nadziei.

— Nemocopperfield