Lab. 12 Układy równań liniowych - metody iteracyjne

Układ równań

Poszukujemy wektora $\mathbf{x}$ który spełnia

$$\mathbf{A} \mathbf{x}=\mathbf{b}$$

$$\left[\begin{array}{cccc} a_{11} & a_{12} & \ldots & a_{1 n} \\ a_{21} & a_{22} & \ldots & a_{2 n} \\ \ldots & \ldots & \ldots & \ldots \\ a_{n 1} & a_{n 2} & \ldots & a_{n n} \end{array}\right] \cdot\left[\begin{array}{l} x_{1} \\ x_{2} \\ \ldots \\ x_{n} \end{array}\right]=\left[\begin{array}{l} b_{1} \\ b_{2} \\ \ldots \\ b_{n} \end{array}\right]$$

Metody

• Dokładne:
  • wyznaczniki Cramera
  • Eliminacja Gaussa
  • Metoda Gaussa-Jordana
  • Metoda Doolittle'a (rozkład LU)
  • Metoda Choleskiego
• Iteracyjne (przybliżone):
  • Metoda Jacobiego
  • Metoda Gaussa-Seidla

Metoda Jacobiego

$$\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}+\mathbf{D}^{-1}(\mathbf{b}- \mathbf{A} \mathbf{x}^{(k)})$$

• macierz $\mathbf{D}$ to macierz diagonalna zawierająca na diagonali elementy $a_{ii}$
• metoda jest dobrze zdefiniowana dla $a_{ii} \neq 0$
• jest zbieżna dla macierzy o dominujących wartościach diagonalnych $\left|a_{i, i}\right|>\sum_{j \neq i}\left|a_{i, j}\right| \quad \text{ lub } \quad\left|a_{i, i}\right|>\sum_{j \neq i}\left|a_{j, i}\right|, \quad i=1,2, \ldots, n$
• złożoność obliczeniowa $ O(n^2)$, stosowana do dużych układów równań

Algorytm

1. ustalamy przybliżenie początkowe $\mathbf{x}^{(0)}$
2. kolejne przybliżenia $\mathbf{x}^{(1)}, \mathbf{x}^{(2)}, \ldots $ uzyskujemy z wzoru
   $$ x_i^{(k+1)} = \left( b_i - \sum_{j=1\\j\neq i}^n a_{ij} x_{j}^{(k)}\right) \bigg/ a_{ii} $$
3. koniec obliczeń gdy $|| \mathbf{x}^{(k+1)} - \mathbf{x}^{(k)}|| < \epsilon$

Przykład w C

k = 0;
do
{
  for (i = 0; i < n; i++) x_prev[i] = x[i];
 
  for (i = 0; i < n; i++)
  {
      x[i] = b[i];
      for (j = 0; j < n; j++)
      {
        if (i!=j) x[i] -= x_prev[j] * A[i][j];
      }
      x[i] /= A[i][i];
  }
  k = k + 1;
}while(dist(x, x_prev, n) > eps && k < max_iter);

Program demo: jacobi.c

Metoda Gaussa-Seidla

Zauważmy, że w metodzie Jacobiego wyznaczając komponent $x_i^{(k+1)}$ możemy wykorzystać wyznaczone w tym samym kroku $(k+1)$ wcześniejsze komponenty $ x_1^{(k+1)}, x_2^{(k+1)}, \ldots, x_{i-1}^{(k+1)}$

Zalety: szybsza zbieżność, do wykonania obliczeń nie trzeba przechowywać dwóch wektorów $\mathbf{x}$

Algorytm

1. ustalamy przybliżenie początkowe $\mathbf{x}^{(0)}$
2. kolejne przybliżenia $\mathbf{x}^{(1)}, \mathbf{x}^{(2)}, \ldots $ uzyskujemy z wzoru
   $$ x_i^{(k+1)} = \left( b_i - \sum_{j=1}^{i-1} a_{ij} x_{j}^{(k+1)} - \sum_{j=i+1}^{n} a_{ij} x_{j}^{(k)} \right) \bigg/ a_{ii} $$
3. koniec obliczeń gdy $|| \mathbf{x}^{(k+1)} - \mathbf{x}^{(k)}|| < \epsilon$

Zadanie

Napisz program, który porówna szybkość zbieżności metod iteracyjnych Jacobiego i Gaussa-Siedla przy rozwiązywaniu układu równań

$$\left[\begin{array}{cccc} n+1 & 1 & \cdots & 1 \\ 1 & n+1 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & n+1 \end{array}\right]\left[\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array}\right]=\left[\begin{array}{c} 2 n \\ 2 n \\ \vdots \\ 2 n \end{array}\right] $$

Dokładne rozwiązanie $ \mathbf{x} = [1, 1, \ldots, 1 ]$.

Program, dla podanej przez użytkownika wartości $n$, rozwiązuje układ równań za pomocą obu metod iteracyjnych w każdym przypadku starając od przybliżenia $\mathbf{x}^{(0)} = [0, 0, \ldots, 0]$ i wykonuje obliczenia aż do uzyskania wyniku o precyzji $\epsilon = 10^{-4}$. Po skończeniu obliczeń program wypisuje dla każdej z metod informację o liczbie wykonanych iteracji.
Zakładamy, że $n$ nie będzie większe od 1000.

Dodatkowe ćwiczenia

Ćw. 1