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

Mam niejednorodny układ liniowy

gdzie jest rzeczywistą macierzą n × n przy n ≤ 4 . Nullspace z A jest gwarancją wymiaru zerowym więc równanie ma unikalną odwrotną x = A - 1 b . 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 A i b . Z powodu tego wymogu i małych wymiarowości Myślałem wdrożyć wyraźne wzory na A - 1$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.

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.

$\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.

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.