Metody Numeryczne - Rozkład LU

Ponieważ iloczyn macierzy L i Umusi dać w wyniku macierz A, to każdy wyraz macierzy A jest sumą iloczynów wyrazów z odpowiedniego wiersza macierzy L i odpowiedniej kolumny macierzy U:

Na przykład wyraz a₁₁ jest sumą iloczynów kolejnych wyrazów z wiersza nr 1 macierzy L przez kolejne wyrazy z kolumny nr 1 macierzy U:

gdzie n jest stopniem macierzy.

Ponieważ obie macierze L i U są macierzami trójkątnymi, to część ich elementów ma wartość zero. Np. w macierzy L w wierszu nr 1 tylko pierwszy wyraz jest niezerowy i równy 1, pozostałe wyrazy tego wiersza są równe 0, zatem:

Pozwala to wyznaczyć pierwszy wyraz macierzy U, który jest taki sam jak wyraz macierzy A.

Rozważmy teraz dwa możliwe przypadki dla wyrazu a_ij:

Indeksy: i ≤ j

Wyraz a_ij leży na głównej przekątnej macierzy A lub ponad nią. Ostatnim niezerowym elementem w macierzy L w wierszu i-tym jest l_ii. Wzór iloczynowy możemy więc zredukować:

Co więcej, element l_ii jest elementem przekątnej głównej macierzy L i ma wartość 1, zatem:

Wzór ten pozwala wyznaczyć element u_ij, jeśli znamy poprzednie elementy macierzy L i U:

indeksy: i > j

Wyraz a_ij leży pod główną przekątną macierzy A. Ostatnim potencjalnie niezerowym elementem w kolumnie j-tej w macierzy U jest element u_jj. Wszystkie pozostałe elementy macierzy U w tej kolumnie mają wartość zero. Wzór iloczynowy redukuje się do:

Wyprowadźmy poza sumę jej ostatni wyraz:

Wzór ten pozwala wyliczyć wyraz l_ij macierzy L, jeśli znane są pozostałe wyrazy macierzy L i U:

Otrzymaliśmy wzory na elementy obu macierzy L i U rozkładu LU macierzy A. Wzory te wykorzystuje znany algorytm Doolitle'a, który wygląda następująco:

Tworzymy dwie macierze L i U o stopniu n macierzy A. Wypełniamy je zerami. W macierzy L główną przekątną wypełniamy elementami o wartości 1.

Dla j od 1 do n wykonujemy dwie poniższe operacje:

Dla i od 1 do j wyznaczamy u_ij wg wzoru:

Dla i od j = 1 do n wyznaczamy l_ij wg wzoru:

Uwaga: obliczenia powinny być wykonywane na liczbach zmiennoprzecinkowych.

Przyjrzyj się drugiemu ze wzorów w algorytmie Doolitle'a. Występuje tutaj dzielenie przez element u_jj. Jeśli u_jj = 0, dzielenie jest niewykonalne i nie można wyznaczyć rozkładu LU macierzy A. Problemem tym zajmiemy się dalej w rozdziale.

Metoda Doolitle'a ma złożoność obliczeniową klasy O(n³), co sprawia iż nadaje się ona do obliczeń wyznaczników nawet dużych macierzy.

Poniższy program wyznacza rozkład LU macierzy wejściowej A, a następnie oblicza jej wyznacznik mnożąc przez siebie elementy głównej przekątnej macierzy U. Macierz A przekazywana jest przez standardowe wejście w następującej postaci:

pierwsza liczba określa stopień macierzy, następne n × n liczb to zawartość kolejnych wierszy macierzy A. Program wyświetla odczytaną macierz oraz wartość jej wyznacznika.

Przykładowe dane (przekopiuj je do schowka i wklej w programie):

Macierz:

Dane wejściowe: 5 5 3 7 4 2 9 2 2 1 1 3 6 2 8 9 9 4 -2 -1 -3 0 5 3 -6 -11

C++

// Macierze
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 0061
//---------------------------

#include <iostream>
#include <windows.h>
#include <iomanip>

using namespace std;

// Definicja klasy,
// może być w pliku
// nagłówkowym
//-----------------

#ifndef _matrix_class
  #define _matrix_class

template <class T>
          class matrix
{
  private:
    // Wyrazy macierzy
    T * A;
  public:
    // Stopień macierzy
    int n;
    // Konstruktor
    matrix();
    // Konstruktor
    matrix(int nn);
    // Destruktor
    ~matrix();
    T  getv(int r, int c);
    void setv(int r, int c, T v);
    // ustawia elementy
    // przekątnej na 1
    void set1();
    // Oblicza wyznacznik
    T det();
};

// Konstruktor #1 tworzy
// macierz ze strumienia cin
//--------------------------
template <class T>
         matrix<T>::matrix()
{
  // Stopień macierzy
  cin >> n;
  // wielkość macierzy
  int nn = n * n;
  // Rezerwujemy pamięć
  // na macierz n x n
  A = new T [nn];

  // Wczytujemy wyrazy macierzy
  for(int i = 0; i < nn; i++)
    cin >> A[i];
}

// Konstruktor #2 tworzy macierz
// stopnia n i wypełnia ją zerami
//-------------------------------
template <class T>
         matrix<T>::matrix(int nn)
{
  // Stopień macierzy
  n = nn;
  nn = n * n;
  // Rezerwujemy pamięć
  // na macierz n x n
  A = new T [nn];

  for(int i = 0; i < nn; i++)
    A[i] = (T)0;
}

// Destruktor usuwa macierz
// i zwraca zajętą przez
// nią pamięć do systemu
//-------------------------
template <class T>
         matrix<T>::~matrix()
{
  delete [] A;
}

// Zwraca wartość elementu
// r - numer wiersza
// c - numer kolumny
//-------------------------
template <class T>
         T matrix<T>
         ::getv(int r, int c)
{
  return A[r * n + c];
}

// Ustawia wartość elementu
// r - numer wiersza
// c - numer kolumny
// v - wartość elementu
//-------------------------
template <class T>
         void matrix<T>
         ::setv(int r, int c, T v)
{
  A[r * n + c] = v;
}

// Ustawia przekatną na 1
//-----------------------
template <class T>
         void matrix<T>::set1()
{
  for(int i = 0; i < n; i++)
    A[i * n + i] = 1;
}

// Oblicza wyznacznik
// macierzy
template <class T>
         T matrix<T>::det()
{
  matrix<T> L(n);
  matrix<T> U(n);
  L.set1();
  int i,j,k;
  T sc;
  for(j = 0; j < n; j++)
  {
    for(i = 0; i <= j; i++)
    {
      sc = (T)0;
      for(k = 0; k < i; k++)
        sc += L.getv(i,k)
              * U.getv(k,j);
      U.setv(i,j,getv(i,j) - sc);
    }
    for(i = j + 1; i < n; i++)
    {
      sc = (T)0;
      for(k = 0; k < j; k++)
        sc += L.getv(i,k)
              * U.getv(k,j);
      if(U.getv(j,j))
        L.setv(i,j,(getv(i,j) - sc)
               / U.getv(j,j));
    }
  }
  sc = (T)1;
  for(i = 0; i < n; i++)
    sc *= U.getv(i,i);
  return sc;
}

#endif // _matrix_class

// Program główny
//---------------
int main()
{

  // Wczytujemy dane
  // dla macierzy

  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(2)
       << fixed;

  cout
  << "Wpisz dane dla macierzy:\n\n";

  matrix<double> A;

  cout
  << "\nWyznacznik macierzy metodą"
     "rozkładu LU Doolitle'a\n"
     "--------------------------"
     "----------------------\n\n";
  for(int i = 0; i < A.n; i++)
  {
    for(int j = 0; j < A.n; j++)
      cout << setw(8)
           << A.getv(i,j);
    cout << endl;
  }
  cout << "\n\nWyznacznik = "
       << A.det();

  cout << endl << endl;
  system("pause");
  return 0;
}

Wynik

Wpisz dane dla macierzy:

5
5 3 7 4 2
9 2 2 1 1
3 6 2 8 9
9 4 -2 -1 -3
0 5 3 -6 -11

Wyznacznik macierzy metodą rozkładu LU Doolitle'a
-------------------------------------------------

    5.00    3.00    7.00    4.00    2.00
    9.00    2.00    2.00    1.00    1.00
    3.00    6.00    2.00    8.00    9.00
    9.00    4.00   -2.00   -1.00   -3.00
    0.00    5.00    3.00   -6.00 -11.00

Wyznacznik = -9368.00

Python (dodatek)

# Macierze
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0061
#---------------------------

# Definicja klasy

class Matrix:

    # Konstruktor tworzy
    # macierz ze strumienia
    # wejściowego lub z wymiaru
    # i wypełnia ją zerami
    #--------------------------
    def __init__(self, nn=None):
        if nn is None:
            # Czytamy stopień
            # macierzy
            self.n = int(input())
            
            # Wczytujemy kolejne
            # wiersze i zmieniamy je
            # w listy liczb float
            self.a = [[float(x)
             for x in input().split()]
             for _ in range(self.n)]
        else:
            # Tworzymy macierz  
            # stopnia nn
            # i wypełniamy ją zerami
            #-----------------------
            self.n = nn
            self.a = [[0.0 
             for _ in range(self.n)]
             for _ in range(self.n)]

    # Ustawia przekatną na 1
    #-----------------------
    def set1(self):
        for i in range(self.n):
            self.a[i][i] = 1

    # Oblicza wyznacznik
    # macierzy
    def det(self):
        m_l = Matrix(self.n)
        m_u = Matrix(self.n)
        m_l.set1()
        for j in range(self.n):
            for i in range(j + 1):
                sc = 0
                for k in range(i):
                    sc += (m_l.a[i][k]
                      * m_u.a[k][j])
                m_u.a[i][j] = \
                  self.a[i][j] - sc
            for i in range(j + 1, self.n):
                sc = 0
                for k in range(j):
                    sc += (m_l.a[i][k]
                    * m_u.a[k][j])
                if m_u.a[j][j]:
                    m_l.a[i][j] = (
                    (self.a[i][j] - sc)
                    / m_u.a[j][j])
        sc = 1
        for i in range(self.n):
            sc *= m_u.a[i][i]
        return sc

# Program główny
#---------------

# Wczytujemy dane
# dla macierzy
print("Wpisz dane dla macierzy:\n")

t = Matrix()

print("Wyznacznik macierzy metodą "
      "rozkładu LU Doolitle'a\n"
      "--------------------------"
      "----------------------\n")
for i in range(t.n):
    for j in range(t.n):
        print(f"{t.a[i][j]:8.2f}",
              end="")
    print()

print(f"\nWyznacznik = {t.det():.2f}")

print()
input("Naciśnij Enter...")

Jeśli przyjrzysz się dokładnie wzorom Doolitle'a, to powinieneś zauważyć, że wyrazy macierzy A są wykorzystywane tylko jeden raz w każdym obiegu wyznaczania wyrazów macierzy L i U. Pozwala to umieścić macierze L i U w macierzy A, przez co algorytm działa w miejscu, tzn. nie zużywa dodatkowej pamięci. Jest to możliwe, ponieważ przekątna główna macierzy dolnej L nie jest używana przez algorytm. Zatem obie macierze L i U wypełnią dokładnie macierz A.

Poniższy program wykorzystuje tę koncepcję do wyliczenia wyznacznika. W przypadku próby dzielenia przez zero ustawiany jest znacznik error na true w macierzy A. Pozwala to programowi sprawdzić, czy obliczony wynik jest poprawny.

Przykładowe dane (przekopiuj je do schowka i wklej w programie) :

Macierz:

Dane wejściowe: 5 5 3 7 4 2 9 2 2 1 1 3 6 2 8 9 9 4 -2 -1 -3 0 5 3 -6 -11

Przykładowy program w języku C++

// Macierze
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 0062
//---------------------------

#include <iostream>
#include <windows.h>
#include <iomanip>

using namespace std;

// Definicja klasy,
// może być w pliku
// nagłówkowym
//-----------------

#ifndef _matrix_class
  #define _matrix_class

template <class T>
         class matrix
{
  private:
    // Wyrazy macierzy
    T * A;
    T EPS;
  public:
    // Stopień macierzy
    int n;
    // Wielkość macierzy
    int s;
    // Błąd obliczenia wyznacznika
    bool error;
    matrix();  // Konstruktor
    ~matrix(); // Destruktor
    T  getv(int r, int c);
    void setv(int r, int c, T v);
    // Oblicza wyznacznik
    T det();
};

// Konstruktor tworzy macierz
// ze strumienia cin
//---------------------------
template <class T>
         matrix<T>::matrix()
{
  error = false;
  EPS = (T)0.000000001;

  // Stopień macierzy
  cin >> n;

  // Wielkość macierzy
  s = n * n;

  // Rezerwujemy pamięć
  // na macierz n x n
  A = new T [s];

  // Wczytujemy wyrazy
  // macierzy
  for(int i = 0; i < s; i++)
    cin >> A[i];
}

// Destruktor usuwa macierz
// i zwraca zajętą przez
// nią pamięć
//-------------------------
template <class T>
         matrix<T>::~matrix()
{
  delete [] A;
}

// Zwraca wartość elementu
// r - numer wiersza
// c - numer kolumny
//-------------------------
template <class T>
         T matrix<T>
         ::getv(int r, int c)
{
  return A[r * n + c];
}

// Ustawia wartość elementu
// r - numer wiersza
// c - numer kolumny
// v - wartość elementu
//-------------------------
template <class T>
         void matrix<T>
         ::setv(int r, int c, T v)
{
  A[r * n + c] = v;
}

// Oblicza wyznacznik
// macierzy
template <class T>
         T matrix<T>::det()
{
  int i,j,k;
  T sc;
  error = false;
  for(j = 0; j < n; j++)
  {
    if((getv(j,j) > - EPS) &&
       (getv(j,j) < EPS))
    {
      // Błąd w obliczaniu
      // wyznacznika
      error = true;
      return 0;
    }
    for(i = 0; i <= j; i++)
    {
      sc = (T)0;
      for(k = 0; k < i; k++)
        sc += getv(i,k)* getv(k,j);
      setv(i,j,getv(i,j) - sc);
    }
    for(i = j + 1; i < n; i++)
    {
      sc = (T)0;
      for(k = 0; k < j; k++)
        sc += getv(i,k) * getv(k,j);
      setv(i,j,(getv(i,j) - sc) /
           getv(j,j));
    }
  }
  sc = (T)1;
  for(i = 0; i < n; i++)
    sc *= getv(i,i);
  return sc;
}

#endif // _matrix_class

// Program główny
//---------------
int main()
{
  // Wczytujemy dane
  // dla macierzy

  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(2)
       << fixed;

  cout << "Wpisz dane dla macierzy:"
       << endl << endl;

  matrix<double> A;

  cout << endl
       << "Wyznacznik macierzy metodą "
          "rozkładu LU Doolitle'a "
          "w miejscu\n"
          "---------------------------"
          "-----------------------"
          "---------\n\n";
  for(int i = 0; i < A.n; i++)
  {
    for(int j = 0; j < A.n; j++)
      cout << setw(8) << A.getv(i,j);
    cout << endl;
  }

  cout << "\nWyznacznik = " << A.det();

  if(A.error)
    cout << " - BŁĄD W OBLICZENIACH !!!";

  cout << endl << endl;
  system("pause");
  return 0;
}

Wynik

Wpisz dane dla macierzy:

5
5 3 7 4 2
9 2 2 1 1
3 6 2 8 9
9 4 -2 -1 -3
0 5 3 -6 -11

Wyznacznik macierzy metodą rozkładu LU Doolitle'a w miejscu
-----------------------------------------------------------

    5.00    3.00    7.00    4.00    2.00
    9.00    2.00    2.00    1.00    1.00
    3.00    6.00    2.00    8.00    9.00
    9.00    4.00   -2.00   -1.00   -3.00
    0.00    5.00    3.00   -6.00 -11.00

Wyznacznik = -9368.00

Python (dodatek)

# Macierze
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0062
#---------------------------

# Definicja klasy,
#-----------------
class Matrix:
        
    # Konstruktor tworzy
    # macierz ze strumienia
    # wejścia
    #----------------------
    def __init__(self):
        # Zerujemy błąd
        self.error = False
    
        # Dokładność
        self.EPS = 0.000000001

        # czytamy stopień macierzy
        self.n = int(input())

        # Wczytujemy kolejne
        # wiersze i zmieniamy je
        # w listy liczb float
        self.a = [[float(x)
         for x in input().split()]
         for _ in range(self.n)]

    # Oblicza wyznacznik macierzy
    def det(self):
        # Zerujemy flagę błędu
        self.error = False

        for j in range(self.n):
            if (abs(self.a[j][j]) <
                self.EPS):
                # Błąd w obliczaniu
                # wyznacznika
                self.error = True
                return 0
            for i in range(j+1):
                sc = sum(self.a[i][k]
                 * self.a[k][j]
                 for k in range(i))
                self.a[i][j] = \
                 self.a[i][j] - sc
            for i in range(j+1,self.n):
                sc = sum(self.a[i][k]
                 * self.a[k][j]
                 for k in range(j))
                self.a[i][j] = (
                 (self.a[i][j] - sc) 
                 / self.a[j][j])
        sc = 1
        for i in range(self.n):
            sc *= self.a[i][i]
        return sc

# Program główny
#---------------

# Wczytujemy dane
# dla macierzy

print("Wpisz dane dla macierzy:")
print()

t = Matrix()

print()
print("Wyznacznik macierzy metodą "
      "rozkładu LU Doolitle'a "
      "w miejscu\n"
      "---------------------------"
      "-----------------------"
      "---------\n")

# Wyświetlamy macierz
for i in range(t.n):
    for j in range(t.n):
        print(f"{t.a[i][j]:8.2f}",
              end=" ")
    print()

print(f"\nWyznacznik = {t.det():.2f}",
      end="")

if t.error:
    print(" - BŁĄD W OBLICZENIACH!!!")

print()
input("Naciśnij Enter...")

Istnieje prostsze rozwiązanie problemu wyznaczania rozkładu LU macierzy. Polega ono na przechodzeniu przez kolejne elementy głównej przekątnej macierzy A od pierwszego do przedostatniego i dla każdego elementu A[k,k] Wykonywaniu kolejno dwóch operacji:

Normalizacji kolumny – elementy w kolumnie pod elementem a_kk są dzielone przez ten element.
Modyfikacji podmacierzy, która zawiera wszystkie elementy macierzy A leżące poniżej i na prawo elementu przekątnej z pominięciem wiersza i kolumny zawierającego element przekątnej. Modyfikacja polega na odjęciu od każdego elementu a_ij tej podmacierzy iloczynu elementów a_ik × a_kj.

Dzięki takiemu podejściu algorytm rozkładu LU znacznie się upraszcza.

Poniższy program wykorzystuje ten algorytm do wyznaczenia rozkładu LU macierzy. Macierz jest wprowadzana tak samo, jak w programie poprzednim. Po dokonaniu rozkładu LU program oblicza wyznacznik macierzy i wyświetla wynik. Jeśli w trakcie pracy program napotka zerowy element na głównej przekątnej, to zakończy działanie z informacją o błędzie.

Przykładowe dane (przekopiuj je do schowka i wklej w programie) :

Macierz:

Dane wejściowe: 5 5 3 7 4 2 9 2 2 1 1 3 6 2 8 9 9 4 -2 -1 -3 0 5 3 -6 -11

Przykładowy program w języku C++

// Macierze
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 0063
//---------------------------

#include <iostream>
#include <windows.h>
#include <iomanip>

using namespace std;

// Definicja klasy,
// może być w pliku
// nagłówkowym
//-----------------

#ifndef _matrix_class
  #define _matrix_class

template <class T>
         class matrix
{
  private:
    // Wyrazy macierzy
    T * A;
    T EPS;
  public:
    // Stopień macierzy
    int n;
    // Błąd obliczenia
    // wyznacznika
    bool error;
    // Konstruktor
    matrix();
    // Destruktor
    ~matrix();
    T getv(int r, int c);
    void setv(int r, int c, T v);
    // Oblicza wyznacznik
    T det();
};

// Konstruktor tworzy
// macierz ze strumienia
// wejścia cin
//----------------------
template <class T>
         matrix<T>::matrix()
{
  // Zerujemy błąd
  error = false;

  // Dokładność
  EPS = (T)0.000000001;

  // Stopień macierzy
  cin >> n;

  // Rozmiar macierzy
  int nn = n * n;

  // Rezerwujemy pamięć
  // na macierz n x n
  A = new T [nn];

  // Wczytujemy kolejne
  // wyrazy macierzy
  for(int i = 0; i < nn; i++)
    cin >> A[i];
}

// Destruktor usuwa macierz
// i zwraca zajętą przez
// nią pamięć
//-------------------------
template <class T>
         matrix<T>::~matrix()
{
  delete [] A;
}

// Zwraca wartość elementu
// r - numer wiersza
// c - numer kolumny
//-------------------------
template <class T>
         T matrix<T>
         ::getv(int r, int c)
{
  return A[r * n + c];
}

// Ustawia wartość elementu
// r - numer wiersza
// c - numer kolumny
// v - wartość elementu
//-------------------------
template <class T>
         void matrix<T>
         ::setv(int r, int c,
                T v)
{
  A[r * n + c] = v;
}

// Oblicza wyznacznik
// macierzy
template <class T>
         T matrix<T>::det()
{
  int i,j,k;
  T sc, akk;

  for(k = 0; k < n - 1; k++)
  {
    // Dzielnik
    akk = getv(k,k);
    if((akk > -EPS) && (akk < EPS))
    {
      error = true;
      return 0;
    }

    // Normalizujemy kolumnę
    for(i = k + 1; i < n; i++)
      setv(i,k,getv(i,k) / akk);

    // Modyfikujemy podmacierz
    for(i = k + 1; i < n; i++)
      for(j = k + 1; j < n; j++)
        setv(i,j,
             getv(i,j) - getv(i,k)
             * getv(k,j));
  }

  // Obliczamy wyznacznik;
  sc = getv(0,0);
  for(k = 1; k < n; k++)
    sc *= getv(k,k);
  return sc;
}

#endif // _matrix_class

// Program główny
//---------------
int main()
{
  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  // Wczytujemy dane
  // dla macierzy
  cout << setprecision(2)
       << fixed;
  cout << "Wpisz dane dla macierzy:"
       << endl << endl;

  matrix<double> A;

  cout << endl
       << "Wyznacznik macierzy "
          "metodą normalizacji "
          "kolumn\n"
          "--------------------"
          "--------------------"
          "------\n\n";
  for(int i = 0; i < A.n; i++)
  {
    for(int j = 0; j < A.n; j++)
      cout << setw(8)<< A.getv(i,j);
    cout << endl;
  }

  cout << endl
       << "Wyznacznik = "
       << A.det();

  if(A.error)
    cout << " - BŁĄD OBLICZENIA !!!";

  cout << endl << endl;
  system("pause");
  return 0;
}

Wynik

Wpisz dane dla macierzy:

5
5 3 7 4 2
9 2 2 1 1
3 6 2 8 9
9 4 -2 -1 -3
0 5 3 -6 -11

Wyznacznik macierzy metodą normalizacji kolumn
----------------------------------------------

    5.00    3.00    7.00    4.00    2.00
    9.00    2.00    2.00    1.00    1.00
    3.00    6.00    2.00    8.00    9.00
    9.00    4.00   -2.00   -1.00   -3.00
    0.00    5.00    3.00   -6.00  -11.00

Wyznacznik = -9368.00

Python (dodatek)

# Macierze
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0063
#---------------------------

# Definicja klasy
#----------------
class Matrix:

    # Konstruktor tworzy
    # macierz ze strumienia
    # wejścia
    #----------------------
    def __init__(self):
        # Zerujemy błąd
        self.error = False
    
        # Dokładność
        self.EPS = 0.000000001

        # czytamy stopień macierzy
        self.n = int(input())

        # Wczytujemy kolejne
        # wiersze i zmieniamy je
        # w listy liczb float
        self.a = [[float(x)
         for x in input().split()]
         for _ in range(self.n)]

    # Oblicza wyznacznik
    # macierzy
    def det(self):
        for k in range(self.n - 1):
            # Dzielnik
            akk = self.a[k][k]
            if abs(akk) < self.EPS:
                self.error = True
                return 0

            # Normalizujemy kolumnę
            for i in range(k+1,self.n):
                self.a[i][k] = \
                 self.a[i][k] / akk

            # Modyfikujemy podmacierz
            for i in range(k+1,self.n):
                for j in range(k+1,self.n):
                    self.a[i][j] = \
                     self.a[i][j] - \
                     self.a[i][k] * \
                     self.a[k][j]

        # Obliczamy wyznacznik
        sc = self.a[0][0]
        for k in range(1, self.n):
            sc *= self.a[k][k]
        return sc

# Program główny
#---------------

# Wczytujemy dane dla macierzy
print("Wpisz dane dla macierzy:")
print()

t = Matrix()

print()
print("Wyznacznik macierzy "
      "metodą normalizacji "
      "kolumn\n"
      "--------------------"
      "--------------------"
      "------\n")

# Wyświetlamy macierz
for i in range(t.n):
    for j in range(t.n):
        print(f"{t.a[i][j]:8.2f}",
              end="")
    print()

print()
print(f"Wyznacznik = {t.det():.2f}",
      end="")

if t.error:
    print(" - BŁĄD OBLICZENIA !!!")

print()
input("Naciśnij Enter...")

do podrozdziału do strony

Podane powyżej metody rozkładu LU zawodzą w przypadku napotkania elementu zerowego przez który należy wykonać dzielenie. W takich przypadkach programy przerywały obliczenia. Tymczasem rozwiązanie istnieje i jest dosyć proste. Wystarczy w takim przypadku tak poprzestawiać wiersze rozkładanej macierzy, aby na jej przekątnej znalazły się wyrazy o największych modułach. Wprowadźmy pojęcie macierzy permutacji P. Jest to macierz jednostkowa o stopniu macierzy A, w której zamieniono miejscami wiersze lub kolumny. Jeśli pomnożymy macierz A przez tak zmienioną macierz permutacji P, to wynikowa macierz będzie miała przestawione wiersze lub kolumny tak samo, jak zostały one przestawione w macierzy P. Obowiązuje zasada: mnożenie P × A przestawia wiersze, A × P przestawia kolumny. Zobaczmy na przykładzie, jak to działa:

Mamy macierz A czwartego stopnia:

Tworzymy dla niej macierz permutacji P:

W takiej postaci permutacja nie zamienia miejscami wierszy/kolumn macierzy A:

Zamieńmy miejscami dwa środkowe wiersze (lub dwie środkowe kolumny, co daje identyczny wynik) w macierzy permutacji P:

Teraz mnożenia (sprawdź to) dadzą wyniki:

Pierwsze mnożenie dało w wyniku macierz A z zamienionymi dwoma środkowymi wierszami, a drugie mnożenie dało w wyniku macierz A z zamienionymi dwiema środkowymi kolumnami. Zatem macierz permutacji P ma za zadanie odpowiednio pozamieniać ze sobą miejscami wiersze lub kolumny. Taka zamiana nazywa się piwotem (ang. pivoting). Termin ten pochodzi z gry w kosza, gdzie oznacza obrót zawodnika z piłką.

Po co to wszystko? Otóż permutacje pozwolą nam tak poukładać wiersze (lub kolumny) macierzy, aby na jej przekątnej głównej znalazły się wyrazy o największym module. Jeśli macierz nie jest zdegenerowana, to jest to zawsze możliwe. Permutacje wykonujemy w trakcie wyznaczania rozkładu LU. Musimy jednak zapamiętywać ich liczbę, ponieważ każda zamiana wierszy/kolumn zmienia znak wyznacznika na przeciwny.

Praktycznie nie tworzy się macierzy P. Permutacje wierszy zapamiętujemy w wektorze wierszowym W, poprzez który odwołujemy się do wierszy macierzy A przy wyznaczaniu rozkładu LU. Podobnie postępowaliśmy w trakcie wyliczania wartości wyznacznika za pomocą rozwinięcia Laplace'a. Każdorazowo w pętli rozkładu LU wyznaczamy największy element w pierwszej kolumnie bieżącej podmacierzy i zamieniamy jego wiersz z pierwszym wierszem poprzez wektor kolumnowy. Algorytm jest następujący:

Algorytm rozkładu LU z piwotowaniem

Dane wejściowe:

ε	–	dokładność przyrównania do zera
n	–	stopień macierzy
A	–	macierz podlegająca rozkładowi LU

Dane wyjściowe

A	–	macierz zawierająca połączone macierze L i U
W	–	wektor wierszowy
md	–	mnożnik do korekcji wyznacznika

Zmienne pomocnicze

i,j,k	–	indeksy
maxw	–	numer wiersza z elementem maksymalnym
maxe	–	moduł elementu maksymalnego
akk	–	element macierzy A leżący na jej przekątnej po piwotowaniu

Lista kroków

K01:	md ← 1	mnożnik korekcyjny ustawiamy na 1
K02:	Dla i = 0,1,...,n – 1: wykonuj: W [ i ] ← i	ustawiamy wektor wierszowy
K03:	Dla k = 0,1,...,n – 1: wykonuj kroki K04...K13	rozpoczynamy rozkład LU
K04:	maxw ← k	wiersz z elementem maksymalnym
K05:	maxe ← \| A [ W [ k ], k ] \|	moduł elementu maksymalnego
K06:	Dla i = k + 1, k + 2,..., n – 1: wykonuj krok K07	szukamy elementu większego
K07:	Jeśli \| A [W [ i ], k ] \| > maxe, to maxw ← i maxe ← \| A [ W [ i ], k ] \|
K08:	Jeśli maxe ≤ ε, to zakończ z wynikiem false	macierz jest zdegenerowana
K09:	Jeśli maxw ≠ W [ k ], to md ← - md W [ k ] ↔ W [ maxw ]	znaleziono element maksymalny
K10:	akk ← A [ W [ k ], k ]	rozkład LU
K11:	Dla i = k + 1, k + 2,...,n – 1: wykonuj: A [ W [ i ], k ] ← A [ W [ i ], k ] / akk	normalizacja kolumny
K12:	Dla i = k + 1,k + 2,...,n – 1: wykonuj krok K13	modyfikujemy podmacierz
K13:	Dla j = k + 1,k + 2,...,n – 1: wykonuj: A [ W [ i ], j ] ← A [ W [ i ], j ] - A [ W [ i ], k ] × A [ W [ k ], j ]
K14:	Zakończ z wynikiem true

Przy obliczaniu wyznacznika do wierszy macierzy A należy odwoływać się poprzez wektor wierszowy W.

Podany algorytm posiada jeszcze jedną zaletę: redukuje błędy zaokrągleń, ponieważ gwarantuje dzielenie przez największy element.

Poniższy program wykorzystuje ten algorytm do wyliczania wartości wyznacznika. Jest to modyfikacja poprzedniego programu.

Macierz:

Dane wejściowe: 5 0 3 7 4 2 9 2 2 1 1 3 6 2 8 9 9 4 -2 -1 -3 0 5 3 -6 -11

C++

// Macierze
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 0064
//---------------------------

#include <iostream>
#include <windows.h>
#include <iomanip>

using namespace std;

// Definicja klasy,
// może być w pliku
// nagłówkowym
//-----------------

#ifndef _matrix_class
  #define _matrix_class

template <class T>
         class matrix
{
  private:
    // Wyrazy macierzy
    T * A;
    T EPS;
  public:
    // Stopień macierzy
    int n;
    // Błąd obliczania wyznacznika
    bool error;
    matrix();  // Konstruktor
    ~matrix(); // Destruktor
    T getv(int r, int c);
    void setv(int r, int c, T v);
    T absT(T x);
    T det(); // Oblicza wyznacznik
};

// Konstruktor tworzy
// macierz ze strumienia
// wejścia cin
//----------------------
template <class T>
         matrix<T>::matrix()
{
  error = false;
  EPS = (T)0.000000001;

  // Stopień macierzy
  cin >> n;

  // Wielkość macierzy
  int s = n * n;

  // Rezerwujemy pamięć
  // na macierz n x n
  A = new T [s];

  // Wczytujemy wyrazy
  // macierzy
  for(int i = 0; i < s; i++)
    cin >> A[i];
}

// Destruktor usuwa macierz
// i zwraca zajętą przez
// nią pamięć
//-------------------------
template <class T>
         matrix<T>::~matrix()
{
  delete [] A;
}

// Zwraca wartość elementu
// r - numer wiersza
// c - numer kolumny
//-------------------------
template <class T>
         T matrix<T>
         ::getv(int r, int c)
{
  return A[r * n + c];
}

// Ustawia wartość elementu
// r - numer wiersza
// c - numer kolumny
// v - wartość elementu
//-------------------------
template <class T>
         void matrix<T>
         ::setv(int r, int c, T v)
{
  A[r * n + c] = v;
}

template <class T>
         T matrix<T>::absT(T x)
{
  if(x < (T)0)
    x = - x;
  return x;
}

// Oblicza wyznacznik macierzy
template <class T>
         T matrix<T>::det()
{
  int i, j, k, * W, md, maxw;
  T sc, akk, maxe;

  // Mnożnik korekcyjny
  // dla wyznacznika
  md = 1;

  // Tworzymy wektor wierszowy
  W = new int [n];

  // Inicjujemy wektor wierszowy
  // numerami wierszy
  for(i = 0; i < n; i++)
    W[i] = i;

  // Rozkład LU
  for(k = 0; k < n; k++)
  {
    // numer wiersza
    // z elementem max
    maxw = k;
    maxe = absT(getv(W[k],k));

    // Szukamy elementu max
    for(i = k + 1; i < n; i++)
      if(absT(getv(W[i],k)) > maxe)
      {
        maxw = i;
        maxe = absT(getv(W[i],k));
      }

     // Jeśli max jest zero,
     // to macierz jest
     // zdegenerowana
     if(maxe <= EPS)
     {
       error = true;
       delete [] W;
       return 0;
     }

     // Jeśli max nie leży
     // na przekątnej,
     // zamieniamy wiersze
     if(maxw != k)
     {
       // Modyfikujemy mnożnik
       md = -md;
       swap(W[k], W[maxw]);
     }

     // Dzielnik
     akk = getv(W[k],k);

     // Normalizujemy kolumnę
     for(i = k + 1; i < n; i++)
        setv(W[i],k,getv(W[i],k)/akk);

     // Modyfikujemy podmacierz
     for(i = k + 1; i < n; i++)
       for(j = k + 1; j < n; j++)
         setv(W[i],j,getv(W[i],j)
              -getv(W[i],k)
              *getv(W[k],j));
  }

  // Obliczamy wyznacznik;
  sc = md;
  for(k = 0; k < n; k++)
    sc *= getv(W[k],k);
  delete [] W;
  return sc;
}

#endif // _matrix_class

// Program główny
//---------------

int main()
{
  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(2)
       << fixed;

  // Wczytujemy dane dla macierzy
  cout << "Wpisz dane dla macierzy:"
       << endl << endl;

  matrix<double> A;

  cout
  << endl
  << "Wyznacznik metodą LU "
     "z częściowym wyborem "
     "elementu głównego\n"
     "---------------------"
     "---------------------"
     "-----------------\n\n";

  for(int i = 0; i < A.n; i++)
  {
    for(int j = 0; j < A.n; j++)
      cout << setw(8)<< A.getv(i,j);
    cout << endl;
  }

  cout << endl
       << "Wyznacznik = "
       << A.det();

  if(A.error)
    cout
    << " - BŁĄD W OBLICZENIACH !!!";

  cout << endl << endl;
  system("pause");
  return 0;
}

Wynik

Wpisz dane dla macierzy:

5
0 3 7 4 2
9 2 2 1 1
3 6 2 8 9
9 4 -2 -1 -3
0 5 3 -6 -11

Wyznacznik metodą LU z częściowym wyborem elementu głównego
-----------------------------------------------------------

    0.00    3.00    7.00    4.00    2.00
    9.00    2.00    2.00    1.00    1.00
    3.00    6.00    2.00    8.00    9.00
    9.00    4.00   -2.00   -1.00   -3.00
    0.00    5.00    3.00   -6.00  -11.00

Wyznacznik = -8598.00

Python (dodatek)

# Macierze
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0064
#---------------------------

# Definicja klasy
#----------------
class Matrix:

    # Konstruktor tworzy
    # macierz ze strumienia
    # wejścia cin
    #----------------------
    def __init__(self):
        self.error = False
        self.EPS = 0.000000001

        # Czytamy stopień macierzy
        # z pierwszego wiersza
        self.n = int(input())

        # Wczytujemy n wierszy
        # i tworzymy z nich
        # macierz - listę list
        self.a = [[float(x)
         for x in input().split()]
         for _ in range(self.n)]

    # Oblicza wyznacznik
    # macierzy
    def det(self):
        # Mnożnik korekcyjny
        # dla wyznacznika
        md = 1

        # Tworzymy wektor
        # wierszowy
        w = [i for i in range(self.n)]

        # Rozkład LU
        for k in range(self.n):
            # numer wiersza
            # z elementem max
            maxw = k
            maxe = abs(self.a[w[k]][k])

            # Szukamy elementu max
            for i in range(k+1,self.n):
                if (abs(self.a[w[i]][k])
                    > maxe):
                    maxw = i
                    maxe = \
                     abs(self.a[w[i]][k])

            # Jeśli max jest zerowe,
            # to macierz jest
            # zdegenerowana
            if maxe <= self.EPS:
                # Sygnalizujemy błąd
                self.error = True
                return 0

            # Jeśli max nie leży
            # na przekątnej, to
            # zamieniamy wiersze
            if maxw != k:
                # Modyfikujemy mnożnik
                md = -md
                w[k], w[maxw] = \
                 w[maxw], w[k]

            # Dzielnik
            akk = self.a[w[k]][k]

            # Normalizujemy kolumnę
            for i in range(k+1,self.n):
                self.a[w[i]][k] /= akk

            # Modyfikujemy podmacierz
            for i in range(k+1,self.n):
                for j in range(k+1,self.n):
                    self.a[w[i]][j] -= \
                     self.a[w[i]][k] \
                     *self.a[w[k]][j]

        # Obliczamy wyznacznik;
        sc = md
        for k in range(self.n):
            sc *= self.a[w[k]][k]
    
        return sc

# Program główny
#---------------

# Wczytujemy dane
# dla macierzy
print("Wpisz dane dla macierzy:\n")

t = Matrix()

print()
print("Wyznacznik metodą LU "
      "z częściowym wyborem "
       "elementu głównego\n"
       "---------------------"
       "---------------------"
       "-----------------\n")

for i in range(t.n):
    for j in range(t.n):
        print(f"{t.a[i][j]:8.2f}",
              end="")
    print()

print()
print(f"Wyznacznik = {t.det():.2f}",
      end="")

if t.error:
    print(" - BŁĄD W OBLICZENIACH !!!")

print()
input("Naciśnij Enter...")

do podrozdziału do strony

Materiały tylko do użytku dydaktycznego. Ich kopiowanie i powielanie jest dozwolone pod warunkiem podania źródła oraz niepobierania za to pieniędzy.
Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl
Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.

Informacje dodatkowe.

Równania

Rozkład LU

Definicje

Rozkład LU Doolitle'a

Rozkład LU z częściowym wyborem elementu głównego

Algorytm rozkładu LU z piwotowaniem

Dane wejściowe:

Dane wyjściowe

Zmienne pomocnicze

Lista kroków