3.2. Алгоритм 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 и, выражая из последнего соотношения _ii, получаем:

.

(3.12)

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

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

.

(3.13)

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

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

Далее (3 шаг) определяются элементы второй строки матрицы 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) соответственно с верхней и нижней треугольной матрицами.

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

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

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 = i+1yi+1 + i+1, i = N-1, N-2, ..., 0

(3.17)

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

(с_i - a_i_i)y_i - b_iy_i+1 = f_i +a_i_i

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

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

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

(с_i - a_i_i)_i+1 - b_i = 0, (с_i - a_i_i)_i+1 - a_i_i = f_i .

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

;	(3.18)
.	(3.19)

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

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

.

(3.20)

Значение y_N определяется следующим образом. Полагаем в формуле (3.17) i=N-1, имеем y_N-1 = _Ny_N + _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. исходя из значений ₁, ₁, вычисленных по формулам (3.20), все остальные коэффициенты _i, _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 = i+1yi + i+1, i = N-1, N-2, …, 0; y0 = 0.	(3.24)

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

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

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

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

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

.

(3.25)

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

Определение. Алгоритм прогонки называется устойчивым, если |_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_i_i  0, |_i|  1, i = 1, 2, …, N, т.е. алгоритм метода правой прогонки корректен и устойчив.

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

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