Równania

// 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ą rozkładu LU
// 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 LUeq(int n, double ** AB, double * X)
{
    int i, j, k, * W, maxw;
    double s, akk, maxe;

    // Tworzymy wektor wierszowy
    W = new int [n];

    // Inicjujemy wektor wierszowy numerami wierszy
    for(i = 0; i < n; i++) W[i] = i;

    // Rozkład LU
    for(k = 0; k < n; k++)
    {
        maxw = k; // numer wiersza z elementem max
        maxe = fabs(AB[W[k]][k]);

        // Szukamy elementu max
        for(i = k + 1; i < n; i++)
            if(fabs(AB[W[i]][k]) > maxe)
            {
                maxw = i;
                maxe = fabs(AB[W[i]][k]);
            }

        // Jeśli max jest zero, to macierz jest zdegenerowana
        if(maxe <= EPS)
        {
            delete [] W;
            return false;
        }

        // Jeśli max nie leży na przekątnej, zamieniamy wiersze
        if(maxw != k) swap(W[k], W[maxw]);

        // Dzielnik
        akk = AB[W[k]][k];

        // Normalizujemy kolumnę
        for(i = k + 1; i < n; i++) AB[W[i]][k] /= akk;

        // Modyfikujemy podmacierz
        for(i = k + 1; i < n; i++)
            for(j = k + 1; j < n; j++) AB[W[i]][j] -= AB[W[i]][k] *  AB[W[k]][j];
    }

    // Teraz rozwiązujemy układ równań w dwóch etapach
    // Najpierw w X wyznaczamy macierz Y stosując podstawienia w przód

    X[0] = AB[W[0]][n];
    for(i = 1; i < n; i++)
    {
        s = AB[W[i]][n];
        for(j = 0; j < i; j++) s -= AB[W[i]][j] * X[j];
        X[i] = s;
    }

    // Wyznaczamy X stosując podstawienia wstecz

    X[n - 1] /= AB[W[n - 1]][n - 1];

    for(i = n - 2; i >= 0; i--)
    {
        s = X[i];
        for(j = n - 1; j > i; j--) s -= AB[W[i]][j] * X[j];
        X[i] = s / AB[W[i]][i];
    }

    delete [] W;
    return true;
}

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

int main()
{
    setlocale(LC_ALL,"");

    cout << setprecision(5) << fixed;

    cout << "Rozwiązywanie układu równań liniowych metodą rozkładu LU" << 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(LUeq(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;
}