Решение системы линейных алгебраических уравнений с помощью LU-разложения

LU-метод основан на том, что если главные миноры матрицы А отличны от нуля, тогда матрицу А можно представить, причем единственным образом, в виде произведения

А=LU,

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

Рассмотрим СЛАУ Ax=f.

A=LU

где

или

Окончательно запишем

Полагая получим следующие рекуррентные формулы для вычисления элементов матрицы L и U

(11)

Если найдены матрицы L и U, то решение исходной системы A=LU сводится к последовательному решению двух систем уравнений с треугольными матрицами

LU – метод позволяет вычислить определитель матрицы А

В среде MathLab разложение матрицы А на треугольные матрицы L и U осуществляется с помощью команды

[L,U,P]=lu(A),

которая кроме нижней треугольной матрицы L и верхней треугольной матрицы U порождает и матрицу перестановок Р.

Метод простой итерации решения систем линейных алгебраических уравнений

Исходную систему (1) преобразуем к виду:

(12)

где , i=1,2,...,n; , a_ii¹0.

Перепишем систему (12) в виде одного матричного уравнения:

(13)

Первая сумма равна нулю, если верхний предел суммирования меньше нижнего.

Так (12) при i=1 имеет вид

Итерационный процесс (12) начинается с начальных значений , которые в общем случае задаются произвольно, но предпочтительнее, если за взять свободные члены исходной системы, т.е . Подставляя начальные значения в правую часть выражения (12) и вычисляя полученное выражение находим первое приближение:

Подставляя приближение в правую часть системы (12) получим:

Продолжая этот процесс получим последовательность , , ,…, ,… приближений вычисляемых по формуле:

, где ` (14)

В развернутой форме записи формула (14) выглядит так:

Условие окончания счета:

где i= .

Теорема

Если какая-либо норма матрицы меньше единицы: , то уравнение (13) имеет единственное решение ξ, к которому стремится последовательность итераций (14) при любом выборе начального приближения .

Число итераций для достижения заданной точности ε определяется из неравенства:

. (15)

Неравенство (15) обычно дает завышенную оценку числа итераций k.

Условие, позволяющее принять приближение в качестве решения с точностью ε, представляется в следующей форме:

. (16)

Метод вращения

Метод вращения, как и метод Гаусса состоит из прямого и обратного ходов.

Прямой ход метода. Исключаем неизвестное из всех уравнений, кроме первого. Для исключения из 2-го уравнения вычисляют числа

, ,

где a и b такие, что , .

Первое уравнение системы заменяем линейной комбинацией первого и второго уравнений с коэффициентами и ,а второе уравнение такой же комбинацией с и - . В результате получим систему

. (17)

Здесь

где j= .

Преобразование системы (1) к системе (17) эквивалентно умножению матрицы A и вектора B слева на матрицу С₁₂ вида

Аналогично для исключения из третьего уравнения вычисляем числа

и ,

такие, что .

Затем первое уравнение системы (17) заменяем линейной комбинацией первого и третьего уравнений с коэффициентами , , а третье уравнение системы (17) тоже заменяем линейной комбинацией с ,– . Это преобразование эквивалентно умножению слева на матрицу

Исключая неизвестное х₁ из всех последующих уравнений получим систему

А⁽¹⁾ х =В,

где матрица A⁽¹⁾=C₁_m…C₁₃C₁₂A,, а вектор правых частей .

Здесь и далее через С_kj обозначена матрица элементарного преобразования, отличающаяся от единичной матрицы Е только четырьмя элементами. Действие матрицы С_kj на вектор x эквивалентно повороту вектора x вокруг оси перпендикулярной плоскости на угол такой, что

, .

Операцию умножения на матрицу С_kj называют плоским вращением или преобразованием Гивенса.

Первый этап состоит из m-1 шагов, в результате чего получается система

(18)

В матричной форме А⁽¹⁾ х =В⁽¹⁾.

На втором этапе, состоящем из m-2 шагов, из уравнений системы (2.7) с номерами 3,4,…,m исключают неизвестное х₂. B результате получим систему

В матричной форме получаем , где , .

После завершения (m-1)-го шага придем к системе с верхней треугольной матрицей вида

где .

Обратный ход метода вращения проводится точно также как и для метода Гаусса.

Метод отражения

Метод отражения основан на преобразованиях Хаусхолдера. Алгоритм Хаусхолдера приводит n´n симметричную матрицу к трехдиагональной форме, за n–2 ортогональных преобразований. Каждое преобразование обнуляет необходимые части выбранного столбца и соответствующей строки. В этом методе QR -разложение матрицы A системы уравнений Ax = b

производится при помощи матриц отражения, которые имеют вид

U = E − 2 ww^T.

Здесь E – единичная матрица; w – n -мерный вектор единичной длины,
(w, w) = 1, а ww^T – квадратная симметричная матрица:

Следовательно, и матрица отражения U является симметричной, U = U^T. Более того, матрица U ортогональна, то есть U^T = U −1, действительно:

UU^T = (E − 2 ww^T)(E − 2 ww^T) ^T = E − 2 ww^T − 2 ww^T + 4 ww^T ww^T = E,

поскольку w^Tw = (w, w) = 1.

Так как U симметрична и ортогональна, U ² = UU^T = E, собственные числа матрицы U удовлетворяют условию , то есть , причем отрицательному собственному значению λ = −1 соответствует собственный вектор w:

Uw = λ w = Ew − 2 ww^T w = w − 2 w = − w = (−1) w → λ = −1.

Положительному собственному значению λ = 1 соответствуют все векторы, ортогональные вектору w. Если произвольный вектор v ортогонален вектору w, то есть (v, w) = 0, то

Uv = Ev − 2 ww^T v = v − 2 w (w, v) = v.

Рассмотрим действие матрицы U на произвольный вектор y, который представим в виде суммы двух ортогональных компонент: y = z + v. Компонента z направлена вдоль вектора w и является проекцией y на w:
z = dw, d = (y, w), а компонента v ортогональна этому вектору: (v, w) = 0, и равна v = y −(y, w) w. Тогда

Uy = U (z + v) = E (z + v) − 2 ww^T (z + v) = z + v − 2 ww^T z − 2 ww^Tv =

= z + v − 2 ww^Tdw = z + v − 2 z = − z + v.

Таким образом, вектор Uy является зеркальным отражением вектора y относительно плоскости, ортогональной вектору w. Используя это свойство матрицы отражения, можно подобрать вектор w таким, чтобы исходный вектор y ≠ 0 в результате отражения Uy получил направление некоторого заданного единичного вектора e. В результате отражения получается вектор

Uy =α e или Uy = −α e, α = (y, y), (19)

поскольку при ортогональных преобразованиях длины векторов сохраняются (U – ортогональная матрица). Направление, перпендикулярное к плоскости отражения, будет определяться вектором (y −αe) или вектором (y +αe).

Таким образом, векторы

или , (20)

где , , являются требуемыми компонентами матрицы отражения. Если векторы y и e параллельны, то отражения делать не надо (при этом или будут равны нулю).

Покажем, что произвольную квадратную матрицу можно представить в виде произведения ортогональной и правой (верхней) треугольной матриц.

Пусть A – квадратная матрица порядка n.

Приведем ее к правой треугольной форме путем последовательного умножения слева на ортогональные матрицы отражения. На первом шаге приведения в качестве вектора y рассмотрим первый столбец матрицы A:

Если (i= 2,3,…,n), следует перейти к следующему шагу приведения, положив A ⁽¹⁾ = A; U = E ₁ и обозначив . В противном случае умножим матрицу A слева на матрицу отражения U₁=E_n−2w₁w₁^T, где w₁ подбирается таким образом, чтобы вектор U₁y₁ стал параллелен вектору e₁= [1, 0, 0,, 0] ^T, то есть в соответствии с формулами (19), (20). Здесь E_n – единичная матрица порядка n, а e ₁, y₁ – n –мерные векторы. В результате такого преобразования в первом столбце матрицы A ⁽¹⁾ все элементы, кроме первого, станут равными нулю.

На втором шаге приведения преобразуемым вектором является второй

столбец матрицы A⁽¹⁾ без первого члена

Преобразование отражения выполняется умножением матрицу A ⁽¹⁾ слева на матрицу

где S_n_-1=E_n_-1-2w_n_-1w_n_-1^T – (n −1) -мерный вектор, вычисляющийся по формулам (19), (20). Тем самым обнуляются элементы второго столбца, расположенные ниже главной диагонали матрицы A⁽²⁾ =U₂ A⁽¹⁾.

Последующие шаги процесса приведения матрицы A проводятся аналогично. После выполнения k -го шага получается матрица A ⁽ ^k ⁾, все элементы которой, находящиеся ниже главной диагонали вплоть до k -го столбца матрицы, равны нулю: при i>j, j=1,2,…,k.

Для выполнения (k +1) -го шага приведения преобразуем вектор

Если компоненты вектора для (для i =(k+ 2), (k +3),…, n), получаем A ^(k⁺¹⁾ = A ^(k); U_k₊₁ = E_n и переходим к следующему шагу. В противном случае строим матрицу отражения

S_k₊₁=E_n_-_k-2w_k₊₁w_k₊₁^T

(вектор w_k₊₁ и матрица S_k₊₁ порядка (n − k)) для преобразования вектора y_k₊₁ в вектор, параллельный вектору e_k₊₁= [1, 0, 0,, 0] ^T (длины (n − k)), и переходим от матрицы A ^(k) к матрице A ^(k⁺¹⁾:

A^(k⁺¹⁾ = U_k₊₁_A^(k), где

Процесс этот всегда осуществим, и после (n −1) -го шага приходим к матрице

имеющей правую треугольную форму. Обозначив через U произведение матриц вращения , это выражение можно записать в виде A ⁽ⁿ⁻¹⁾= UA, или A = QR, где Q = U^T – ортогональная матрица, а R = A ⁽ⁿ⁻¹⁾ – правая треугольная матрица.

Решение системы Ax = b посредством метода отражения выполняется следующим образом. Умножив систему слева на последовательность матриц

отражения, сводим ее к виду с верхней треугольной матрицей

UAx = Ub Þ Rx = Ub, или A ⁽ⁿ⁻¹⁾ x = Ub.

Если все диагональные элементы матрицы A⁽ⁿ⁻¹⁾ отличны от нуля, то неизвестные x_i для i = n, (n −1),…,1 находятся, как и в методе Гаусса, обычным обратным ходом. Если же хотя бы один из диагональных элементов

матрицы равен нулю, то система A ⁽ⁿ⁻¹⁾ x = Ub вырождена, и в силу эквивалентности вырождена и исходная система.

Этот метод в настоящее время считается одним из наиболее устойчи-

вых к вычислительной погрешности, но более трудоемок в сравнении с методом Гаусса. Для получения QR -разложения квадратной матрицы A порядка n общего вида требуется около (4 / 3) n ³ арифметических операций.

Метод квадратного корня

Метод квадратного корня по своему идейному содержанию близок к LU-методу. Основное отличие в том, что он дает упрощение для симметричных матриц.ID_1

Этот метод основан на разложении матрицы А в произведение

где S–верхняя треугольная матрица с положительными элементами на главной

(22)

Из условия (2) получаются уравнения

Так как матрица А симметричная, не ограничивая общности, можно считать, что в системе (23) выполняется неравенство i≤j. Тогда (23) можно переписать в виде

В частности, при i=j получится

(24)

Далее, при i<j получим

По формулам (24) и (25) находятся рекуррентно все ненулевые элементы матрицы S.

Обратный ход метода квадратного корня состоит в последовательном решении двух систем уравнений с треугольными матрицами.

Решения этих систем находятся по рекуррентным формулам

Всего метод квадратного корня при факторизации A=S^TS требует примерно операций умножения и деления и m операций извлечения квадратного корня.