Serwis Edukacyjny w I-LO w Tarnowie Materiały dla uczniów liceum |
Wyjście Spis treści Wstecz Dalej Autor artykułu: mgr Jerzy Wałaszek |
©2024 mgr Jerzy Wałaszek |
ProblemDla danej macierzy An × n należy znaleźć macierz odwrotną wykorzystując rozkład LU macierzy A. |
Macierz odwrotna (ang. the inverse of a matrix) (A-1) do macierzy A spełnia równanie:
A × A-1 =
I = A-1
× A gdzie I oznacza macierz jednostkową. |
Oznaczmy X = A-1. Wtedy powyższe równanie możemy zapisać jako:
A × X = I |
Niech:
Weźmy następnie pierwszą kolumnę macierzy X:
Zgodnie z regułami mnożenia macierzy mamy:
W wyniku otrzymujemy pierwszą kolumnę macierzy jednostkowej, gdyż z założenia X = A-1. Wyrazy kolumny X1 możemy potraktować jako niewiadome, a cały iloczyn jako układ n równań liniowych względem tych niewiadomych. Układ ten rozwiązujemy metodą LU, którą opisaliśmy w poprzednim rozdziale. Następnie tak samo postępujemy z drugą kolumną X 2 (rozkład LU wykonujemy tylko jeden raz na początku obliczeń – później wykorzystujemy go do rozwiązania kolejnych układów równań).
Postępując podobnie dla pozostałych kolumn X3, X4 … Xn, wyznaczymy całą macierz X, która jest odwrotnością macierzy A. W algorytmie macierz X może jednocześnie być macierzą jednostkową, której kolejne kolumny przekształcamy w kolumny Xi. Macierz jednostkowa na końcu zostanie więc zastąpiona macierzą odwrotną do A.
n | : | stopień macierzy, n ∈ N |
A | : | macierz kwadratowa stopnia n. Elementy A ∈ R |
r | = | true, macierz kwadratowa stopnia n X zawiera macierz odwrotną do A. Elementy X ∈ R |
r | = | false, błąd dzielenia przez zero |
i | : | indeks, i ∈ N |
K01: | r ← false | zakładamy porażkę |
K02: | Wykonaj rozkład LU dla macierzy A |
|
K03: | Jeśli rozkład się nie powiódł, to zakończ |
|
K04: | Ustaw w X macierz jednostkową | |
K05: | Dla i = 1, 2, …n: wykonuj kroki K06…K07 |
obliczamy kolejne kolumny |
K06: | Przyjmij za wyrazy
wolne
i-tą kolumnę macierzy X. Rozwiąż układ równań metodą LU. Zapisz wektor niewiadomych w i-tej kolumnie X |
macierzy odwrotnej |
K07: | Jeśli
rozwiązanie się nie powiodło, to zakończ |
|
K08: | r ← true | zaznaczamy sukces |
K09: | Zakończ |
Powyższy algorytm nie jest nowym algorytmem, wykorzystuje on dwa poprzednio opisane algorytmy rozkładu LU oraz rozwiązywania układu równań liniowych z wykorzystaniem macierzy L i U otrzymanych w rozkładzie LU. Jedyną różnicą będzie odpowiednie potraktowanie kolumn macierzy X jako wektora wyrazów wolnych B oraz na wyjściu jako wektora niewiadomych. Zadanie to prosto zrealizujemy, modyfikując funkcję lusolve(), tak aby przyjmowała jako parametr numer kolumny macierzy X.
Uwaga:
Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich. |
Program wczytuje dane definiujące
macierz
A i wylicza macierz do niej odwrotną. Jeśli
obliczenie nie może zostać wykonane, program wypisuje komunikat
"DZIELNIK ZERO". Dane dla macierzy A są zdefiniowane następująco: Pierwsza liczba określa stopień n. |
Dane przykładowe (przekopiuj do schowka i
wklej do okna konsoli):
3 25 5 1 64 8 1 144 12 1 |
Pascal// Rozkład LU - macierz odwrotna // Data: 20.10.2010 // (C)2020 mgr Jerzy Wałaszek //----------------------------- program luinv; type NReal = array of double; // typ tablic wierszy MReal = array of NReal; // typ tablicy wskaźników const eps = 1e-12; // Funkcja dokonuje rozkładu LU macierzy A //---------------------------------------- function ludist (n : integer; A : MReal) : boolean; var i, j, k : integer; begin for k := 0 to n - 2 do begin if abs (A [k][k]) < eps then Exit (false); for i := k + 1 to n - 1 do A [i][k] := A [i][k] / A [k][k]; for i := k + 1 to n - 1 do for j := k + 1 to n - 1 do A [i][j] := A [i][j] - A [i][k] * A [k][j]; end; ludist := true; end; // Funkcja wyznacza wektor X na podstawie A i B //--------------------------------------------- function lusolve (k, n : integer; A, X : MReal) : boolean; var i, j : integer; s : double; begin for i := 1 to n - 1 do begin s := 0; for j := 0 to i - 1 do s := s + A [i][j] * X [j][k]; X [i][k] := X [i][k] - s; end; if abs (A [n-1][n-1]) < eps then Exit (false); X [n-1][k] := X [n-1][k] / A [n-1][n-1]; for i := n - 2 downto 0 do begin s := 0; for j := i + 1 to n - 1 do s := s + A [i][j] * X [j][k]; if abs (A [i][i]) < eps then Exit (false); X [i][k] := (X [i][k] - s) / A [i][i]; end; lusolve := true; end; // Program główny var A, X : MReal; n, i, j : integer; ok : boolean; begin // Odczytujemy wymiar macierzy A read (n); // Tworzymy macierze A i X SetLength (A, n); SetLength (X, n); for i := 0 to n - 1 do begin SetLength (A [i], n); SetLength (X [i], n); end; // Odczytujemy dane dla macierzy A for i := 0 to n - 1 do for j := 0 to n - 1 do read (A [i][j]); // Wykonujemy rozkład LU macierzy A if ludist (n, A) then begin // Tworzymy macierz jednostkową w X for i := 0 to n - 1 do begin for j := 0 to n - 1 do X [i][j] := 0; X [i][i] := 1; end; // Wyznaczamy kolejne kolumny macierzy odwrotnej w X ok := true; for i := 0 to n - 1 do if not lusolve (i, n, A, X) then begin ok := false; break; end; end else ok := false; // Jeśli obliczenia się powiodły, wyświetlamy X if ok then begin for i := 0 to n - 1 do begin for j := 0 to n - 1 do write (X [i][j] :10:5, ' '); writeln; end; end else writeln('DZIELNIK ZERO'); // Usuwamy macierze z pamięci for i := 0 to n - 1 do begin SetLength (A [i], 0); SetLength (X [i], 0); end; SetLength (A, 0); SetLength (X, 0); end. |
// Rozkład LU - macierz odwrotna // Data: 20.10.2010 // (C)2020 mgr Jerzy Wałaszek //----------------------------- #include <iostream> #include <iomanip> #include <cmath> using namespace std; const double eps = 1e-12; // Funkcja dokonuje rozkładu LU macierzy A //---------------------------------------- bool ludist (int n, double ** A) { int i, j, k; for(k = 0; k < n - 1; k++) { if(fabs (A [k][k]) < eps) return false; for(i = k + 1; i < n; i++) A [i][k] /= A [k][k]; for(i = k + 1; i < n; i++) for(j = k + 1; j < n; j++) A [i][j] -= A [i][k] * A [k][j]; } return true; } // Funkcja wyznacza wektor X na podstawie A i X [i] //--------------------------------------------------- bool lusolve (int k, int n, double ** A, double ** X) { int i, j; double s; for(i = 1; i < n; i++) { s = 0; for(j = 0; j < i; j++) s += A [i][j] * X [j][k]; X [i][k] -= s; } if(fabs (A [n-1][n-1]) < eps) return false; X [n-1][k] /= A [n-1][n-1]; for(i = n - 2; i >= 0; i--) { s = 0; for(j = i + 1; j < n; j++) s += A [i][j] * X [j][k]; if(fabs (A [i][i]) < eps) return false; X [i][k] = (X [i][k] - s) / A [i][i]; } return true; } // Program główny int main() { double **A, **X; int n, i, j; bool ok; cout << setprecision (5) << fixed; // Odczytujemy wymiar macierzy A cin >> n; // Tworzymy macierze A i X A = new double * [n]; X = new double * [n]; for(i = 0; i < n; i++) { A [i] = new double [n]; X [i] = new double [n]; } // Odczytujemy dane dla macierzy A for(i = 0; i < n; i++) for(j = 0; j < n; j++) cin >> A [i][j]; // Wykonujemy rozkład LU macierzy A if(ludist (n, A)) { // Tworzymy macierz jednostkową w X for(i = 0; i < n; i++) { for(j = 0; j < n; j++) X [i][j] = 0; X [i][i] = 1; } // Wyznaczamy kolejne kolumny macierzy odwrotnej w X ok = true; for(i = 0; i < n; i++) if(! lusolve (i, n, A, X)) { ok = false; break; } } else ok = false; // Jeśli obliczenia się powiodły, wyświetlamy X if(ok) { for(i = 0; i < n; i++) { for(j = 0; j < n; j++) cout << setw (10) << X [i][j] << " "; cout << endl; } } else cout << "DZIELNIK ZERO\n"; // Usuwamy macierze z pamięci for(i = 0; i < n; i++) { delete [] A [i]; delete [] X [i]; } delete [] A; delete [] X; return 0; } |
Basic' Rozkład LU - macierz odwrotna ' Data: 20.10.2010 ' (C)2020 mgr Jerzy Wałaszek '----------------------------- Const eps As Double = 1e-12 ' Funkcja dokonuje rozkładu LU macierzy A '---------------------------------------- Function ludist (n As Integer, A As Double Ptr Ptr) As Integer Dim As Integer i, j, k For k = 0 To n - 2 If Abs (A [k][k]) < eps Then Return 0 For i = k + 1 To n - 1: A [i][k] /= A [k][k] : Next For i = k + 1 To n - 1 For j = k + 1 To n - 1: A [i][j] -= A [i][k] * A [k][j] : Next Next Next ludist = 1 End Function ' Funkcja wyznacza wektor X na podstawie A i X [i] '--------------------------------------------------- Function lusolve (k As Integer, n As Integer, A As Double Ptr Ptr, X As Double Ptr Ptr) As Integer Dim As Integer i, j Dim As Double s For i = 1 To n - 1 s = 0 For j = 0 To i - 1: s += A [i][j] * X [j][k] : Next X [i][k] -= s Next If Abs (A [n - 1][n - 1]) < eps Then Return 0 X [n - 1][k] /= A [n - 1][n - 1] For i = n - 2 To 0 Step -1 s = 0 For j = i + 1 To n - 1: s += A [i][j] * X [j][k] : Next If Abs (A [i][i]) < eps Then Return 0 X [i][k] = (X [i][k] - s) / A [i][i] Next lusolve = 1 End Function Dim As Double Ptr Ptr A, X Dim As Integer i, j, ok, n Open Cons For Input As #1 ' Odczytujemy wymiar macierzy A Input #1, n ' Tworzymy macierze A i X A = New Double Ptr [n] X = New Double Ptr [n] For i = 0 To n - 1 A [i] = New Double [n] X [i] = New Double [n] Next ' Odczytujemy dane dla macierzy A For i = 0 To n - 1 For j = 0 To n - 1 Input #1, A [i][j] Next Next Close #1 ' Wykonujemy rozkład LU macierzy A If ludist (n, A) = 1 Then ' Tworzymy macierz jednostkową w X For i = 0 To n - 1 For j = 0 To n - 1: X [i][j] = 0: Next X [i][i] = 1 Next ' Wyznaczamy kolejne kolumny macierzy odwrotnej w X ok = 1 For i = 0 To n - 1 If lusolve (i, n, A, X) = 0 Then ok = 0 Exit For End If Next Else ok = 0 End If ' Jeśli obliczenia się powiodły, wyświetlamy X If ok = 1 Then For i = 0 To n - 1 For j = 0 To n - 1: Print Using "####.##### ";X [i][j];: Next Print Next Else Print "DZIELNIK ZERO" End If ' Usuwamy macierze dynamiczne For i = 0 To n - 1 Delete [] A [i] Delete [] X [i] Next Delete [] A Delete [] X End |
Wynik: |
3 0.04762 -0.08333 0.03571 -0.95238 1.41667 -0.46429 4.57143 -5.00000 1.42857 |
Zespół Przedmiotowy Chemii-Fizyki-Informatyki w I Liceum Ogólnokształcącym im. Kazimierza Brodzińskiego w Tarnowie ul. Piłsudskiego 4 ©2024 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:
Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.
Informacje dodatkowe.