Układy liniowe z macierzami rzadkimi
<<< Powrót do strony głównej
przedmiotu Metody numeryczne
Oglądaj wskazówki i rozwiązania __SHOWALL__
Ukryj wskazówki i rozwiązania __HIDEALL__
Ćwiczenie: Metoda Richardsona
Rozwiązanie
Niech ekstremalnymi wartościami własnymi macierzy
będą
. Ponieważ macierz iteracji metody Richardsona
, to jej wartości własne muszą leżeć w przedziale
i --- oczywiście --- aby iteracja miała sens,
(dlaczego?).
Co więcej,
wtedy i tylko wtedy, gdy
, a najmniejszą normę spektralną macierzy
uzyskamy, gdy
i wówczas
(Przy okazji zauważ, że gdyby macierz
nie była określona, tzn. miałaby zarówno dodatnie, jak i ujemne wartości własne, to metoda Richardsona mogłaby w ogóle nie być zbieżna (dla pewnych wektorów startowych))
Jeśli chodzi o mnożenie przez macierz jednowymiarowego laplasjanu, to najprościej wcale nie używać struktury macierzy! Rzeczywiście, jedyne, czego potrzebujemy to operacja mnożenia wektora przez macierz
, a to realizuje pętla:
void LapMult(double *x, int N, double *y)
/*
Mnożenie wektora przez macierz 1-wym laplasjanu L wymiaru N
Wejście:
x - wektor, który mnożymy przez L
N - jego długość
Wyjście:
y - wynik: y = Lx
*/
{
int i;
for (i = 1; i < N-1; i++)
y[i] = -x[i-1] + 2.0*x[i] - x[i+1];
y[0] = 2.0*x[0] - x[1];
y[N-1] = -x[N-2] + 2.0*x[N-1];
}
Ćwiczenie
Zaimplementuj operacje:
- mnożenia macierzy
przez wektor
,
- wyłuskania wartości elementu
,
- zmiany wartości pewnego zerowego wyrazu macierzy na niezerową,
jeśli macierz jest zadana w formacie
Przetestuj dla kilku macierzy z kolekcji
MatrixMarket.
Rozwiązanie
Dodawanie nowego elementu w formacie AIJ jest bardzo łatwe. W przeciwieństwie do
pozostałych, w których trzeba zachowywać uporządkowanie. Mnożenie przez wektor
najwygodniejsze jest w CSR, bo dodatkowo narzuca zasadę lokalności w przestrzeni. Wyłuskanie wartości jest najmniej efektywne w
formacie AIJ. Szczegóły opisane są w rozdziale 3.5 książki
Zobacz także implementacje w Fortranie, w pakiecie SPARSKIT, będącym czymś w rodzaju odpowiednika BLAS dla macierzy rozrzedzonych.
Ćwiczenie: Konwersja formatu macierzy rzadkiej
Napisz procedurę aij2csr
, konwertującą macierz w formacie AIJ do CSR i csr2aij
, działającą w drugą stronę.
Rozwiązanie
Zobacz, jak to zrobiono w pakiecie SPARSKIT.
Ćwiczenie
Jak tanio rozwiązywać układ równań z macierzą cykliczną trójdiagonalną, tzn.
Dla uproszczenia załóżmy, że macierz jest dodatnio określona i symetryczna.
Zaimplementuj opracowaną metodę, korzystając z BLASów i LAPACKa.
Wskazówka
Dziel i rządź! (Tylko: jak?)
Ćwiczenie: CGNE
Ktoś mógłby sugerować, że skoro CG działa tylko dla macierzy symetrycznych, to dowolny układ
z macierzą nieosobliwą można transformować do równoważnego mu układu równań normalnych
którego macierz
jest już oczywiście macierzą symetryczną i dodatnio określoną.
Wskaż potencjalne wady tej metody i podaj sposób jej implementacji.
Wskazówka
Jakie jest uwarunkowanie macierzy

?
Rozwiązanie
Nietrudno sprawdzić, że dla normy spektralnej macierzy,
a więc w przypadku macierzy źle uwarunkowanych należy spodziewać się patologicznie dużej liczby iteracji. Chociaż dobre (symetryczne, dodatnio określone) imadło
mogłoby pomóc, np.
to jednak znacznie lepiej stosować metody opracowane specjalnie dla macierzy niesymetrycznych, np. GMRES (oczywiście z nieodzownym ściskaniem macierzy, gdy jest źle uwarunkowana...).
Implementacja metody iteracyjnej to tylko decyzja, jak realizować mnożenie przez
. Niedobra metoda to
B = <math>\displaystyle A^TA</math>;
...
while ...
y = B*x;
end
gdyż
będzie bardziej wypełniona niż A. Znacznie lepiej
...
while ...
y = A*x;
y = (y'*A)';
end
co realizuje się kosztem równym dwukrotnemu mnożeniu przez macierz
(w formacie AIJ) i nie wymaga dodatkowej pamięci.