Równania

Opisana w poprzednim rozdziale metoda Cramera wyznaczania rozwiązań układu równań liniowych posiada podstawową wadę: dla n równań należy wyliczyć n + 1 wyznaczników stopnia n. Jeśli do obliczenia wyznaczników zastosowany będzie rozkład LU o złożoności obliczeniowej O(n³), to otrzymamy metodę o złożoności O(n⁴). Okazuje się, że układ równań liniowych można rozwiązać ze złożonością O(n³) przy zastosowaniu tzw. eliminacji Gaussa (ang. Gaussian Elimination).

Metoda ta jest dwuetapowa. W pierwszym etapie przekształcamy układ, tak aby wyeliminować z kolejnych równań niewiadome występujące w równaniach poprzednich. W ten sposób w ostatnim równaniu otrzymamy tylko jedną niewiadomą. W drugim etapie równanie to rozwiązujemy w prosty sposób. Następnie wykorzystujemy otrzymane rozwiązania idąc w górę do równania pierwszego do wyliczenia kolejnych niewiadomych.

Przykład:

Mamy następujący układ równań:

Do równania (2) dodajemy stronami równanie (1) pomnożone przez ujemny iloraz współczynników przy niewiadomej x w równaniach (2) i  (1). W tym przypadku iloraz ten jest równy:

Otrzymujemy:

Do równania (3) dodajemy stronami równanie (1) pomnożone przez ujemny iloraz współczynników przy niewiadomej x w równaniach (3) i (1) :

Otrzymujemy nowy układ równań:

Zwróć uwagę, iż w wyniku wykonanych operacji z równań (2') i (3') zniknęła niewiadoma x.

Teraz wyeliminujemy z równania (3') niewiadomą y. W tym celu do równania (3') dodajemy równanie (2') pomnożone przez ujemny iloraz współczynników przy niewiadomej y w równaniach (3') i (2') :

Układ równań wygląda następująco:

Dokonaliśmy eliminacji niewiadomych. Z równania (3") możemy wyliczyć bezpośrednio wartość niewiadomej z:

Wyliczoną wartość niewiadomej z wstawiamy do równania (2') i wyliczamy wartość niewiadomej y:

Mamy już policzone wartości niewiadomych y i z. Wstawiamy te wartości do równania (1) i wyliczamy x:

Otrzymaliśmy pełne rozwiązanie układu równań:

Metoda eliminacji Gaussa w podanej tutaj postaci nie jest niestety niezawodna. Przy odejmowaniu równań mnożymy współczynniki odejmowanego równania przez wartość będącą ilorazem odpowiednich współczynników układu równań. Może się zdarzyć, iż dzielnik ilorazu wynosi 0  (wpada w otoczenie 0 o promieniu ε). Wtedy iloraz nie istnieje i wykonywanie obliczeń należy przerwać. Rozwiązanie tego problemu podajemy w następnym rozdziale.

Algorytm eliminacji Gaussa

Dane wejściowe:

Dane wyjściowe

Zmienne pomocnicze

Lista kroków

K01:	Dla i = 0,1,...,n - 2: wykonuj kroki K02...K05	Dokonujemy eliminacji w równaniach
K02:	Jeśli | AB [ i,i ] | < ε,     to zakończ z wynikiem false	Sprawdzamy, czy dzielnik jest różny od zera
K03:	Dla j = i + 1,..., n - 1:     wykonuj kroki K04...K05	Wyliczamy nowe współczynniki
K04:		Ujemny iloraz
K05:	Dla k = i + 1,..., n:         wykonuj:             AB [ j,k ] ← AB [ j,k ] + m × AB [ i,k ]	Modyfikujemy współczynniki
K06:	Dla i = n - 1,...,0: wykonuj kroki K07...K10	Teraz wyliczamy niewiadome w kierunku odwrotnym
K07:	Jeśli | AB [ i,i ] | < ε,     to zakończ z wynikiem false	Sprawdzamy, czy dzielnik jest różny od zera
K08:	s ← AB [ i,n ]	Do s trafia wyraz wolny
K09:	Dla j = n -1,...,i + 1:     wykonuj: s ← s - AB [ i,j ] × X [ j ]
K10:		Wyliczamy wartość zmiennej
K11:	Zakończ z wynikiem true

Poniższy program jest przykładową realizacją opisanego algorytmu. Dane do programu mają następującą budowę:

Układ równań:

// Eliminacja Gaussa
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 0067
//---------------------------

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

using namespace std;

const double EPS = 0.000000001;

// Funkcja rozwiązuje układ równań
// metodą eliminacji Gaussa
// n - stopień macierzy
// AB - wskaźnik macierzy n*(n+1)
// X - wskaźnik macierzy rozwiązań n*1
// Jeśli układ został rozwiązany,
// to wynikiem jest true.
// Inaczej wynikiem jest false.
//------------------------------------
bool Gauss(int n,
           double ** AB,
           double * X)
{
  int i,j,k;
  double m,s;

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

  // Obliczanie wartości
  // niewiadomych
  for(i = n - 1; i >= 0; i--)
  {
    if(fabs(AB[i][i]) < EPS)
      return false;
    s = AB[i][n];
    for(j = n - 1; j > i; j--)
      s -= AB[i][j] * X[j];
    X[i] = s / AB[i][i];
  }
  return true;
}

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

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

  cout << setprecision(5)
       << fixed;

  cout <<
  "Rozwiązywanie układu równań "
  "liniowych metodą eliminacji "
  "Gaussa\n"
  "----------------------------"
  "----------------------------"
  "------\n\n"
  "Wpisz dane wejściowe:\n\n";

  // Liczba równań w układzie
  int n;

  cin >> n;

  int i,j;

  // Tworzymy dynamicznie
  // potrzebne macierze

  // Macierz współczynników
  // przy niewiadomych x
  // oraz wyrazów wolnych
  double ** AB ;
  // Macierz niewiadomych
  double * X ;

  AB = new double * [n];
  X  = new double [n];

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

  // Odczytujemy dane
  // układu równań
  for(i = 0; i < n; i++)
    for(j = 0; j <= n; j++)
      cin >> AB[i][j];

  cout << endl;

  // Rozwiązujemy układ równań
  if(Gauss(n,AB,X))
  {
    for(i = 0; i < n; i++)
      cout << "x" << i << " = "
           << setw(12) << X[i]
           << endl;
  }
  else
    cout <<
    "Błąd w trakcie obliczeń!!!\n";

  cout << endl;

  // Usuwamy tablice dynamiczne
  for(i = 0; i < n; i++)
    delete [] AB[i];
  delete [] AB;
  delete [] X;

  system("pause");
  return 0;
}

# Eliminacja Gaussa
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0067
#---------------------------

EPS = 0.000000001

# Funkcja rozwiązuje układ równań
# metodą eliminacji Gaussa
# n  - stopień macierzy
# ab - macierz n x (n +1)
# x  - macierz rozwiązań n x 1
# Jeśli układ został rozwiązany,
# to wynikiem jest true.
# Inaczej wynikiem jest false.
#-------------------------------
def gauss(n,ab,x):
    # Eliminacja
    for i in range(n - 1):
        if abs(ab[i][i]) < EPS:
            return False
        for j in range(i + 1, n):
            m = -ab[j][i] / ab[i][i]
            for k in range(i+1, n+1):
                ab[j][k] += m*ab[i][k]

    # Obliczanie wartości
    # niewiadomych
    for i in reversed(range(n)):
        if abs(ab[i][i]) < EPS:
            return False
        s = ab[i][n]
        for j in reversed(range(i+1,n)):
            s -= ab[i][j] * x[j]
        x[i] = s / ab[i][i]
    return True

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

print("Rozwiązywanie układu równań "
      "liniowych metodą eliminacji "
      "Gaussa\n"
      "----------------------------"
      "----------------------------"
      "------\n\n"
      "Wpisz dane wejściowe:\n")

# Liczba równań w układzie
n = int(input())

# Macierz niewiadomych
x = [0.0 for _ in range(n)]

# Odczytujemy dane
# układu równań
ab = [[float(s)
       for s in input().split()]
       for _ in range(n)]
            
print()

# Rozwiązujemy układ równań
if gauss(n,ab,x):
    for i in range(n):
        print(f"x{i} = {x[i]:12.5f}")
else:
    print("Błąd w trakcie obliczeń!!!")

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