Eliminacja Gaussa-Crouta


Tematy pokrewne
Macierze
Podstawowe pojęcia dotyczące macierzy
Reprezentacja macierzy w pamięci komputera
Mnożenie macierzy przez skalar
Dodawanie macierzy
Transponowanie macierzy
Mnożenie macierzy
Potęgowanie macierzy
Eliminacja Gaussa
Eliminacja Gaussa-Crouta
Wyznacznik macierzy
Rozkład LU
Rozkład LU – rozwiązywanie układu równań liniowych
Rozkład LU – macierz odwrotna
Rozkład LU – wyznacznik macierzy

Problem

Znaleźć rozwiązanie układu równań liniowych postaci:

 

a1,1x1 + a1,2x2 + ... + a1,nxn = b1
a2,1x1 + a2,2x2 + ... + a2,nxn = b2
...
an,1x1 + an,2x2 + ... + an,nxn = bn

 

Wykorzystując metodę eliminacji Gaussa-Crouta.

 

 

Opisany w poprzednim rozdziale algorytm eliminacji Gaussa posiada dosyć istotną wadę. Bazuje on na dzieleniu przez elementy przekątnej głównej macierzy współczynników. Otóż może się zdarzyć, iż element przekątnej głównej jest równy zero lub zostanie wyzerowany w wyniku obliczeń. W takim przypadku algorytm zostanie przerwany z błędem dzielenia przez zero. Jako przykład niech posłuży poniższy układ równań:

 

0x1 + 2x2 + 3x3 + 4x4 = 49
1x1 + 0x2 + 3x3 + 4x4 = 45
1x1 + 2x2 + 0x3 + 4x4 = 36
1x1 + 2x2 + 3x3 + 0x4 = 23

 

Dane dla programów dla tego układu są następujące:

 

4
0 2 3 4 49
1 0 3 4 45
1 2 0 4 36
1 2 3 0 23

 

Po wprowadzeniu tych danych do przykładowego programu umieszczonego w poprzednim rozdziale otrzymujemy komunikat "DZIELNIK ZERO". Tymczasem rozwiązanie równania istnieje i jest następujące:

 

x1 = 2
x2 = 3
x3 = 5
x4 = 7

 

Ulepszeniem algorytmu eliminacji Gaussa jest metoda Crouta, która polega na tym, iż na początku eliminacji wyszukujemy w wierszu macierzy AB element o największym module, po czym zamieniamy miejscami kolumnę ze znalezionym elementem z kolumną zawierającą element głównej przekątnej. W ten sposób dzielnik będzie posiadał największą na moduł wartość i pozbędziemy się sytuacji, gdy może on posiadać wartość 0 (chyba że cały wiersz macierzy AB jest zerowy – ale wtedy równania i tak się nie da rozwiązać).

Tak zmodyfikowany algorytm eliminacji Gaussa wymaga zapamiętania, które kolumny zostały zamienione, ponieważ na etapie wyznaczania niewiadomych musimy znać numery wyliczanych niewiadomych. Osiągniemy to wprowadzając do algorytmu dodatkowy wektor W przechowujący informację o numerach kolumn – zamiana miejscami kolumn nie wymaga faktycznego ich przemieszczania w pamięci – wystarczy odwoływać się do nich poprzez wektor kolumn. Jeśli w wektorze W zamienimy ze sobą numery dwóch kolumn (są to tylko dwa elementy, zatem wymiana jest w czasie stałym), to efekt będzie taki, jakby te kolumny zostały zamienione miejscami w macierzy AB:

 

 

Na powyższym rysunku widzimy odwołanie do kolumny 2 macierzy poprzez wektor numerów kolumn, oznaczony tutaj kolorem żółtym. Gdy kolumny nie były zamieniane, to odwołanie prowadzi do tej samej kolumny – na pozycji 2 w wektorze kolumn jest numer 2. Jednakże po zamianie (tutaj kolumny nr 2 z nr 5) odwołanie do kolumny nr 2 przenosi nas do kolumny nr 5, chociaż w macierzy docelowej kolumny te znajdują się dalej na swoich miejscach. Rozwiązanie to zaoszczędzi nam czasu przy przenoszeniu elementów. Dodatkowo zawartość wektora numerów kolumn może być wykorzystana w końcowej fazie algorytmu przy wyznaczaniu wartości zmiennych.

 

Algorytm eliminacji Gaussa-Crouta

Wejście
n    liczba niewiadomych
AB  – macierz n × (n + 1) zawierająca współczynniki ai,j, i = 1...n; j = 1...n oraz w kolumnie n + 1 współczynniki bi układu równań
X  – wektor n elementowy, w którym zostaną umieszczone niewiadome xi.
Wyjście:

Jeśli zmienna r ma wartość true, to wektor X zawiera rozwiązanie układu równań. W przeciwnym razie nastąpiło dzielenie przez zero i algorytm nie może znaleźć poprawnych rozwiązań układu równań. Pierwotna macierz współczynników AB zostaje zniszczona – zastąpiona macierzą trójkątną, powstałą w czasie redukcji współczynników.

Elementy pomocnicze:
W  – wektor n + 1 elementowy zawierający numery kolumn macierzy AB. W N
i,j,k    zmienne sterujące pętli. i,j,k N
m   mnożnik, przez który będą mnożone elementy macierzy AB. m R
s  – zlicza sumę iloczynów. s R
ε  – określa dokładność porównania z zerem; ε = 10-12, ε R.
Lista kroków:
K01: r ← false ; zakładamy porażkę
K02: ε ← 10-12 ; określamy dokładność porównania z zerem
K03: Dla i = 1,2,...,n+1: wykonuj W[i] ← i ; do wektora wpisujemy kolejne numery kolumn
K04: Dla i = 1,2,...,n-1 wykonuj K05...K11 ; dokonujemy eliminacji współczynników
K05:     ki  
K06:     Dla j = i+1, i+2, ...,n wykonuj:
        jeśli | AB[i,W[k]] | < | AB[i,W[j]] |, to kj
; wyszukujemy element o największym module
K07:     W[k] ↔ W[i] ; w wektorze zamieniamy numery kolumn wg
; znalezionego elementu
K08:     Jeśli | AB[i,W[i]] | < ε, to zakończ ; jeśli znaleziony element jest zerem, kończymy
K09:     Dla j = i+1, i+2, ..., n wykonuj K10...K11  
K10:
        m ← -    AB[j,W[i]]
AB[i,W[i]]
; wyznaczamy mnożnik
K11:         Dla k = i+1, i+2, ..., n+1 wykonuj:
            AB[j,W[k]] ← AB[j,W[k]] + m × AB[i,W[k]]
; sumujemy wiersz j-ty z wierszem i-tym
; przemnożonym przez m
K12: Dla i = n, n-1, ..., 1 wykonuj K13...K16 ; wyliczamy kolejne niewiadome
K13:     Jeśli | AB[i,W[i]] | < ε, to zakończ  
K14:     sAB[i,n+1] ; do s trafia współczynnik b
K15:     Dla j = n, n-1, ..., i+1 wykonuj: ss - AB[i,W[j]] × X[W[j]] ; zliczamy sumę iloczynów współczynników przez niewiadome
K16:
    X[W[i]] ← -   s
AB[i,W[i]]
; obliczamy kolejną niewiadomą
K17: r ← true ; zaznaczamy sukces
K18: Zakończ  

 

Program

Ważne:

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 współczynników AB i wylicza niewiadome X. Jeśli obliczenie nie może zostać wykonane, program wypisuje komunikat "DZIELNIK ZERO".

Dane dla macierzy AB są zdefiniowane następująco:

Pierwsza liczba określa liczbę niewiadomych n.

Następne n × (n + 1) liczb określa zawartość macierzy AB wierszami – pierwsze n liczb każdego wiersza odnosi się do współczynników a, a (n + 1)-sza liczba odnosi się do wyrazu wolnego b. Poniżej mamy przykładowe dane wejściowe dla programu:

Równanie:

 

2x2 + 3x3 + 4x4 = 49
x1 + 3x3 + 4x4 = 45
x1 + 2x2 + 4x4 = 36
x1 + 2x2 + 3x3 = 23

 

Dane dla programu:

 

4
0 2 3 4 49
1 0 3 4 45
1 2 0 4 36
1 2 3 0 23

 

Lazarus
// Eliminacja Gaussa-Crouta
// Data: 3.03.2010
// (C)2012 mgr Jerzy Wałaszek
//-----------------------------

program gausscrout;

type
  NInt  = array of integer; // typ wektora kolumn
  NReal = array of double;  // typ tablic wierszy
  MReal = array of NReal;   // typ tablicy wskaźników

const
  eps = 1e-12;

// Funkcja realizuje algorytm eliminacji Gaussa-Crouta
//----------------------------------------------------
function Gauss(n : integer; AB : MReal; X : NReal; W : NInt) : boolean;
var
  i,j,k : integer;
  m,s   : double;
begin
  
  for i := 0 to n do W[i] := i;

  // eliminacja współczynników

  for i := 0 to n - 2 do
  begin
    k := i;
    for j := i + 1 to n - 1 do
      if abs(AB[i][W[k]]) < abs(AB[i][W[j]]) then k := j;
    j := W[k]; W[k] := W[i]; W[i] := j;
    for j := i + 1 to n - 1 do
    begin
      if abs(AB[i][W[i]]) < eps then exit(false);
      m := -AB[j][W[i]]/ AB[i][W[i]];
      for k := i + 1 to n do
        AB[j][W[k]] := AB[j][W[k]] + m * AB[i][W[k]];
    end;
  end;

  // wyliczanie niewiadomych

  for i := n - 1 downto 0 do
  begin
    if abs(AB[i][W[i]]) < eps then exit(false);
    s := AB[i][n];
    for j := n - 1 downto i + 1 do
      s := s - AB[i][W[j]] * X[W[j]];
    X[W[i]] := s / AB[i][W[i]];
  end;
  Gauss := true;
end;

// Program główny

var
  AB    : MReal;
  X     : NReal;
  W     : NInt;
  n,i,j : integer;

begin

  // odczytujemy liczbę niewiadomych

  read(n);

  // tworzymy macierze AB, X i W o odpowiednich rozmiarach

  SetLength(AB,n);
  SetLength(X,n);
  SetLength(W,n+1);

  for i := 0 to n - 1 do SetLength(AB[i],n + 1);

  // odczytujemy dane dla macierzy AB

  for i := 0 to n - 1 do
    for j := 0 to n  do read(AB[i][j]);

  if Gauss(n,AB,X,W) then
  begin
    for i := 0 to n - 1 do writeln('x',i+1,' = ',X[i]:9:4);
  end
  else
    writeln('DZIELNIK ZERO');

  // usuwamy macierze z pamięci

  for i := 0 to n - 1 do SetLength(AB[i],0);
  SetLength(AB,0);
  SetLength(X,0);
  SetLength(W,0);

end.
Code::Blocks
// Eliminacja Gaussa-Crouta
// Data: 3.03.2010
// (C)2012 mgr Jerzy Wałaszek
//-----------------------------

#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

const double eps = 1e-12;

// Funkcja realizuje algorytm eliminacji Gaussa-Crouta
//----------------------------------------------------
bool gauss(int n, double ** AB, double * X, int * W)
{
  int    i,j,k;
  double m,s;

  for(i = 0; i <= n; i++) W[i] = i;
  
  // eliminacja współczynników

  for(i = 0; i < n - 1; i++)
  {
    k = i;
    for(j = i + 1; j < n; j++)
      if(fabs(AB[i][W[k]]) < fabs(AB[i][W[j]])) k = j;
    swap(W[k],W[i]);
    for(j = i + 1; j < n; j++)
    {
      if(fabs(AB[i][W[i]]) < eps) return false;
      m = -AB[j][W[i]] / AB[i][W[i]];
      for(k = i + 1; k <= n; k++)
        AB[j][W[k]] += m * AB[i][W[k]];
    }
  }

  // wyliczanie niewiadomych

  for(i = n - 1; i >= 0; i--)
  {
    if(fabs(AB[i][W[i]]) < eps) return false;
    s = AB[i][n];
    for(j = n - 1; j >= i + 1; j--)
      s -= AB[i][W[j]] * X[W[j]];
    X[W[i]] = s / AB[i][W[i]];
  }
  return true;
}

// Program główny

int main()
{
  double **AB, *X;
  int n,i,j,*W;

  cout << setprecision(4) << fixed;
  
  // odczytujemy liczbę niewiadomych

  cin >> n;

  // tworzymy macierze AB, X i W

  AB = new double * [n];
  X  = new double [n];
  W  = new int [n + 1];

  for(i = 0; i < n; i++) AB[i] = new double[n + 1];

  // odczytujemy dane dla macierzy AB

  for(i = 0; i < n; i++)
    for(j = 0; j <= n; j++) cin >> AB[i][j];

  if(gauss(n,AB,X,W))
  {
    for(i = 0; i < n; i++)
      cout << "x" << i + 1 << " = " << setw(9) << X[i]
           << endl;
  }
  else
    cout << "DZIELNIK ZERO\n";

  // usuwamy macierze z pamięci

  for(i = 0; i < n; i++) delete [] AB[i];
  delete [] AB;
  delete [] X;
  delete [] W;

  return 0;
}
Free Basic
' Eliminacja Gaussa-Crouta
' Data: 3.03.2010
' (C)2012 mgr Jerzy Wałaszek
'-----------------------------

Const eps As Double = 1e-12

' Funkcja realizuje algorytm eliminacji Gaussa-Crouta
'----------------------------------------------------
Function Gauss(n As Integer, AB As Double Ptr Ptr, X As Double Ptr, W As Integer Ptr) As Integer

  Dim As Integer i,j,k
  Dim As Double m,s
  
  For i = 0 To n
    W[i] = i
  Next

  ' eliminacja współczynników

  For i = 0 To n - 2
    k = i
    For j = i + 1 To n - 1
      If Abs(AB[i][W[k]]) < Abs(AB[i][W[j]]) Then k = j
    Next
    Swap W[k], W[i] 
    For j = i + 1 To n - 1
      If Abs(AB[i][W[i]]) < eps Then Return 0
      m = -AB[j][W[i]]/ AB[i][W[i]]
      For k = i + 1 To n
        AB[j][W[k]] += m * AB[i][W[k]]
      Next
    Next
  Next
  
  ' wyliczanie niewiadomych

  For i = n - 1 To 0 Step -1
    If Abs(AB[i][W[i]]) < eps Then Return 0
    s = AB[i][n]
    For j = n - 1 To i + 1 Step -1
      s -= AB[i][W[j]] * X[W[j]]
    Next
    X[W[i]] = s / AB[i][W[i]]
  Next
  Gauss = 1
End Function

' Program główny
'---------------

Dim AB As Double Ptr Ptr
Dim X As Double Ptr
Dim W As Integer Ptr
Dim As Integer n,i,j

Open Cons For Input As #1

' odczytujemy liczbę niewiadomych

Input #1,n

' tworzymy macierze AB, X i W o odpowiednich rozmiarach

AB = New Double Ptr [n]
X  = New Double [n]
W  = New Integer [n + 1]

For i = 0 To n - 1
  AB[i] = New Double [n + 1]
Next

' odczytujemy dane dla macierzy AB

For i = 0 To n - 1
  For j = 0 To n
    Input #1,AB[i][j]
  Next
Next

Close #1

If Gauss(n,AB,X,W) = 1 Then
  For i = 0 To n - 1
    Print Using "x# = ####.####"; i + 1; X[i]
  Next
Else
  Print "DZIELNIK ZERO"
End If

' usuwamy macierze z pamięci

For i = 0 To n - 1
  Delete [] AB[i]
Next
Delete [] AB
Delete [] X
Delete [] W

End
Wynik
4
0 2 3 4 49
1 0 3 4 45
1 2 0 4 36
1 2 3 0 23
x1 =    2.0000
x2 =    3.0000
x3 =    5.0000
x4 =    7.0000

 



List do administratora Serwisu Edukacyjnego Nauczycieli I LO

Twój email: (jeśli chcesz otrzymać odpowiedź)
Temat:
Uwaga: ← tutaj wpisz wyraz  ilo , inaczej list zostanie zignorowany

Poniżej wpisz swoje uwagi lub pytania dotyczące tego rozdziału (max. 2048 znaków).

Liczba znaków do wykorzystania: 2048

 

W związku z dużą liczbą listów do naszego serwisu edukacyjnego nie będziemy udzielać odpowiedzi na prośby rozwiązywania zadań, pisania programów zaliczeniowych, przesyłania materiałów czy też tłumaczenia zagadnień szeroko opisywanych w podręcznikach.



   I Liceum Ogólnokształcące   
im. Kazimierza Brodzińskiego
w Tarnowie

©2017 mgr Jerzy Wałaszek

Dokument ten rozpowszechniany jest zgodnie z zasadami licencji
GNU Free Documentation License.