6.7.Решение систем линейных уравнений с использованием lu-разложения матрицы системы

LU-разложением (LU-факторизацией) называют представление матрицы A в виде произведения LU, где L — нижняя треугольная матрица с единичной диагональю, а U — верхняя треугольная матрица:

_(6.15)

Метод решения систем линейных уравнений, основанный на данном разложении, является разновидностью метода Гаусса.

Обозначим элементы L через l_ij, элементы U - через u_ij, (i,j=1,...,n). По свойствам 9 и 10 определителей (Глава 5):

det(L)=1; det(L) = l₁₁  l₂₂ …l_nn = 1;

det(A) = det(L)det(U) = 1det(U) = u₁₁  u₂₂ …u_nn.

6.7.1. Идея использования lu-разложения и алгоритм прямого метода расчета

Идея LU-разложения основана на матричном представлении прямого хода метода Гаусса без выбора главных элементов строк. Обозначив для общности матрицу начальной системы уравнений А через А⁽⁰⁾, представим исходную систему уравнений в виде:

А⁽⁰⁾X =B. (6.16)

При выполнении прямого хода вначале обнуляются элементы столбца 1 матрицы системы под главной диагональю путем сложения строк расширенной матрицы с номерами 2,3,…,n с первой строкой, умноженной, соответственно, на коэффициенты (-а⁽⁰⁾₂₁/а⁽⁰⁾₁₁), (-а⁽⁰⁾₃₁/а⁽⁰⁾₁₁),…, (-а⁽⁰⁾_n1/а⁽⁰⁾₁₁). Это преобразование эквивалентно умножению левой и правой частей уравнения (6.16) на нижнюю треугольную матрицу с единичной диагональю вида:

С₁ (А⁽⁰⁾)А⁽⁰⁾X = С₁(А⁽⁰⁾)B.

Обозначим произведение С₁(А⁽⁰⁾)А⁽⁰⁾ через А⁽¹⁾ и представим преобразованную систему в виде А⁽¹⁾X = С₁(А⁽⁰⁾)B. Для того, чтобы в данной системе обнулить элементы столбца 2 матрицы системы А⁽¹⁾под главной диагональю, складываем строки расширенной матрицы с номерами 3,4,…,n со второй строкой, умноженной, соответственно, на коэффициенты (-а⁽¹⁾₃₂/а⁽¹⁾₂₂), (-а⁽¹⁾₄₂/а⁽¹⁾₂₂),…, (-а⁽¹⁾_n2/а⁽¹⁾₂₂). Это преобразование эквивалентно умножению левой и правой частей преобразованного уравнения на нижнюю треугольную матрицу с единичной диагональю вида:

С₂(А⁽¹⁾)А⁽¹⁾X = С₂(А⁽¹⁾)С₁(А⁽⁰⁾)B.

Затем в полученной новой системе уравнений умножаем обе части на матрицы С₃(А⁽²⁾), С₄(А⁽³⁾),… С_n-1(А^(n-2)). В итоге всех умножений, вводя обозначение С = С _n-1(А^(n-2)) С _n-2(А^(n-3))…С₂(А⁽¹⁾)С₁(А⁽⁰⁾), получим систему вида:

С А⁽⁰⁾X = С B, (6.17)

у которой матрица СА⁽⁰⁾ после завершения прямого хода имеет верхний треугольный вид. Матрица С является произведением (n-1) нижних треугольных матриц с единичной диагональю. Несложно показать, что С также будет нижней треугольной с единичной диагональю. Так как det(C) = 1, то обратная матрица C^-1 существует. Нетрудно доказать, что она также является нижней треугольной с единичной диагональю. Умножая обе части системы (6.17) на C^-1, получим исходную систему уравнений (6.16) в следующей форме:

С^-1  С А⁽⁰⁾X =B,

Матрица С^-1 удовлетворяет требованиям к матрице L, а произведение С А⁽⁰⁾ - к U, т.е.: L=С^-1, U=С А⁽⁰⁾ – матрица исходной системы после прямого хода.

Для определения коэффициентов матриц L и U используют прямой метод расчета, основанный на правиле умножения матриц и учете структуры матриц L и U. По нему на каждом шаге k (k =1,2,…,n) по ранее найденным элементам матриц определяются строка матрицы U и столбец матрицы L с номером k. Для примера рассмотрим вывод формул для k =1,2.

k =1. Из структуры матриц L и U следует, что элементы а_1j первой строки матрицы А образуются умножением на 1 (l₁₁) элементов u_1j: а_1j = l₁₁ u_1j = u_1j. Поэтому в первой строке матрицы U: u_1j = а_1j, (j = 1,...,n). Элементы а_i1 первого столбца матрицы А образуются умножением l_i1 на u₁₁: а_i1 = l_i1u₁₁. Поэтому в первом столбце матрицы L: l_i1 = а_i1/u₁₁, (i = 2,...,n). В итоге расчетные формулы шага i =1 с учетом l₁₁ =1 имеют вид:

u_1j = а_1j, (j = 1,...,n); l₁₁ =1; l_i1 = а_i1/ u₁₁, (i = 2,...,n). (6.18.1)

k =2. Элементы а_2j второй строки матрицы А рассчитываются по формуле: а_2j= l₂₁ u_1j+ l₂₂  u_2j. Учитывая, что l₂₂=1, u₂₁=0 и первая строка U уже найдена, формула для расчета элементов ее строки 2 имеет вид: u_2j = а_2j - l₂₁ u_1j, (j = 2,...,n). Аналогично из формулы для элементов второго столбца матрицы А: а_i2 = l_i1 u₁₂+ l_i2u₂₂, учитывая, что l₁₂=0, l₂₂=1 и первый столбец L уже найден, расчетная формула элементов ее столбца 2 имеет вид: l_i2 = (а_i2 - l_i1 u₁₂)/(u₂₂), (i = 3,...,n).. В итоге расчетные формулы шага i =2 имеют вид:

u₂₁=0; u_2j = а_2j - l₂₁ u_1j, (j = 2,...,n);

l₁₂=0; l₂₂=1; l_i2 = (а_i2 - l_i1 u₁₂)/(u₂₂). (i = 3,...,n). (6.18.2)

По аналогии можно вывести расчетные формулы для произвольного шага k (k = 1,...,n), который заключается в расчете строки матрицы U и столбца L с номерами k:

u_k1=u_k2=…=u_k(k-1)=0;

u_kj = а_kj – ( l_k1 u_1j + l_k2 u_2j +…+ l_k(k-1) u_(k-1)j), (j = k,...,n); (6.18.k.1)

l_1k=l_2k=…=l_(k-1)k =0; l_kk = 1;

l_ik = [а_jk – ( l_j1 u_1k + l_j2 u_2k +…+ l_j(k-1) u_(k-1)k)]/u_kk, (i = (k +1),...,n). (6.18.k.2)

Вычисления по формулам (6.18.k) начинаются с расчета величины диагонального элемента u_kk матрицы U.

Если u_kk  0, то при (k = 1,...,n-1), можно продолжать дальнейшие расчеты по формулам (6.18.k.1)-(6.18.k.2), поскольку все величины в них определены.

Если u_kk = 0 при (k = 1,...,n-1), то из-за деления на нуль возникает неопределенность в расчете коэффициентов l_ik. Следовательно, для заданной исходной матрицы А LU-разложение не существует и следует прекратить выполнение алгоритма.