Метод LU – розкладу
(метод розкладу на трикутні матриці або метод LU – факторизації )
Алгоритми цього методу досить близькі до методу Гауса. Основна перевага методу LU – факторизації в порівнянні з методом Гауса є можливість більш простого одержання розв’язку для різних векторів b системи лінійних алгебраїчних рівнянь
A ? x = b (1)
A – матриця розміру n EMBED Equation.3 n із сталими коефіцієнтами;
b – n-мірний вектор вільних членів;
х – n-мірний вектор невідомих.
Матриця системи А представляється у вигляді добутку двох трикутних матриць L та U
А =L ? U, де
EMBED Equation.3 EMBED Equation.3 ,
L – нижня трикутна матриця, U – верхня трикутна матриця.
Зауважимо, що на головній діагоналі матриці U стоять одиниці. Це в свою чергу означає, що визначник матриці А дорівнює добутку діагональних елементів lii матриці L. (det U = 1; det L = l11l22…lnn; det A = det L).
Отже, систему (1) можна записати у вигляді
L ? U ? x = b. (2)
Введемо допоміжний вектор у як
U ? x = у (3)
Тоді можна записати, що
L ? у = b. (4)
Розв’язування системи (1) або (2) відбувається в два етапи: спочатку розв’язується система L ? у = b, а потім U ? x = у.
Така послідовність розв’язку дає перевагу методу (в порівнянні з методом Гауса) – якщо потрібно розв’язати декілька систем з однією і тією ж матрицею А, задача суттєво спрощується, оскільки зберігаються матриці L та U.
Розв’язок системи L ? у = b називається прямим ходом, а системи U ? x = у – оберненим ходом.
Розглянемо прямий хід.Завдяки спеціальній формі матриці L вектор у можна легко визначити. Для цього (4) перепишемо у вигляді системи рівнянь
EMBED Equation.3
Звідси можна одержати, що в загальному вигляді
EMBED Equation.3 (5)
EMBED Equation.3 для EMBED Equation.3
При оберненому ході вектор х визначається з системи рівнянь (3)
x1 + u12x2 + u13x3 + … + u1nxn = y1
x2 + u23x3 + … + u2nxn = y2
x3 + … + u3nxn = y3
………………………………
xn = yn
Починаючи з останнього рівняння, можна послідовно знайти компоненти вектора х. В загальному вигляді вони визначаються за формулами
EMBED Equation.3 для EMBED Equation.3 (6)
Тепер розглянемо LU – розклад матриці А.
Елементи матриці L та U визначаються за наступними формулами:
EMBED Equation.3 де EMBED Equation.3 EMBED Equation.3 , де ( EMBED Equation.3 )
EMBED Equation.3 EMBED Equation.3 , де ( EMBED Equation.3 )
EMBED Equation.3 , де ( EMBED Equation.3 ).
Ці ж самі формули можна записати у більш зручному вигляді (метод Краута).
У варіанті методу, що називається методом Краута, використовується наступна послідовність знаходження елементів матриць L та U ( EMBED Equation.3 ) – де k – крок.
Для всіх EMBED Equation.3
EMBED Equation.3 де EMBED Equation.3 (7)
EMBED Equation.3 де EMBED Equation.3 (8)
Домовимось, як звичайно, вважати значення суми рівним нулю, якщо верхня границя (межа) сумування менша від нижньої.
Крім того, EMBED Equation.3 .
Тобто, при k =1
EMBED Equation.3 для всіх EMBED Equation.3
EMBED Equation.3 для всіх EMBED Equation.3 .
EMBED Equation.3
На місце матриці А записується матриця такого вигляду:
EMBED Equation.3
Алгоритм
1) k =1
Обчислюємо
EMBED Equation.3 для всіх EMBED Equation.3
EMBED Equation.3 для всіх EMBED Equation.3 .
2) k = k + 1
EMBED Equation.3 де EMBED Equation.3
EMBED Equation.3 де EMBED Equation.3 .
Продовжуємо до тих пір, доки k = n.
Зауваження щодо алгоритму побудови програми:
Необхідні в цих співвідношеннях значення елементів матриць L та U обчислюються на попередніх етапах процесу.
Кожний елемент aij матриці А потрібен тільки для обчислення відповідних елементів матриці L та U ( тобто, в подальшому процесі aij не потрібні). Оскільки нульові елементи матриць L та U, а також одиничну діагональ матриці U запам’ятовувати не потрібно, тобто в процесі обчислення матриці L та U можуть бути записані на місце матриці А. Причому матриця L розташована в нижньому трикутнику (i ? j), U – відповідно в верхньому трикутнику (i < j) матриці А.
Розклад матриці А на матриці L та U , як правило, об’єднують з прямим ходом. Може бути використана наступна послідовність обчислень: спочатку обчислюється перший стовпець матриці L , потім перший рядок матриці U та перший елемент вектора у; далі – другий стовпець матриці L , другий рядок матриці U та другий елемент вектора у і т.д.
Матрицю А можна розкласти на дві трикутні матриці L та U при умові, що головні діагональні мінори матриці А відмінні від нуля. Якщо ця умова не виконується, наприклад, для k-го головного діагонального мінора, то при обчисленні елемента lkk матриці L за формулою (7) він стане рівним нулю. Це в свою чергу призведе до неможливості обчислення елементів ukj матриці U (ділення на нуль). Усунути таку ситуацію можна було б шляхом перевірки умов нерівності нуля головних діагональних мінорів до розкладу А = L U . Однак така перевірка пов’язана із значними затратами машинного часу, що перевищує затрати часу пов’язані із розкладом А = L U (крім того перевірка – це обчислення тих самих визначників).
Тому доцільніше проводити перевірку умови lkk = 0 безпосередньо в процесі обчислення елементів lij та uij. Тоді при виконанні цієї умови потрібно переставляти відповідний k-й рядок матриці А з наступними k + s рядками ( EMBED Equation.3 ), до тих пір, доки не виконається умова EMBED Equation.3