Алгоритм LU-разложения

Данный алгоритм можно рассматривать как конкретную форму метода Гаусса. Алгоритм LU-разложения используется не только для решения СЛАУ, но и также для обращения матрицы, т.е. вычисления матрицы, обратной данной.

Пусть и будем искать представление A в виде:

,

(3.10)

где L и U – соответственно нижняя и верхняя треугольные матрицы вида

.

Известно, что если все угловые миноры матрицы А отличны от нуля, т.е.

то разложение вида (3.10) существует и единственно. Для того чтобы получить расчётные формулы, поступим следующим образом. Обозначим a_ij произведение i-ой строки матрицы L на j-ый столбец матрицы R, причём будем считать вначале, что i<j.

Тогда .

Выразим из последней формулы u_ij.

.

(3.11)

Как это принято, будем считать в формуле (3.11) и далее, что сумма вида равна нулю, если значение верхней границы индекса суммирования меньше нижней границы.

В случае i = j имеем

Учитывая, что u_ii º 1 и, выражая из последнего соотношения l_ii, получаем:

.

(3.12)

Наконец, при i > j получаем

откуда, с учетом того, что u_jj º 1 приходим к формуле

.

(3.13)

Итак, расчетные формулы (3.11) – (3.13) получены. Для того чтобы при их применении не использовались неизвестные (не вычисленные) величины, необходимо выбрать соответствующий порядок вычисления элементов матриц L и U.

Например, можно рекомендовать порядок расчета элементов матриц L и U, схематически изображенный на рис. 3.1. На нем цифры слева для матрицы L и сверху - для матрицы U означают, что на первом шаге рассчитывается l₁₁ по формуле (3.12), затем вычисляется элемент u₁₂ по формуле (3.11).

Далее (3 шаг) определяются элементы второй строки матрицы L в порядке, указанном стрелкой: l₂₁ и l₂₂ (по формулам (3.13) и (3.12) соответственно).

На 4 шаге выполняется расчет элементов 3 столбца матрицы U в порядке, обозначенном стрелкой: u₁₃, u₂₃ (формулы (3.11)) и т.д.

Рис. 3.1.

Рассмотрим теперь применение LU-разложения для решения СЛАУ вида

Ах = b,

где

A = LU.

Введем вспомогательный вектор у,

у = U × x.

(3.14)

Тогда исходную систему можно записать так

L × Ux = b ó Ly = b.

(3.15)

В силу формул (3.14) и (3.15) решение исходной СЛАУ сводится к последовательному решению систем (3.15) и (3.14) соответственно с верхней и нижней треугольной матрицами.

Метод прогонки.

Метод прогонки представляет собой вариант метода Гаусса, примененный к специальным системам линейных алгебраических уравнений, и учитывающий ленточную структуру матрицы системы. Пусть имеем СЛАУ со специальной трехдиагональной формой матрицы

c0y0 – b0y1 = f0, -aiyi-1 + ciyi – biyi+1 = fi, 1 £ i £ N-1; -aNyN-1 + cNyN = fN,

(3.16)

или в матричной форме: AY = F, где Y = (y₀, y₁,..., y_n)^T - вектор неизвестных; F = (f₀,f₁,...,f_n)^T - вектор правых частей; А - квадратная (N+1)´(N+1) матрица

Системы вида (3.16) возникают при конечно-разностной аппроксимации краевых задач математической физики, описываемых обыкновенными дифференциальными уравнениями второго порядка с постоянными и переменными коэффициентами, а также уравнениями в частных производных. Ставится задача разработать экономичные методы решения задач вида (3.16), число арифметических операций для которых пропорционально числу неизвестных. Таким методом для системы (3.16) является метод прогонки. Специфика матрицы А состоит в расположении ненулевых элементов, матрица А - разреженная матрица, из (N+1)² элементов которой ненулевыми являются не более 3N+1 элементов. Это позволяет получить для решения СЛАУ простые расчетные формулы.

Будем искать решение (3.16) в виде

yi = ai+1yi+1 + bi+1, i = N-1, N-2,..., 0

(3.17)

c неопределенными коэффициентами a_i, b_i. Выражение y_i-1 = a_iy_i + b_i подставим в (3.16)

(с_i - a_ia_i)y_i - b_iy_i+1 = f_i +a_ib_i

c учетом (3.17) имеем

[(с_i - a_ia_i)a_i+1 - b_i]y_i+1 + (c_i - a_ia_i)b_i+1 - a_ib_i = f_i .

Это равенство имеет место для любых y_i, если

(с_i - a_ia_i)a_i+1 - b_i = 0, (с_i - a_ia_i)b_i+1 - a_ib_i = f_i .

Отсюда получаем рекуррентные формулы для определения a_i+1, b_i+1

;	(3.18)
.	(3.19)

Коэффициенты a_i, b_i, i = 1, 2,..., N называются прогоночными.

Если коэффициенты a_i и b_i известны, а также известно y_N то, двигаясь справа налево (от i+1 к i) последовательно определяем все y_i. Задача нахождения a_i, b_i по формулам (3.18), (3.19) решается слева направо (от i к i+1). Начальные значения прогоночных коэффициентов a_i, b_i можно определить следующим образом. Полагаем в формуле (3.17) i=0, имеем y₀=a₁y₁+b₁, а из первого уравнения (3.16) c₀y₀ - b₀y₁ = f₀, откуда

.

(3.20)

Значение y_N определяется следующим образом. Полагаем в формуле (3.17) i=N-1, имеем y_N-1 = a_Ny_N + b_N, а из последнего уравнения (3.18) -a_Ny_N-1 + c_Ny_N = f_N откуда

.

(3.21)

Расчетные формулы (3.17) - (3.21) можно получить также из (3.16), если применить метод исключения Гаусса. Прямой ход метода заключается в том, что на первом шаге из всех уравнений системы (3.16) при помощи первого уравнения исключается y₀, затем из преобразованных уравнений для i=2,..., N при помощи уравнения, соответствующего i=1, исключается y_i и т.д. В результате получим одно уравнение относительно y_N. На этом прямой ход метода прогонки заканчивается. На обратном ходе для i=N-1, N-2,..., 0 находятся y_i.

Порядок счета в методе прогонки следующий:

1) исходя из значений a₁, b₁, вычисленных по формулам (3.20), все остальные коэффициенты a_i, b_i для i=2, 3,..., N-1 определяются последовательно по формулам (3.18) и (3.19);

2) исходя из значения y_N, рассчитанного по формуле (3.21), все остальные неизвестные y_i, i=N -1, N-2,..., 0 определяются последовательно по формуле (3.17).

Изложенный метод поэтому называется правой прогонкой.

Аналогично выводятся формулы левой прогонки:

	(3.22)
	(3.23)
yi+1 = xi+1yi + hi+1, i = N-1, N-2, …, 0; y0 = h0.	(3.24)

Здесь y_i находятся последовательно для значений i=1, 2,..., N; ход вычислений - слева направо.

В случае, если необходимо найти только одно неизвестное, например, y_m (0 £ m £ N) или группу идущих подряд неизвестных, целесообразно комбинировать правую и левую прогонки. При этом получается метод встречных прогонок.

Произведем подсчет числа арифметических действий для метода правой прогонки. Анализ формул (3.17)-(3.21) показывает, что общее число арифметических операций есть Q=8N+1. Коэффициенты a_i не зависят от правой части СЛАУ (3.16) и определяются только элементами a_i, b_i, c_i матрицы А. Поэтому, если требуется решить серию задач (3.16) с различными правыми частями, то прогоночные коэффициенты a_i вычисляются только для первой серии. Для каждой последующей серии задач определяются только коэффициенты b_i и решение y_i, причем используются ранее найденные a_i.

На решение первой из серии задач расходуется 8N+1 операций, а на решение каждой следующей задачи 5N+3 операций. Число арифметических операций, необходимое для решения СЛАУ (3.16) методом левой прогонки и методом встречных прогонок такое же, т.е. Q»8N. Метод правой прогонки будем называть корректным, если c_i-a_ia_i ¹ 0 при i=1,2,...,N.

Решение y_i находится по рекуррентной формуле (3.17). Эта формула может порождать накопление ошибок округления результатов арифметических операций. Пусть прогоночные коэффициенты b_i и a_i найдены точно, а при вычислении y_N допущена ошибка Dy_N, т.е. . При вычислениях с помощью формулы (3.17) мы получаем

.

(3.25)

Вычитая из (3.25) значение y_i по формуле (3.17) имеем для погрешности с заданным Dy_N. Отсюда ясно, что если a_i по модулю больше единицы, и если N достаточно велико, то вычисленное значение будет значительно отличаться от искомого решения y_i. В этом случае говорят, что алгоритм прогонки неустойчив.

Определение. Алгоритм прогонки называется устойчивым, если |a_i| £ 1, i = 1, 2, …, N.

Условия корректности и устойчивости алгоритма правой прогонки определяются следующей теоремой.

Теорема. Пусть коэффициенты системы (3.16) действительны и удовлетворяют условиям:

|b₀| ³ 0, |a_N| ³ 0, |c₀| > 0, |c_N| > 0, |a_i| > 0, |b_i| ³ 0, i = 1, 2, …, N-1;

|ci| ³ |ai| + |bi|, i = 1, 2, …, N-1;	(3.26)
|c0| ³ |b0|, |cN| ³ |aN|	(3.27)

причем хотя бы в одном из неравенств (3.26) и (3.27) выполняется строгое неравенство, т.е. матрица А имеет диагональное преобладание. Тогда для алгоритма (3.17)-(3.21) имеют место неравенства: c_i - a_ia_i ¹ 0, |a_i| £ 1, i = 1, 2, …, N, т.е. алгоритм метода правой прогонки корректен и устойчив.

Условия (3.26) и (3.27) теоремы обеспечивают также корректность и устойчивость алгоритмов левой и встречных прогонок. Эти условия сохраняются и для случая системы (3.16) с комплексными коэффициентами a_i, b_i, c_i.

Легко показать, что при выполнении условий (3.26)-(3.27) теоремы система (3.16) имеет единственное решение при любой правой части.