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:

ε	–	stała określająca dokładność porównania z zerem.
n	–	liczba równań i niewiadomych w układzie.
AB	–	macierz o rozmiarze n × (n +1), w której znajdują się współczynniki przy niewiadomych oraz wyrazy wolne w kolumnie n + 1.

Dane wyjściowe

true/false	–	określa, czy operacja powiodła się, czy nie.
X	–	macierz n × 1 z wartościami niewiadomych, jeśli wynikiem jest true.

Zmienne pomocnicze

i,j,k	–	indeksy
m	–	wartość ujemnego ilorazu, przez który będą mnożone współczynniki równań przy eliminacjach
s	–	suma

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

#include <iostream>
#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 x (n +1)
// X - wskaźnik macierzy rozwiązań n x 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()
{
    setlocale(LC_ALL,"");

    cout << setprecision(5) << fixed;

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

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

    cin >> n;

    int i,j;

    // Tworzymy dynamicznie potrzebne macierze

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

    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ń!!!" << endl;

    cout << endl;

    // Usuwamy tablice dynamiczne

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

    return 0;
}