Układy równań liniowych

(1)

Układy równań liniowych

Marcin Orchel

1 Zadania

1.1 Zadania na 3.0

Napisać skrypt w R. W skrypcie:

• Porównać jakość rozwiązania dwóch wybranych układów równań liniowych dwoma wybranymi metodami, w tym co najmniej jedną metodą iteracyjną.

• przetestować metodę rozkładu LU

• przetestować metodę dekompozycji QR lub metodę Cholesky decomposition

• sprawdzić czy testowane układy równań mają jedno rozwiązanie, brak rozwiązań lub nieskończenie wiele rozwiązań za pomocą sprawdzenia rzędu macierzy A, rzędu macierzy uzupełnionej i liczby zmiennych

• wykonać testy dla macierzy osobliwej

• sprawdzić uwarunkowanie macierzy, wykonać testy dla macierzy źle uwarunkowa- nych

• wykonać testy dla różnych wielkości macierzy A

• za pomocą zapisu macierzowego rozwiązać równania z tą samą macierzą A i róż- nymi wektorami b

• dodać błąd normalny do współczynników i wykonać testy jakości dla różnych od- chyleń standardowych dodawanego błędu

• pozostawić pewną ilość współczynników zerową, przykładowo 30% wygenerowa- nych współczynników jest zerami, wykonać testy dla różnej ilości wyzerowanych współczynników

• przetestować szybkość wybranych metod, wyświetlić na wykresie zależność czasu od rozmiaru macierzy A

• Dodać komentarz do skryptu opisujący krótko na czym polegają użyte metody oraz wnioski z badań.

(2)

Wskazówki

• Jakość możemy sprawdzić tak, że po znalezieniu wektora x obliczyć Ax i porównać wyliczone b z oryginalnym b wyliczając błąd RMS. Oryginalne b to b przed doda- niem błędu normalnego w przypadku macierzy z dodanym błędem normalnym.

• jedną z wybranych macierz może być macierz losowa. Współczynniki A i b mogą być losowane według rozkładu jednostajnego.

• alternatywnie można wybrać macierz A i wektor x, następnie wyliczyć b, i wyli- czamy wtedy błąd na wektorze x

• Ilość rozwiązań układu Ax = b. Brak rozwiązań wtw gdy rank(A) < rank[A|b].

Jedno rozwiązanie wtw gdy rank(A) = rank(A|b) = n, gdzie n to liczba zmien- nych. Nieskończenie wiele rozwiązań, wtw gdy rank[A] = rank[A|b] < n.

• Warunek rzadkości w compressed sensing. Dla równania 2x + y = 1, mamy dwa rozwiązania rzadkie z jednym zerem x = 0 i y = 1 oraz x = 1/2 i y = 0.

Wskazówki do R

• https://stat.ethz.ch/R-manual/R-devel/library/base/html/solve.html, po- danie b,http://www.netlib.org/lapack/explore-html/d7/d3b/group__double_

g_esolve.html#ga5ee879032a8365897c3ba91e3dc8d512

• https://stat.ethz.ch/R-manual/R-devel/library/base/html/qr.html

• backsolve, forwardsolve,https://stat.ethz.ch/R-manual/R-devel/library/base/

html/backsolve.html, możliwość podania wielu prawych stron równania

• https://www.rdocumentation.org/packages/matrixcalc/topics/lu.decomposition

• pakiet limSolve

• https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/ginv.html,https:

//www.rdocumentation.org/packages/ibdreg/topics/Ginv

• https://stat.ethz.ch/R-manual/R-devel/library/Matrix/html/sparseQR-class.

html

• https://www.rdocumentation.org/packages/systemfit/topics/systemfit

• https://stat.ethz.ch/R-manual/R-devel/library/stats/html/toeplitz.html, https://www.rdocumentation.org/packages/pracma/topics/toeplitz,https:

//www.rdocumentation.org/packages/matrixcalc/topics/toeplitz.matrix

• https://stat.ethz.ch/R-manual/R-devel/library/Matrix/html/Hilbert.html, https://www.rdocumentation.org/packages/matrixcalc/topics/hilbert.matrix

(3)

• https://www.rdocumentation.org/packages/matrixcalc/topics/is.singular.

matrix

• https://stat.ethz.ch/R-manual/R-devel/library/Matrix/html/rankMatrix.

html

• https://www.rdocumentation.org/packages/Matrix/topics/KNex Wskazówki do Matlaba

• do rozwiązywania układów równań można użyć operatora \, możemy sprawdzić w dokumentacji jakie metody rozwiązywania układów równań są zaimplementowane dla tego operatora,http://www.mathworks.com/help/matlab/ref/mldivide.html, lub metody linsolve http://www.mathworks.com/help/matlab/ref/linsolve.

html lub odwracania macierzy http://www.mathworks.com/help/matlab/ref/

inv.html, a także różnych metod iteracyjnych, przykładowohttp://www.mathworks.

com/help/matlab/ref/lsqr.html,http://www.mathworks.com/help/matlab/ref/

minres.html,http://www.mathworks.com/help/matlab/ref/gmres.html,http:

//www.mathworks.com/help/matlab/ref/tfqmr.html,http://www.mathworks.

com/help/matlab/ref/bicgstab.html,http://www.mathworks.com/help/matlab/

ref/pcg.html,http://www.mathworks.com/help/matlab/ref/bicg.html

• sprawdzenie rzędu macierzy wykonuje się za pomocą funkcjihttp://www.mathworks.

com/help/matlab/ref/rank.html

• sprawdzenie uwarunkowania macierzy wykonuje się za pomocą funkcji cond,http:

//www.mathworks.com/help/symbolic/cond.html

• macierz jednostkowa polecenie eyehttp://www.mathworks.com/help/matlab/ref/

eye.html

• macierz ze wszystkimi jedynkami oneshttp://www.mathworks.com/help/matlab/

ref/ones.html, zerami zeroshttp://www.mathworks.com/help/matlab/ref/zeros.

html

• do rozkładu LU można użyć poleceniahttp://www.mathworks.com/help/matlab/

ref/lu.html

• przykład macierzy rzadkiej z matlaba west0479, np. używany tutaj http://www.

mathworks.com/help/matlab/math/accessing-sparse-matrices.html

• można dokonać ręcznej permutacji wierszy lub kolumn np. poleceniem colamd http://www.mathworks.com/help/matlab/ref/colamd.htmllub symamdhttp:

//www.mathworks.com/help/matlab/ref/symamd.html

• do tworzenia macierzy diagonalnych można użyć spdiagshttp://www.mathworks.

com/help/matlab/ref/spdiags.html

(4)

• generacja macierzy rzadkich losowych np.http://www.mathworks.com/help/matlab/

ref/sprand.html,http://www.mathworks.com/help/matlab/ref/speye.html

• sprawdzenie jak bardzo macierz jest rzadka nnz http://www.mathworks.com/

help/matlab/ref/nnz.html, sprawdzenie czy macierz jest symetryczna issymme- trichttp://www.mathworks.com/help/matlab/ref/issymmetric.html, czy jest diagonalna isdiaghttp://www.mathworks.com/help/matlab/ref/isdiag.html

• konwersja macierzy rzadkiej do pełnej full http://www.mathworks.com/help/

matlab/ref/full.html

• zerowanie elementów macierzyhttp://www.mathworks.com/matlabcentral/fileexchange/

2360-the-matrix-computation-toolbox/content/matrixcomp/sparsify.m

• predefiniowane macierze w Matlabie,http://www.mathworks.com/help/matlab/

ref/gallery.html, magichttp://www.mathworks.com/help/matlab/ref/magic.

html, hadamardhttp://www.mathworks.com/help/matlab/ref/hadamard.html, macierz źle uwarunkowana hilb http://www.mathworks.com/help/matlab/ref/

hilb.html

Przykłady w R

• A<- matrix(c(1, 3, 3, -1, 8, -1, 4, -9, -7), nrow=3, byrow=TRUE) B<- matrix(c(1, 2, 3, 4, 5, 6), nrow=3, byrow=TRUE)

solve(A, B) A %*% solve(A,B)

1.2 Zadania na 4.0

• Zaimplementować metodę liniową najmniejszych kwadratów ze wzoru

x =A^TA⁻¹A^Tb (1)

• Dla jakiego przypadku tez wzór jest prawidłowy.

• Porównać rezultat z wynikiem działania operatora \

• Przetestować metodę dla przypadku gdy macierz A ma większą liczbę wierszy niż kolumn.

Wskazówki

• http://en.wikipedia.org/wiki/Overdetermined_system

(5)

1.3 Zadania na 5.0

• Rozwiązać zadanie na 4.0 z dekompozycją QR, to znaczy po podstawieniu QR do (1) za macierz A otrzymujemy

Rx = Q^Tb , (2)

oraz

x = R⁻¹Q^Tb , (3)

gdzie macierz R jest macierzą górnotrójkątną oraz porównać wyniki z wynikami z zadania na 4.0.

• Rozwiązać metodą postulowanego rozwiązania wybrane równanie różniczkowe.