Algorytmy i Struktury Danych - Eliminacja 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ń:

0x₁+2x₂+3x₃+4x₄ = 49
1x₁+0x₂+3x₃+4x₄ = 45
1x₁+2x₂+0x₃+4x₄ = 36
1x₁+2x₂+3x₃+0x₄ = 23

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

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:

x₁ = 2
x₂ = 3
x₃ = 5
x₄ = 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 a_i,j, i = 1…n; j = 1…n oraz w kolumnie n+1 współczynniki b_i układu równań.
X : wektor n elementowy, w którym zostaną umieszczone niewiadome x_i.

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 A; 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: ; do wektora wpisujemy kolejne numery kolumn
     wykonuj: W[i] ← i
K04: Dla i = 1, 2, …, n-1: ; dokonujemy eliminacji współczynników
     wykonuj kroki K05…K11
K05:   k ← i
K06:   Dla j = i+1, i+2, …, n:
       wykonuj: ; wyszukujemy element o największym module
         jeśli |AB[i, W[k]]| < |AB[i, W[j]]|, 
         to k ← j
K07:   W[k] ↔ W[i] ; w wektorze zamieniamy numery
       ; kolumn wg znalezionego elementu

K08:   Jeśli |AB[i, W[i]]| < ε, ; jeśli znaleziony element jest zerem, 
       to zakończ              ; kończymy
K09:   Dla j = i+1, i+2, …, n :
       wykonuj kroki 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: ; sumujemy wiersz j-ty z wierszem i-tym
           AB[j, W[k]] ← AB[j, W[k]]+m×AB[i, W[k]] ; przemnożonym przez m
K12: Dla i = n, n-1, …, 1: ; wyliczamy kolejne niewiadome
     wykonuj kroki K13…K16
K13:   Jeśli |AB[i, W[i]]| < ε, 
       to zakończ
K14:   s ← AB[i, n+1] ; do s trafia współczynnik b
K15:   Dla j = n, n-1, …, i+1:
       wykonuj:
         s ←  s-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

Przykładowe programy

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

2x₂+3x₃+4x₄ = 49
 x₁+3x₃+4x₄ = 45
 x₁+2x₂+4x₄ = 36
 x₁+2x₂+3x₃ = 23

Dane przykładowe (przekopiuj do schowka i wklej do okna konsoli):

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

Pascal

// Eliminacja Gaussa-Crouta
// Data: 3.03.2010
// (C)2020 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.

C++

// Eliminacja Gaussa-Crouta
// Data: 3.03.2010
// (C)2020 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;
}

Basic

' Eliminacja Gaussa-Crouta
' Data: 3.03.2010
' (C)2020 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

Python (dodatek)

# Eliminacja Gaussa-Crouta
# Data: 19.04.2024
# (C)2024 mgr Jerzy Wałaszek
#---------------------------

EPS = 1e-12


# Funkcja realizuje algorytm
# eliminacji Gaussa-Crouta
#---------------------------
def gauss(ab, x):
    n = len(ab)
    w = [i for i in range(n+1)]
     # eliminacja współczynników
    for i in range(n-1):
        k = i
        for j in range(i+1, n):
           if abs(ab[i][w[k]]) < abs(ab[i][w[j]]):
               k = j
        w[k], w[i] = w[i], w[k]
        for j in range(i+1, n):
            if abs(ab[i][w[i]]) < EPS:
                return False
            m = -ab[j][w[i]]/ab[i][w[i]]
            for k in range(i+1, n+1):
                ab[j][w[k]] += m*ab[i][w[k]]
    # wyliczanie niewiadomych
    for i in reversed(range(n)):
        if abs(ab[i][w[i]]) < EPS:
            return False
        s = ab[i][n]
        for j in reversed(range(i+1, n)):
            s -= ab[i][w[j]]*x[w[j]]
        x[w[i]] = s/ab[i][w[i]]
    return True


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

# odczytujemy liczbę niewiadomych
n = int(input())
# tworzymy macierze AB i X
ab = []
x  = [0.0]*n
# odczytujemy dane dla macierzy AB
for i in range(n):
    arr = input().split()
    arr = [float(arr[j]) for j in range(n+1)]
    ab.append(arr)
if gauss(ab, x):
    for i in range(n):
        print("x%d = %9.4f" % (i+1, x[i]))
else:
    print("DZIELNIK ZERO")
print()
# usuwamy macierze z pamięci
ab = None
x = None

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

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.

Eliminacja Gaussa-Crouta