~~NOCACHE~~
~~REVEAL theme=simple&size=1024x800~~

====== Lab. 11 Rozwiązywanie układów równań liniowych c. d.======

===== 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 Gausa-Seidla

===== Rozkład LU =====

Dekompozycja macierzy kwadratowej $\mathbf{A}$ na iloczyn macierzy trójkątnej dolnej $\mathbf{L}$ oraz macierzy trójkątnej górnej $\mathbf{U}$

$$ \mathbf{A}= \mathbf{L} \mathbf{U}$$

gdzie
$$\mathbf{L}=\left[\begin{array}{cccc}
1 & 0 & . . & 0 \\
l_{21} & 1 & . . & 0 \\
. . & . . & . . & . . \\
l_{n 1} & l_{n 2} & . . & 1
\end{array}\right] \quad \mathbf{U}=\left[\begin{array}{cccc}
u_{11} & u_{12} & . . & u_{1 n} \\
0 & u_{22} & . . & u_{2 n} \\
. . & . . & . . & . . \\
0 & . . & 0 & u_{n n}
\end{array}\right]$$

  * metoda Doolitle'a (gdy $l_{ii}=1$)
  * metoda Choleskiego-Banachiewicza (gdy $l_{ii}=u_{ii}$) 

Tadeusz Banachiewicz 1938

===== Rozwiązanie układu równań =====

$$ \mathbf{A}\mathbf{x} = \mathbf{b}$$
po dekompozycji
$$ \mathbf{LU}\mathbf{x} = \mathbf{b}$$
sprowadza się do rozwiązania dwóch prostszych problemów
$$ \mathbf{L}\mathbf{y} = \mathbf{b}$$
$$ \mathbf{U}\mathbf{x} = \mathbf{y}$$

===== Metoda Doolitle'a =====

$$\mathbf{A}=\left[\begin{array}{llll}
a_{11} & a_{12} & a_{13} & a_{14} \\
a_{21} & a_{22} & a_{23} & a_{24} \\
a_{31} & a_{32} & a_{33} & a_{34} \\
a_{41} & a_{42} & a_{43} & a_{44}
\end{array}\right]=\left[\begin{array}{cccc}
1 & 0 & 0 & 0 \\
l_{21} & 1 & 0 & 0 \\
l_{31} & l_{32} & 1 & 0 \\
l_{41} & l_{42} & l_{43} & 1
\end{array}\right]\left[\begin{array}{rrrr}
u_{11} & u_{12} & u_{13} & u_{14} \\
0 & u_{22} & u_{23} & u_{24} \\
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{array}\right]$$

$$=\left[\begin{array}{llll}
u_{11} & u_{12} & u_{13} & u_{14} \\
l_{21} u_{11} & l_{21} u_{12}+u_{22} & l_{21} u_{13}+u_{23} & l_{21} u_{14}+u_{24} \\
l_{31} u_{11} & l_{31} u_{12}+l_{32} u_{22} & l_{31} u_{13}+l_{32} u_{23}+u_{33} & l_{31} u_{14}+l_{32} u_{24}+u_{34} \\
l_{41} u_{11} & l_{41} u_{12}+l_{42} u_{22} & l_{41} u_{13}+l_{42} u_{23}+l_{43} u_{33} & l_{41} u_{14}+l_{42} u_{24}+l_{43} u_{34}+u_{44}
\end{array}\right]$$

===== =====

{{  zajecia:mn1_2021_2:lu_step1.png?800  |}}

wiersz 1  $$u_{11}=a_{11} \quad u_{12}=a_{12}  \quad u_{13}=a_{13} \quad u_{14}=a_{14}$$ 
kolumna 1 $$l_{21}=\frac{a_{21}}{u_{11}} \quad l_{31}=\frac{a_{31}}{u_{11}} \quad l_{41}=\frac{a_{41}}{u_{11}}$$

===== =====

{{  zajecia:mn1_2021_2:lu_step2.png?800  |}}

wiersz 2
$$u_{22}=a_{22}-l_{21} u_{12} \quad u_{23}=a_{23}-l_{21} u_{13} \quad u_{24}=a_{24}-l_{21} u_{14}$$

kolumna 2 
$$l_{32}=\left(a_{32}-l_{31} u_{12}\right) / u_{22} \quad l_{42}=\left(a_{42}-l_{41} u_{12}\right) / u_{22}$$

===== =====

{{  zajecia:mn1_2021_2:lu_step3.png?800  |}}


wiersz 3
$$u_{33}=a_{33}-l_{31} u_{13}-l_{32} u_{23} \quad u_{34}=a_{34}-l_{31} u_{14}-l_{32} u_{24}$$

kolumna 3
$$l_{43}=\left(a_{43}-l_{41} u_{13}-l_{42} u_{23}\right) / u_{33}$$

wiersz 4
$$u_{44}=a_{44}-l_{41} u_{14}-l_{42} u_{24}-l_{43} u_{34}$$

===== Algorytm rozkładu LU =====

  * Dla $i = 1, \ldots, n$
    * $ u_{ij}=0, \, l_{ji}=0, \qquad j = 1, \ldots, i-1$
    * $l_{ii} = 1 $
    * $u_{ij}=a_{ij}-\sum_{k=1}^{i-1} l_{ik} u_{kj}, \qquad j = i, \ldots, n$
    * $l_{ji}=\frac{a_{j i}-\sum_{k=1}^{i-1} l_{j k} u_{k i}}{u_{i i}}, \qquad j = i+1, \ldots, n$

Każdy element $a_{ij}$ wykorzystywany jest tylko raz, więc można wykorzystać miejsce w macierzy $\mathbf{A}$ do przechowania wartości $l_{ij}$ i $u_{ij}$

===== Macierze: Przykład w C =====

<code C>
#define MAX 100
typedef double matrix[MAX][MAX];

void print_matrix(matrix A, int w, int k)
{
   printf("\nMacierz %dx%x\n", w, k);

   for(int i=0; i< w;i++)
   {
      for(int j=0; j< k; j++)
         printf(" %5.1lf", A[i][j]);
      printf("\n");
   }
}
</code>

===== Mnożenie macierzy: Przykład w C =====

$$ \mathrm{C} = \mathrm{A} \mathrm{B}$$

<code C>
void mul(const matrix A, const matrix B, matrix C, int nwa, int nka, int awb, int nkb)
{
   for (int i = 0; i < nwa; i++)
   {
      for (int j  = 0; j < nkb; j++)
      {
         C[i][j] = 0.0;
         for (int k = 0; k < nka; k++)
         {
            C[i][j] += A[i][k] * B[k][j];
         }
      }
   }
}
</code>


===== Rozkład LU: Przykład w C =====

<code C>
void rozklad_lu(const matrix A, matrix L, matrix U, int n)
{
   for (int i = 0; i < n; i++)
   {
      for (int j = 0; j < i; j++)
      {
         U[i][j] = 0;
         L[j][i] = 0;
      }

      L[i][i] = 1.0;

      for (int j = i; j < n; j++)
      {
         U[i][j] = A[i][j]; 
         for (int k = 0; k <= i-1; k++)
         {
            U[i][j] -= L[i][k] * U[k][j]; 
         }         
      }

      for (int j = i+1; j < n; j++)
      {
         L[j][i] = A[j][i]; 
         for (int k = 0; k <= i-1; k++)
         {
            L[j][i] -= L[j][k] * U[k][i]; 
         }         
         L[j][i] /= U[i][i];
      }
   }
}
</code>

Przykładowy kod {{zajecia:mn1_2021_2:lu.c|}}

===== Rozwiązanie układu równań =====

$$\mathbf{A}\mathbf{x} = \mathbf{b}$$ 
zamieniamy na
$$\mathbf{L}\mathbf{U}\mathbf{x} = \mathbf{b}$$ 
stąd
$$ \mathbf{L}\mathbf{y} = \mathbf{b}$$
$$ \mathbf{U}\mathbf{x} = \mathbf{y}$$

===== =====

Z pierwszego układu $ \mathbf{L}\mathbf{y} = \mathbf{b}$ 

$$ \left[\begin{array}{cccc}
1 & 0 & . . & 0 \\
l_{21} & 1 & . . & 0 \\
. . & . . & . . & . . \\
l_{n 1} & l_{n 2} & . . & 1
\end{array}\right] 
\cdot\left[\begin{array}{l}
y_{1} \\
y_{2} \\
\ldots \\
y_{n}
\end{array}\right]=\left[\begin{array}{l}
b_{1} \\
b_{2} \\
\ldots \\
b_{n}
\end{array}\right]$$

otrzymujemy wektor $\mathbf{y}$

$$\left\{\begin{array}{l}
y_{1}=b_{1} \\
y_{i}=b_{i}-\sum_{k=1}^{i-1} l_{i k} y_{k} \quad i=2,3, \ldots, n
\end{array}\right.$$


===== =====


Z drugiego układu $ \mathbf{U}\mathbf{x} = \mathbf{y}$ 
$$
\mathbf{U}=\left[\begin{array}{cccc}
u_{11} & u_{12} & . . & u_{1 n} \\
0 & u_{22} & . . & u_{2 n} \\
. . & . . & . . & . . \\
0 & . . & 0 & u_{n n}
\end{array}\right]
\cdot\left[\begin{array}{l}
x_{1} \\
x_{2} \\
\ldots \\
x_{n}
\end{array}\right]=\left[\begin{array}{l}
y_{1} \\
y_{2} \\
\ldots \\
y_{n}
\end{array}\right]$$

uzyskujemy rozwiązanie $\mathbf{x}$

$$\left\{\begin{array}{l}
x_{n}=\frac{y_{n}}{u_{n n}} \\
x_{i}=\frac{y_{i}-\sum_{k=i+1}^{n} u_{i k} x_{k}}{u_{i i}} \quad i=n-1, n-2, \ldots, 1
\end{array}\right.$$



===== Zadanie 10 =====

Treść zadania znajduje się w Moodle pod adresem [[https://moodle.umk.pl/WFAIIS/mod/assign/view.php?id=4880|Zadanie 10]]

===== Dodatkowe ćwiczenia =====


**Ćw. 1** 
Wiedząc, że $\mathbf{A}\mathbf{A}^{-1} = \mathbf{I}$ napisz program, który wykorzystując rozkład LU wyznaczy macierz odwrotną $\mathbf{A}^{-1}$ dowolnej (podanej przez użytkownika) macierzy kwadratowej $\mathbf{A}$. Rozwiązanie zadania sprowadza się do rozwiązanie $k$ układów równań, gdzie wektorem wyrazów wolnych jest wektor posiadający $1$ na $k$-tej pozycji (kolumna macierzy jednostkowej) a rozwiązaniem jest $k$-ta kolumna macierzy odwrotnej.  

**Ćw. 2**

Metodą rozkładu LU można uprościć obliczenia wyznacznika, gdyż
$$\det{(\mathbf{A})} = \det{(\mathbf{LU})} = \det{(\mathbf{L})}\det{(\mathbf{U})} = 1 \cdot\det{(\mathbf{U})} = u_{11} \cdot u_{22} \cdot \ldots \cdot u_{nn}$$
Napisz program, który wykorzystując rozkład LU obliczy wyznacznik dowolnej macierzy kwadratowej $\mathbf{A}$.
Współczynniki macierzy $\mathbf{A}$ podaje użytkownik.