Вы здесь

LU разложение. Решение систем линейных уравнений

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

Теория

В данной статье описан и реализован алгоритм решения систем линейных уравнений при помощи LU разложения матрицы (алгоритм разложения смотрите в статье «LU разложение матрицы»). Систему уравнений можно представить в матричном виде — L*U*x=b, где L*U=A. Для удобства расчётов введём вектор дополнительных переменных y=(y1,y2,y3,..yn) так, что U*x=y. Уравнение L*U*x=b можно будет переписать в виде системы уравнений (рис 1):


рис. 1. Решение систем уравнений при помощи LU разложения матрицы.

Представим уравнение L*y=b в развёрнутом виде (рис. 2):


рис.2. Решение систем уравнений при помощи LU разложения матрицы. Представление уравнения L*y=b в развёрнутом виде.

Из уравнения L*y=b последовательно получаем значения y. Умножаем первую строку матрицы L на вектор y и получим следующее равенство 1*y1+0*y2+0*y3+...+0*yn=b1 значение y1=b1. Вторую строку матрицы L умножаем на вектор y, получим следующее уравнение L21*y1+1*y2+0*y3+0*y4+...+0*yn=b2, отсюда выразим y2=b2-L21*y1. Далее yn=bn-(Ln1*y1+Ln2*y2...). Зная значения y, аналогично посчитаем значения вектора x в обратном порядке (рис. 3):


рис.3 Решение систем уравнений при помощи LU разложения матрицы. Представление уравнения U*x=y в развёрнутом виде.

Умножаем последнюю строку матрицы U на вектор x, получим следующее уравнение 0*x1+0*x2+0*x3+...+Unn*xn=yn, отсюда xn=yn/Unn. Вторую строку матрицы U умножим на вектор x, равенство приобретает вид 0*x1+U22*x2+U23*x3+...+U2n*xn=y2, из этого равенства выразим x2=y2-(U23*x3+.....+U2n*xn)/U22. Далее x1=y1-(U12*x2+U13*x3+.....+U1n*xn)/U11. В конце теоретической части приведу формулы расчёта x и y (рис.4):


рис.4. Решение систем уравнений при помощи LU разложения матрицы. Формулы расчёта вектора x и y.

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

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

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