Serwis Edukacyjny
w I-LO w Tarnowie
obrazek

Materiały dla uczniów liceum

  Wyjście       Spis treści       Wstecz       Dalej  

obrazek

Autor artykułu: mgr Jerzy Wałaszek

©2019 mgr Jerzy Wałaszek
I LO w Tarnowie

obrazek

Równania

Rozkład LU

SPIS TREŚCI
Podrozdziały

Definicje

Macierz trójkątna (ang. triangular matrix) jest specjalnym rodzajem macierzy kwadratowej. Występuje ona w dwóch odmianach:

Macierz trójkątna dolna  (ang. lower triangular matrix, L) posiada wszystkie wyrazy leżące ponad główną przekątna równe zero:

Macierz trójkątna górna (ang. upper triangular matrix, U) posiada wszystkie wyrazy leżące poniżej przekątnej głównej równe zero:

Rozkład LU (ang. LU decomposition) polega na rozłożeniu macierzy kwadratowej A na iloczyn dwóch macierzy trójkątnych, dolnej L o wyrazach przekątnej głównej równych 1 oraz górnej U:

Zgodnie z własnościami wyznaczników, wyznacznik iloczynu macierzy jest równy iloczynowi ich wyznaczników:

Z kolei wyznacznik macierzy trójkątnej jest równy iloczynowi wyrazów na przekątnej głównej, co wynika bezpośrednio z rozwinięcia Laplace'a. Zatem:

Wyznacznik macierzy A jest równy wyznacznikowi macierzy U. Z kolei wyznacznik macierzy U, która jest macierzą trójkątną, jest równy iloczynowi wyrazów na głównej przekątnej macierzy U:

I ostatecznie:

Wynika z tego, iż jeśli uda nam się znaleźć macierz U dla macierzy A, to wyznacznik macierzy A policzymy w prosty sposób przez wymnożenie przez siebie wyrazów głównej przekątnej macierzy U. Zajmijmy się zatem sposobem wyznaczenia macierzy U.

Na początek:  podrozdziału   strony 

Rozkład LU Doolitle'a

Ponieważ iloczyn macierzy L i U musi 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 a11 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 aij:

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

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

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

Indeksy: i > j

Wyraz aij leży pod główną przekątną macierzy A. Ostatnim potencjalnie niezerowym elementem w kolumnie j-tej w macierzy U jest element ujj. 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 lij 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 uij wg wzoru:

   

    Dla i od j + 1 do n wyznaczamy lij 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 ujj. Jeśli ujj = 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(n3), 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
Przykładowy program w języku C++
// Macierze
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//---------------------------

#include <iostream>
#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:
            T * A;                   // Wyrazy macierzy
        public:
            int n;                   // Stopień macierzy
            matrix();                // Konstruktor
            matrix(int nn);          // Konstruktor
            ~matrix();               // Destruktor
            T  getv(int r, int c);
            void setv(int r, int c, T v);
            void set1();             // ustawia elementy przekątnej na 1
            T det();                 // Oblicza wyznacznik
    };

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

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

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

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

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

    // Ustawia przekatną na 1
    //-----------------------
    template <class T> void matrix<T>::set1()
    {
        int i;
        for(i = 0; i < n; i++) setv(i,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 sum;
        for(j = 0; j < n; j++)
        {
            for(i = 0; i <= j; i++)
            {
                sum = (T) 0;
                for(k = 0; k < i; k++) sum += L.getv(i,k) * U.getv(k,j);
                U.setv(i,j,getv(i,j) - sum);
            }
            for(i = j + 1; i < n; i++)
            {
                sum = (T) 0;
                for(k = 0; k < j; k++) sum += L.getv(i,k) * U.getv(k,j);
                if(U.getv(j,j)) L.setv(i,j,(getv(i,j) - sum) / U.getv(j,j));
            }
        }
        sum = (T) 1;
        for(i = 0; i < n; i++) sum *= U.getv(i,i);
        return sum;
    }
#endif // _matrix_class

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

int main()
{

    // Wczytujemy dane dla macierzy

    setlocale(LC_ALL,"");

    cout << setprecision(2) << fixed;

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

    matrix<double> A;

    cout << endl
         << "Wyznacznik macierzy metodą rozkładu LU Doolitle'a" << endl
         << "-------------------------------------------------" << endl << endl;
    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 << "det A = " << A.det();

    cout << endl << endl;

    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

det A = -9368.00

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
//---------------------------

#include <iostream>
#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:
            T * A;                   // Wyrazy macierzy
            T eps;
        public:
            int n;                   // Stopień macierzy
            bool error;              // Błąd obliczania wyznacznika
            matrix();                // Konstruktor
            ~matrix();               // Destruktor
            T  getv(int r, int c);
            void setv(int r, int c, T v);
            T det();                 // Oblicza wyznacznik
    };

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

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

    // 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 sum;
        for(j = 0; j < n; j++)
        {
            if((getv(j,j) > - eps) && (getv(j,j) < eps))
            {
                error = true; // Błąd w obliczaniu wyznacznika
                return 0;
            }
            for(i = 0; i <= j; i++)
            {
                sum = (T) 0;
                for(k = 0; k < i; k++) sum += getv(i,k) * getv(k,j);
                setv(i,j,getv(i,j) - sum);
            }
            for(i = j + 1; i < n; i++)
            {
                sum = (T) 0;
                for(k = 0; k < j; k++) sum += getv(i,k) * getv(k,j);
                setv(i,j,(getv(i,j) - sum) / getv(j,j));
            }
        }
        sum = (T) 1;
        for(i = 0; i < n; i++) sum *= getv(i,i);
        return sum;
    }
#endif // _matrix_class

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

int main()
{

    // Wczytujemy dane dla macierzy

    setlocale(LC_ALL,"");

    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" << endl
         << "-----------------------------------------------------------" << endl << endl;
    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 << "det A = " << A.det();

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

    cout << endl << endl;

    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

det A = -9368.00

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:

  1. Normalizacji kolumny – elementy w kolumnie pod elementem akk są dzielone przez ten element.
  2. 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 aij tej podmacierzy iloczynu elementów aik × akj.

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
//---------------------------

#include <iostream>
#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:
            T * A;                   // Wyrazy macierzy
            T eps;
        public:
            int n;                   // Stopień macierzy
            bool error;              // Błąd obliczania wyznacznika
            matrix();                // Konstruktor
            ~matrix();               // Destruktor
            T  getv(int r, int c);
            void setv(int r, int c, T v);
            T det();                 // Oblicza wyznacznik
    };

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

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

    // 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 sum, 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;
        sum = getv(0,0);
        for(k = 1; k < n; k++) sum *= getv(k,k);
        return sum;
    }
#endif // _matrix_class

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

int main()
{

    // Wczytujemy dane dla macierzy

    setlocale(LC_ALL,"");

    cout << setprecision(2) << fixed;

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

    matrix<double> A;

    cout << endl
         << "Wyznacznik macierzy metodą normalizacji kolumn" << endl
         << "----------------------------------------------" << endl << endl;
    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 << "det A = " << A.det();

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

    cout << endl << endl;

    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

det A = -9368.00
Na początek:  podrozdziału   strony 

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

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
Przykładowy program w języku C++
// Macierze
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//---------------------------

#include <iostream>
#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:
            T * A;                   // Wyrazy macierzy
            T eps;
        public:
            int n;                   // Stopień macierzy
            bool error;              // Błąd obliczania wyznacznika
            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 cin
    //-----------------------------------------------
    template <class T> matrix<T>::matrix()
    {
        error = false;
        eps = (T) 0.000000001;
        cin >> n;       // Stopień macierzy
        int nn = n * n;
        A = new T [nn]; // Rezerwujemy pamięć na macierz n x n

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

    // 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 < 0) x = - x;
        return x;
    }

    // Oblicza wyznacznik macierzy
    template <class T> T matrix<T>::det()
    {
        int i, j, k, * W, md, maxw;
        T sum, 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++)
        {
            maxw = k; // numer wiersza z elementem max
            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 = false;
                return 0;
            }

            // Jeśli max nie leży na przekątnej, zamieniamy wiersze
            if(maxw != k)
            {
                md = -md; // Modyfikujemy mnożnik
                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;
        sum = md * getv(W[0],0);
        for(k = 1; k < n; k++) sum *= getv(W[k],k);
        delete [] W;
        return sum;
    }
#endif // _matrix_class

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

int main()
{

    // Wczytujemy dane dla macierzy

    setlocale(LC_ALL,"");

    cout << setprecision(2) << fixed;

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

    matrix<double> A;

    cout << endl
         << "Wyznacznik metodą LU z częściowym wyborem elementu głównego" << endl
         << "-----------------------------------------------------------" << endl << endl;
    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 << "det A = " << A.det();

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

    cout << endl << endl;

    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

det A = -8598.00
Na początek:  podrozdziału   strony 

Zespół Przedmiotowy
Chemii-Fizyki-Informatyki

w I Liceum Ogólnokształcącym
im. Kazimierza Brodzińskiego
w Tarnowie
ul. Piłsudskiego 4
©2019 mgr Jerzy Wałaszek

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.