Stabilne numerycznie jednoznaczne rozwiązanie małego układu liniowego

11

Mam niejednorodny układ liniowy

A x = b

$Ax=b$

gdzie jest rzeczywistą macierzą przy . Nullspace z jest gwarancją wymiaru zerowym więc równanie ma unikalną odwrotną . Ponieważ wynik wchodzi w prawą stronę ODE, którą zamierzam rozwiązać za pomocą metody adaptacyjnej, ważne jest, aby rozwiązanie było płynne w odniesieniu do niewielkich odmian elementów i . Z powodu tego wymogu i małych wymiarowości Myślałem wdrożyć wyraźne wzory na $A$ $n\times n$ $n\leq 4$ $A$ $x=A^{-1} b$ $A$ $b$ $A^{-1} b$ . Elementy mogą mieć dokładnie zero lub przyjmować bardzo różne wartości. Moje pytanie brzmi, czy ma to dla ciebie sens i czy są znane stabilne wyrażenia tego. Koduję w C dla systemów x86.

linear-algebra ode

— wysoki
źródło

Wiem, że przychodzi bardzo późno, ale oto moja sugestia: ponieważ wiadomo, że eliminacja Gaussa z całkowitym obrotem jest stabilna, sensowne może być zaprogramowanie algorytmu dla małych rozmiarów. Obracanie komplikuje materię, ponieważ istnieją

sposoby wybierania kolejnych osi, prowadząc do

różnych zestawów wzorów; możesz zmniejszyć tę złożoność, zamieniając to, co należy zamienić, zmniejszając liczbę przypadków do

.

(n!)^{2}

$(n!)^2$

(n!)^{2}

$(n!)^2$

1^{2} + 2^{2} + \dots n^{2}

$1^2+2^2+\cdots n^2$

— Yves Daoust,

6

Przed wdrożeniem wyraźnych formuł zadałbym sobie pytanie: „czy warto?”:

Czy warto poświęcić czas na pisanie, debugowanie i sprawdzanie tych wyraźnych formuł, podczas gdy można łatwo połączyć się z BLAS + LAPACK, które wykorzystują klasyczną eliminację Gaussa?
Czy spodziewasz się stabilności? Nie sądzę, aby programowanie jawnych formuł (takich jak reguła Cramera) zapewniło lepszą stabilność, wręcz przeciwnie.
Czy spodziewasz się przyspieszyć? Czy profilowałeś już cały swój program? Jaką część czasu spędza na rozwiązywaniu hese systemów 4x4?
Jakie jest prawdopodobieństwo, że za rok poprawisz swój model i potrzebujesz 5 równań zamiast 4?

Moja rada: najpierw użyj kombinacji BLAS / LAPACK, sprawdź, czy to działa, profiluj cały program, poproś ucznia o wdrożenie jawnych formuł (przepraszam, że jesteś tutaj sarkastyczny) i dokonaj porównania szybkości i niezawodności.

— GertVdE
źródło

Wysiłek, jaki zajmuje mi wdrożenie, to około 15 minut, ponieważ po prostu wpisuję ogólną macierz 1x1, 2x2, 3x3 i 4x4 do CAS (dla mnie Klon) i odwracam ją. Zwraca wyraźny (podobny do C) wynik (podobno oparty na zasadzie Cramera). Twój drugi punkt jest dokładnie moją troską. W rezultacie pojawią się produkty wyższego rzędu elementów macierzy. Oczywiście może to powodować błędy wynikające z „prawie anulowania” różnych warunków. Ale pytanie brzmi, czy można zapisać wynik w formie, w której tak się nie dzieje. Prędkość nie jest głównym problemem w tym miejscu.

— highsciguy

6

$\mathcal{O}(n^{3})$

$A$ $A$ $\det(A) \neq 0$ $x$ $b$ $x$ $A$

$A$

$n = 2$

Spodziewałbym się, że problemem z użyciem BLAS + LAPACK do przeprowadzenia GEPP w rozwiązaniu ODE byłoby dowolne skończone różnicowanie zastosowane w niejawnej metodzie ODE. Wiem, że ludzie rozwiązali programy liniowe w ramach oceny po prawej stronie, a ponieważ zrobili to naiwnie (po prostu podłączyli program liniowy do rozwiązania po prawej stronie, wywołując algorytm simpleksowy), znacznie zmniejszyli dokładność swoich obliczone rozwiązanie i znacznie wydłużyło czas potrzebny na rozwiązanie problemu. Mój kolega z laboratorium wymyślił, jak rozwiązać takie problemy w znacznie bardziej wydajny i dokładny sposób; Będę musiał sprawdzić, czy jego publikacja została już wydana. Możesz mieć podobny problem, niezależnie od tego, czy zdecydujesz się użyć GEPP, czy reguły Cramera.

Jeśli istnieje jakikolwiek sposób, aby obliczyć analityczną matrycę jakobską dla swojego problemu, możesz to zrobić, aby zaoszczędzić sobie liczbowych bólów głowy. Będzie to tańsze i prawdopodobnie dokładniejsze niż przybliżenie różnic skończonych. Wyrażenia dla pochodnej odwrotności macierzy można znaleźć tutaj, jeśli są potrzebne. Ocena pochodnej odwrotności macierzy wygląda na to, że zajęłoby to co najmniej dwa lub trzy liniowe rozwiązania układu, ale wszystkie miałyby tę samą macierz i różne prawe boki, więc nie byłoby znacznie droższe niż pojedynczy układ liniowy rozwiązać.

A jeśli istnieje jakikolwiek sposób, aby porównać obliczone rozwiązanie do rozwiązania ze znanymi wartościami parametrów, zrobiłbym to, abyś mógł zdiagnozować, czy napotkałeś którykolwiek z tych pułapek numerycznych.

— Geoff Oxberry
źródło

Kiedy piszesz tutaj gładko, masz na myśli, że jest on również gładki, gdy jest oceniany ze skończoną precyzją, tj. Stabilny (tak próbowałem powiedzieć). Zobacz także mój komentarz do odpowiedzi GertVdE. Myślę, że mogę wykluczyć prawie pojedyncze macierze (przypuszczam, że w takich przypadkach analiza mojego problemu fizycznego musi zostać przeformułowana).

— highsciguy

1

A

$A$

det (A) \neq 0

$\det(A) \neq 0$

n

$n$

A

$A$

-2

Nie jestem pewien, czy to może pomóc, ale myślę, że kiedy mówisz o stabilnym rozwiązaniu, mówisz o metodach aproksymacyjnych. Gdy obliczasz rzeczy wprost, stabilność nie ma sensu. Oznacza to, że musisz zaakceptować przybliżone rozwiązanie, jeśli chcesz zyskać na stabilności.

— ctNGUYEN
źródło

5

Przybliżenie zmiennoprzecinkowe (zaokrąglenie, anulowanie itp.) Liczy się, jeśli chodzi o stabilność. Nawet jeśli masz wzór na odpowiedź, musisz dowiedzieć się, czy można ją dokładnie obliczyć w arytmetyce o skończonej precyzji.

— Bill Barth

Nie uważam tej odpowiedzi za tak negatywną, jak inni ją widzą. Oczywiście problem stabilności istnieje również w przypadku wyraźnych wyników. Uważam jednak, że ctNGUYEN chciał tylko powiedzieć, że przybliżone rozwiązanie, takie jak ekspansja w małej ilości, może faktycznie być bardziej precyzyjne niż pełny, wyraźny wynik, który moim zdaniem jest poprawny. W pewnym sensie proszę o wyraźne rozwiązania, które traktują tak trudne przypadki, aby formuła zawsze była stabilna.

— highsciguy