→ Slide 1

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

→ Slide 2

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]$$

→ Slide 3

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

→ Slide 4

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

→ Slide 5

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}$$

→ Slide 6

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]$$

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}}$$

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}$$

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}$$

→ Slide 7

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}$

→ Slide 8

Macierze: Przykład w 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");
   }
}

→ Slide 9

Mnożenie macierzy: Przykład w C

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

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];
         }
      }
   }
}

→ Slide 10

Rozkład LU: Przykład w 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];
      }
   }
}

Przykładowy kod lu.c

→ Slide 11

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.$$

→ Slide 12

Zadanie 10

Treść zadania znajduje się w Moodle pod adresem Zadanie 10

→ Slide 13

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.