Równania

// Rozkład LU w równaniach
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 0069
//---------------------------

#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ą 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++)
  {
    // numer wiersza
    // z elementem max
    maxw = k;
    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()
{
  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(5)
       << fixed;

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

  cout << endl;

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

  system("pause");
  return 0;
}