Решение СЛАУ методом LU-разложения

LU-разложение — это представление матрицы A в виде A=L•U, где L — нижнетреугольная матрица с еденичной диагональю, а U — верхнетреугольная матрица. LU-разложение является модификациеё метода Гаусса. Основные применения данного алгоритма — решение систем алгебраических уравнений, вычисление определителя, вычисление обратной матрицы и др.

Рассмотрим алгоритм на примере матрицы

Алгоритм

1. Создаем матрицы
   
   и
   
   
2. Для каждого столбца j = 1… 3 матрицы будем вычислять как
   
   
   Для каждой строки вычислим
   
3. Выполняем шаг 2 пока j<=3
4. Получем
   
   и
   
   Такие, что A=L•U
   

Результаты после каждого шага:

1. 
   
   
   
   
2. 
   
   
   
   
3. 
   
   
   
   Проверкой удостоверяемся, что L*U=A
   

Реализация алгоритма на C++

#include <iostream>
#include <vector>

using namespace std;

void LU(vector <vector <double>> A, vector <vector <double>> &L, 
		vector <vector <double>> &U, int n)
{
	U=A;

	for(int i = 0; i < n; i++)
		for(int j = i; j < n; j++)
			L[j][i]=U[j][i]/U[i][i];
	
	for(int k = 1; k < n; k++)
	{
		for(int i = k-1; i < n; i++)
			for(int j = i; j < n; j++)
				L[j][i]=U[j][i]/U[i][i];

		for(int i = k; i < n; i++)
			for(int j = k-1; j < n; j++)
				U[i][j]=U[i][j]-L[i][k-1]*U[k-1][j];
	}

}

void proisv(vector <vector <double>> A, vector <vector <double>> B, 
			vector <vector <double>> &R, int n)
{
	for(int i = 0; i < n; i++)
		for(int j = 0; j < n; j++)
			for(int k = 0; k < n; k++)
				R[i][j] += A[i][k] * B[k][j];
}

void show(vector <vector <double>> A, int n)
{
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < n; j++)
		{
			cout <<"\t"<< A[i][j] << "\t";
		}
		cout << endl;
	}
}

int main()
{
	const int n=4;
	vector <vector <double>> A (n,0), L (n,0), U(n,0), R(n,0);
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < n; j++)
		{
			A[i].push_back(rand()%20-10);
			L[i].push_back(0);
			U[i].push_back(0);
			R[i].push_back(0);
		}
	}
	LU(A,L,U,n);
	cout << "Fisrt matrix" << endl;
	show(A,n);
	cout << "U matrix" << endl;
	show(U,n);
	cout << "L matrix" << endl;
	show(L,n);
	proisv(L,U,R,n);
	cout << "L*U matrix" << endl;
	show(R,n);
	return 0;
}
 

Особенности LU-разложения

1. Легко вычисляется определитель матрицы:
   
   
   
2. Однажды найдя LU-разложение для матрицы мы можем очень быстро решать системы линейных алгебраических уравнений с различной правой частью.
   
   
   Пусть
   Тогда
   Так как — нижнетреугольная матрица, то очень легко находим
   Решаем
   Легко находим , так как — вехнетреугольная матрица
   
   
3. Сложность алгоритма:
   LU-разложение:
   
   Последующее решение систем:
   
   

P.S. Небольшая неровность текста может быть вызвана моей отсутствием большого опыта вставки картинок.