LU - разложение, для решения СЛАУ;

[Marat10.04]

Пусть дано CЛАУ (Сист.Лин.Алгебр.Уравн.)
    Ax = b;
Можем матрицу A[ ][ ] представить в виде произведения двух матриц.
   A = LU
где L - нижн. треугольная матрица,
      U - верхн. треугольная, c единичками на диагонали.

Потом подставив LU вместо A получем:
  LUx = b
Обозначив y = Ux получем две простых систем
  Ly = b  - здесь найдем y
  Ux = y - отсюда x

код на C++ для LU - разложения:

Код: [Выделить]

#include <iostream>
using namespace std;

int main ()
{
  int n;
  double sum = 0;

  cout << "Введите порядок матрицы\n n = ";
  cin >> n;

  double A [ n ][ n ];
  double L [ n ][ n ];
  double U [ n ][ n ];

//задаем матрицу A[][] ...

  for (int i  = 0; i < n; i++)
  {
    for (int j = 0; j < n; j++)
    {
      cout << "\na[" << i << "][" << j << "] = ";
      cin >> A [i][j];

      L [i][j] = 0;
      U [i][j] = 0;

      if (i == j)
        U [i][j] = 1;
    }
  }

//==============================================

//находим первый столбец L[][] и первую строку U[][]

  for (int i = 0; i < n; i++)
  {
    L [i][0] = A [i][0];
    U [0][i] = A [0][i] / L [0][0];
  }

//дальше вычисляем L[][], U[][] по формуле

  for (int i = 1; i < n; i++)
  {
      for (int j = 1; j < n; j++)
      {
          if (i >= j) //нижний треугольник
          {
              sum = 0;
              for (int k = 0; k < j; k++)
                  sum += L [i][k] * U [k][j];

              L [i][j] = A [i][j] - sum;
          }
          else // верхний
          {
              sum = 0;
              for (int k = 0; k < i; k++)
                  sum += L [i][k] * U [k][j];

              U [i][j] = (A [i][j] - sum) / L [i][i];
          }
      }
   }

//====================================================
   cout << "\n\n";
   for (int i = 0; i < n; i++)
   {
     for (int j = 0; j < n; j++)
       cout << "  " << L [i][j] << "   ";
     cout << "\n\n";
   }

   cout << "\n\n";

   for (int i  = 0; i < n; i++)
   {
     for (int j = 0; j < n; j++)
       cout << "  " << U [i][j] << "   ";
     cout << "\n\n";
   }

   return 0;
}