Вы здесь

LU разложение.Нахождение (вычисление) обратной матрицы

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

Теория

В данной статье рассмотрим алгоритм нахождения обратной матрицы при помощи разложения матрицы A на верхнюю треугольную U и на нижнюю треугольную L (алгоритм разложения смотрите в статье «LU разложения матрицы»). Если матрицы A,U,L обратимы, следовательно, справедливо следующее утверждение (рис. 1):


рис. 1. Вычисление обратной матрицы при помощи LU разложения.

Умножив равенство (рис. 1) на U, затем на L (рис. 1), мы получим два справедливых уравнения (рис 2.):


рис. 2. Нахождение обратной матрицы при помощи LU разложения.


Рис. 3. Нахождение обратной матрицы при помощи LU разложения. Развернутое представление матриц.

В первом уравнении (рис. 3) нам не известны элементы X и нижняя часть матрицы L-1. Используя известные элементы матрицы L-1 (0 и 1), мы легко можем выразить часть неизвестных элементов матрицы X. Умножая последнюю строку матрицы U на последний столбец матрицы X, получим 0*X[1][n]+0*X[2][n]+0*X[3][n]+...U[n][n]*X[n][n]=1, отсюда легко получаем X[n][n]=1/U[n][n]. Далее из уравнения 0*X[1][n]+0*X[2][n]+0*X[3][n]+...+U[n-1][n-1]*X[n-1][n]+U[n-1][n]*X[n][n]=0, получим X[n-1][n]=-(U[n-1][n]*X[n][n])/U[n-1][n-1]. Подобным образом мы можем посчитать часть неизвестных элементов матрицы X. Нужно обратить внимание, что поиск неизвестных элементов производится с конца. Из второго уравнения мы легко получим вторую часть неизвестных матрицы X, используя только известные значения матрицы U-1(только 0). Умножая вторую строку матрицы L на первый столбец матрицы X, получим следующее уравнение: L[2][1]*X[1][1]+X[2][1]*1+X[3][1]*0+...+X[n][1]*0=0, из него находим значение X[2][1]=-L[2][1]*X[1][1], далее L[3][1]*X[1][1]+L[3][2]*X[2][1]+1*X[3][1]+0*X[4][1]+0*X[5][1]+...+0*X[n][1]=0, аналогично получим X[3][1]=-(L[3][1] * X[1][1] + L[3][2]*X[2][1]). Подобным образом мы получим недостающую часть неизвестных значений матрицы X. Для расчёта элементов матрицы X можно воспользоваться следующими формулами (рис. 4):


рис. 4. Вычисление значений элементов обратной матрицы при помощи LU разложения.

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

В данном примере я привел код рабочей функции, который раскладывает матрицу A на две - L и U(более подробно о разложении матрицы смотрите в статье «LU разложение матрицы») и находит обратную матрицу A-1.

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