\documentclass[10pt,a4paper,polish]{article}

%% Polonizacja
\usepackage[T1]{fontenc}
\usepackage[latin2]{inputenc}
\usepackage{polski}
\usepackage{indentfirst}

%% Wzorki matematyczne
\usepackage{amsmath,amsfonts,amssymb}


%% Twierdzenia
\newtheorem{thm}{Twierdzenie}[section]
\newtheorem{df}{Definicja}[section]

\DeclareMathOperator{\sign}{sign}

\numberwithin{equation}{subsection}

\newcommand*{\n}[1][n]{$#1$\nobreakdash-\hspace{0pt}}
\newcommand*{\abs}[1]{\ensuremath{\lvert #1 \rvert}}
\newcommand*{\norm}[1]{\ensuremath{\lVert #1 \rVert}}
\newcommand*{\vect}[1]{\ensuremath{\mathbf{#1}}}
\newcommand*{\matr}[1]{\ensuremath{\mathbf{#1}}}

\def\v{\vspace{-1.4 ex}}



\begin{document}

\newcounter{seria}
\setcounter{seria}{3}
% 1 - sortowanie
% 2 - układy równań
% 3 - układy równań: iteracyjnie
\newcommand*{\preparedby}{Jakub Narębski}
\input{naglowek2.tex}


\section{Metody iteracyjne dla układu równań liniowych}
\label{sec:iter-lin}

\subsection{Iteracyjne poprawianie rozwiązania}
\label{sec:iter-poprawianie}

Rozwiązanie układu równań $\matr{A}\vect{x} = \vect{b}$ którąś z~metod
dokładnych (np. metodą Gaussa) zazwyczaj obarczone jest pewnym
błędem.  Błąd ten możemy oszacować \emph{a posteriori} obliczając
tzw. \textbf{wektor reszt}
\begin{equation}
  \label{eq:wektor-reszt}
  \vect{r} = \matr{A}\,\vect{x} - \vect{b}.
\end{equation}
\emph{Uwaga:} wektor (jego składowe) ten należy obliczać ze zwiększoną
dokładnością, np. w~arytmetyce podwójnej precyzji, gdyż w~przeciwnym
przypdaku możemy błędnie wnioskować o~tym, czy otrzymane rozwiązanie
spełnia układ równań, czy też nie.

Załóżmy, że dokładny wynik $\vect{\hat{x}}$ różni się od otrzymanego
o~małą poprawkę: $\vect{\hat{x}} = \vect{x} - \delta\vect{x}$.  Łatwo
zauważyć, że $\delta\vect{x}$ spełnia równanie 
\begin{equation}
  \label{eq:dx-rown}
  \matr{A}\, \delta\vect{x} = \matr{A}\,\vect{x} - \vect{b} = \vect{r}.
\end{equation}
Możemy zatem obliczyć $\delta\vect{x}$ rozwiązując powyższe równanie.
Oczywiście wynik znów będzie opatrzony pewnym błędem, tak więc
$\vect{x} - \delta\vect{x}$ nie będzie ścisłym rozwiązaniem, a~pewnym
(lepszym niż $\vect{x}$) jego przybliżeniem.

Z~powyższych uwag wynika następujący przepis na \textbf{iteracyjne
  poprawianie rozwiązania} układu równań liniowych:
\begin{enumerate}
\item Obliczamy wektor reszt $\vect{r}^{(i)} = \matr{A}\,\vect{x}^{(i)} - \vect{b}$.
\v
\item Sprawdzamy warunek stopu.
\v
\item Obliczamy poprawkę ${\delta\vect{x}}^{(i)}$ rozwiązując układ równań
  $\matr{A}\, {\delta\vect{x}}^{(i)} = \vect{r}^{(i)}$.
\v
\item Jako następne przybliżenie rozwiązania bierzemy $\vect{x}^{(i+1)}
  = \vect{x}^{(i)} - {\delta\vect{x}}^{(i)}$.
\end{enumerate}

Jako warunek stopu możemy przyjąć $\norm{\vect{r}^{(i)}}_2 < \varepsilon$,
lub
\begin{math}
  \norm{\vect{r}^{(i)}}_{\infty} < \varepsilon \norm{\matr{A}\,\vect{x}^{(i)}}_{\infty},
\end{math}
albo
\begin{math}
  \norm{\vect{r}^{(i)}}_{\infty} < \varepsilon \norm{\vect{b}}_{\infty},
\end{math}
(o ile $\vect{x}^{(i)}$ jest dostatecznie bliski rozwiązania).  Norma
wektorowa $\norm{\bullet\null\,}_{\infty}$ zdefiniowana jest wzorem
$\norm{\vect{b}}_{\infty} = \max(\abs{b_1},
\abs{b_2},\dots,\abs{b_n})$.



\subsection{Metody iteracyjne rozwiązujące $\matr{A}\,\vect{x} =
  \vect{b}$ oparte na podziale macierzy}
\label{sec:iter-rozw}

W~wielu przypadkach pojawia się potrzeba rozwiązywania wielkich
układów równań z~macierzą \emph{rzadką} tzn. o~liczbie elementów
niezerowych $\mathcal{O}(n)$.  Klasyczne metody bezpośrednie
(np. metoda Gaussa) zazwyczaj niszczą własność rzadkości, co powoduje
zwiększone zapotrzebowanie na pamięć.  Działają w~czasie
$\mathcal{O}(n^3)$ co jest zbyt kosztowne.  Jeśli macierz rzadka ma
szczególną postać to można użyć odpowiedniego algorytmu.

Innym sposobem rozwiązywania układu równań liniowych jest zastosowanie
metody iteracyjnej.  Wymagają one w~każdym kroku ilości operacji
proporcjonalnej do kosztu obliczenia $\matr{A}\,\vect{x}$, dla
macierzy rzadkich równemu $\mathcal{O}(n)$.  Metody te nadają się także
do wyznaczenia rozwiązania z~zadaną niewielką dokładnością niewielkim
kosztem.

Najprostszym rodzajem metod iteracyjnych są metody oparte na podziale
macierzy
\begin{equation}
  \label{eq:podzial-macierzy}
  \matr{A} = \matr{M} - \matr{N}.
\end{equation}
Rozwiązanie układu $\matr{A}\,\vect{x} = \vect{b}$ jest punktem stałym
$\matr{M}\,\vect{x} = \matr{N}\,\vect{x} + \vect{b}$.  Punkt stały
możemy znaleźć za pomocą \emph{iteracji prostej}
\begin{equation}
  \label{eq:iteracja-prosta}
  \matr{M}\,\vect{x}^{(i+1)} = \matr{N}\,\vect{x}^{(i)} + \vect{b}.
\end{equation}
Podział macierzy~\eqref{eq:podzial-macierzy} trzeba tak wybrać, by
macierz $\matr{M}$ była nieosobliwa i~by układ równań z~tą macierzą
było łatwo rozwiązywać.  Ponadto aby metoda iteracji prostej była
zbieżna każda wartość własna macierzy $\matr{M}^{-1}\matr{N}$ powinna
być co do modułu mniejsza od jedności (promień spektralny tej macierzy
powinien być mniejszy od jedności $\varrho(\matr{M}^{-1}\matr{N})
< 1$).



\subsubsection{Metoda Richardsona}
\label{sec:metoda-richardsona}

Najprostszą metodą iteracyjną jest \textbf{metoda Richardsona}
określona wzorem
\begin{equation}
  \label{eq:metoda-Richardsona}
  \vect{x}^{(i+1)} = \vect{x}^{(i)} - \tau (\matr{A}\vect{x}^{(i)} - \vect{b}).
\end{equation}
Aby ta metoda była zbieżna wszystkie wartości własne macierzy
$\matr{A}$ muszą być tego samego znaku.  Zachodzi to np. dla macierzy
symetrycznych dodatnio określonych.

Niech $\matr{A} = \matr{A}^T > 0$.  Oznaczmy przez
$\lambda_{\text{min}}$ najmniejszą, a~przez $\lambda_{\text{max}}$
największą wartość własną macierzy $\matr{A}$: 
\(0 < \lambda_{\text{min}} \leq \lambda_{\text{max}}\).  Metoda Richardsona
jest zbieżna gdy \(0 < \tau < 2/\lambda_{\text{max}}\).  Najszybszą
zbieżność osiągniemy dla 
\begin{equation}
  \label{eq:tau-opt}
  \tau_{\text{opt}} = \frac{2}{\lambda_{\text{min}} + \lambda_{\text{max}}}.
\end{equation}



\subsubsection{Metoda Jacobiego}
\label{sec:metoda-jacobiego}

Przedstawmy macierz $\matr{A}$ wyjściowego układu $\matr{A}\,\vect{x}
= \vect{b}$ w~postaci
\begin{equation}
  \label{eq:Jacobi-rozklad}
  \matr{A} = \matr{L} + \matr{D} + \matr{U},
\end{equation}
gdzie $\matr{L}$ jest macierzą poddiagonalną, $\matr{D}$ jest macierzą
diagonalną zaś $\matr{U}$ macierzą naddiagonalną.  W~metodzie
Jacobiego zakłada się, że macierz $\matr{D}$ jest odwracalna tzn. ma
wszyskie elementy na diagonali niezerowe.  Można to uzyskać
zamieniając miejscami wiersze macierzy $\matr{A}$ (podobnie jak
w~metodzie Gaussa z~częsciowym wyborem elementu podstawowego).  Jeśli
takich zamian nie trzeba dokonywać to $\matr{D} =
\mathrm{diag}(a_{ii})$, $\matr{L} = (a_{ij})_{i>j}$, zaś $\matr{U} =
(a_{ij})_{i<j}$.

Metoda ta zadana jest \textbf{wzorem Jacobiego}:
\begin{equation}
  \label{eq:metoda-Jacobiego}
  \matr{D}\,\vect{x}^{(i+1)} = 
  -(\matr{L} + \matr{U})\,\vect{x}^{(i)} + \vect{b}.
\end{equation}
Macierz odwrotna do macierzy $\matr{D}$ dana jest wzorem
\begin{equation}
  \label{eq:D-odwrotna}
  \matr{D}^{-1} = \mathrm{diag}(a^{-1}_{ii}).
\end{equation}

Warunkiem dostatecznym zbieżności tej metody jest by macierz
$\matr{A}$ (ewentualnie z~przestawionymi wierszami) była \emph{silnie
  diagonalnie dominująca} lub \emph{symetryczna dodatnio określona}.



\subsubsection{Metoda Gaussa--Seidla}
\label{sec:metoda-gaussa-seidla}

Przedstawmy macierz $\matr{A}$ wyjściowego układu $\matr{A}\,\vect{x}
= \vect{b}$ w~postaci
\begin{equation}
  \label{eq:GS-rozklad}
  \matr{A} = \matr{L} + \matr{D} + \matr{U}.
\end{equation}

\textbf{Metoda Gaussa--Seidla} zadana jest wzorem:
\begin{equation}
  \label{eq:metoda-GS}
  (\matr{L} + \matr{D})\,\vect{x}^{(i+1)} = 
  -\matr{U}\,\vect{x}^{(i)} + \vect{b}.
\end{equation}
Należy zauważyć, że macierz $(\matr{L} + \matr{U})$ jest macierzą
trójkatną górną.  Układ równań z~taka macierzą możemy rozwiązać
w~czasie $\mathcal{O}(n^2)$ czyli takim jak mnożenie macierzy przez
wektor.

Warunkiem koniecznym zbieżności tej metody jest by macierz $\matr{A}$
(ewentualnie z~przestawionymi wierszami) była \emph{diagonalnie
  dominująca} lub \emph{symetryczna dodatnio określona}.  Metoda
Gaussa--Seidla jest zazwyczaj szybciej zbieżna niż metoda Jacobiego.



\subsubsection{Metoda SOR}
\label{sec:metoda-SOR}

Metodę Gaussa--Seidla można przyspieszyć mnożąc poprawkę przez
odpowiednio dobraną liczbę $\omega$.  Można to sprowadzić do wzoru:
\begin{equation}
  \label{eq:metoda-SOR}
  (\omega\matr{L} + \matr{D})\,\vect{x}^{(i+1)} =
  (1 - \omega) \matr{D}\,\vect{x}^{(i)} - \omega\matr{U}\,\vect{x}^{(i)} + \omega\vect{b}.
\end{equation}
Dla $\omega > 1$ metoda ta nosi nazwę \textbf{metody kolejnych
  nadrelaksacji} (w~skrócie SOR, od ang. succesive overrelaxation).
Optymalną wartością parametru $\omega$ dla macierzy $\matr{A}$
symetrycznej, dodatnio określonej i~mającej własność $A$ jest
\begin{equation}
  \label{eq:omega-opt}
  \omega_{\text{opt}} = \frac{2}{1 + \sqrt{1 - \lambda^2_{\text{max}}}},
\end{equation}
gdzie $\lambda_{\text{max}}$ jest największą co do modułu wartością
własną $\matr{D}^{-1}(\matr{L} + \matr{U})$.  Aby metoda była zbieżna

dla kazdych wartości parametrów początkowych musi zachodzić $0 <
\omega < 2$.



\section{Uzupełnienie zadania domowego}
\label{sec:uzup-zad-dom}

Napisać funkcję dokonującą jednego kroku iteracyjnego poprawienia
rozwiązania równania $\matr{A}\,\vect{x} = \vect{b}$, funkcję
poprawiającą uzyskane rozwiązanie tego równania do osiągnięcia zadanej
dokładności (lub przekroczenia zadanej maksymalnej ilości iteracji).
Funkcja jako argumenty może przyjmować macierz $\matr{A}$, jej rozkład
$\matr{L}\matr{U}$, wektor $\vect{b}$ oraz wektor $\vect{x}$
zawierający początkowe przybliżenie.  Funkcje mogą nadpisywać wektor
$\vect{x}$.  Jeśli funkcja potrzebuje dodatkowego wektora pomocniczego
można go przekazać jako parametr.  Jako warunek stopu można przyjąć
\begin{equation}
  \label{eq:warunek-stopu}
  \norm{\vect{r}} < \varepsilon_{\text{abs}} + 
  \varepsilon_{\text{rel}} \norm{\matr{A}\,\vect{x}}.
\end{equation}

Napisane funkcje przetestować za pomocą programu korzystającego
z~rozkładu $\matr{L}\matr{U}$ macierzy do rozwiązywania układu równań
(zadanie 2) oraz powyższych funkcji do poprawiania rozwiązania.
\emph{Uwaga:} jeśli funkcja dokonująca rozkładu $\matr{L}\matr{U}$
nadpisuje macierz $\matr{A}$ należy najpierw zrobić kopię macierzy do
wykorzystania w~metodach poprawiających rozwiązanie.

Za uzupełnienie zadania domowego można otrzymać maksymalnie $1/3$ pkt.

%%% 
\input{podpis.tex}

\end{document}
