Вы здесь

LU разложение матрицы

Исходный код на C++: 

Теория

Пусть дана матрица A размером n*n, если матрица A не вырождена, то существует разложение матрицы A на верхнюю треугольную U и на нижнюю треугольную L. Если все элементы диагонали одной из матриц U или L фиксированы и не равны нулю, то такое разложение единственно. Найдём такие элементы матрицы L и матрицы U (рис. 1):

lu факторизация
рис. 1. LU разложение матрицы.

Умножим поэлементно матрицу L на U и получим:
Элемент матрицы A11=1*U11+0*0+0*0+...0*0;
Элемент матрицы A12=1*U12+0*U22+0*0+...0*0;
Элемент матрицы A13=1*U13+0*U23+0*U33+...0*0;
Элемент матрицы A1n=1*U1n+0*U2n+0*U3n+...0*Unn;
U11=A11, U12=A12, U13=U13,...U1n=A1n;
Элемент матрицы A21=L21*U11+1*0+0*0+...0*0;
Элемент матрицы A22=L21*U12+1*U22+0*0+...0*0;
Элемент матрицы A23=L21*U13+1*U23+0*U33+...0*0;
Элемент матрицы A2n=L21*U1n+1*U2n+0*U3n+...0*Unn;
L21*U11=A21, L21*U12+U22=A22, L21*U13+U23=A23,...L21*U1n+U2n=A2n;
…....
Элемент матрицы An1=Ln1*U11+Ln2*0+0*0+...0*0;
Элемент матрицы An2=Ln1*U12+Ln2*U22+Ln3*0+...0*0;
Элемент матрицы An3=Ln1*U13+Ln2*U23+Ln3*U33+...0*0;
Элемент матрицы Ann=Ln1*U1n+Ln2*U2n+Ln3*U3n+...1*Unn;
Ln1*U11=An1, Ln1*U12+Ln2*U22=An2, Ln1*U13+Ln2*U23+Ln3*U33=An3,...Ln1*U1n +Ln2*U2n+ Ln3*U3n+...1*Unn = Ann.
На основании этого легко получить элементы матриц L и U. Первая строка матрицы U равна первой строке матрицы A (U11=A11,U12=A12,...U1n=A1n). Используя уравнение L21*U11=A21, мы легко найдем значение L21=A21/U11. Зная L21 и U11, U12...U1n, найдём вторую строку матрицы U (U22=A22-L21*U12, U23=A23-L21*U13,....U2n=A2n-L21*U1n и т.д). Используя формулы, можно посчитать элементы матриц L (рис. 2) и U (рис. 3).

Вычисление  элементов нижней треугольной матрицы L
рис. 2. LU разложение матрицы. Формула расчёта элементов нижней треугольной матрицы L при i>j

Вычисление  элементов верхней треугольной матрицы U
рис. 3. LU разложение матрицы. Формула расчёта элементов верхней треугольной матрицы U при i<=j

В данной реализации матрицы L и U хранятся в виде матрицы (рис. 4), что позволяет экономить память.

Визуализация матрицы LU
Рис. 4. Представление LU матрицы.

Программная реализация

int LU_fun(int cnt_str,double **mass,double **&amp;LU)
{
int i,j,k;
LU=new double *[cnt_str];//создаём массив под матрицу LU
for(i=0;i&lt;cnt_str;i++)
LU[i]=new double [cnt_str];
for(i=0;i&lt;cnt_str;i++)
{
for(j=0;j&lt;cnt_str;j++)
{
double sum=0;
if(i&lt;=j)
{
for(k=0;k&lt;i;k++)
sum+=LU[i][k]*LU[k][j];
LU[i][j]=mass[i][j]-sum;//вычисляем элементы верхней треугольной матрицы
}
else
{
for(k=0;k&lt;j;k++)
sum+=LU[i][k]*LU[k][j];
if(LU[j][j]==0)
return 0;
LU[i][j]=(mass[i][j]-sum)/LU[j][j];//вычисляем элементы нижней треугольной матрицы
}
}
}
return 1;
}