ИТЕРАЦИОННИ МЕТОДИ ЗА ЛИНЕЙНИ УРАВНЕНИЯ

Нека да сме благодарни на проф. Михаил Константинов за разрешението да четем "Елементи на линейната алгебра: вектори и матрици" online.

12.1 Уводни бележки

Редица практически задачи изискват решаването на големи линейни системи
Ах = b, A F^nxn,
например с n = 1000 или n = 10 000, в които матрицата А има специална структура. Типичен пример тук са т.нар. разредени матрици, които често се оказват и симетрични. В такива случай директните методи, основани на декомпозиция на матрицата А на множители, могат да се окажат неефикасни, тъй като те разрушават специалната структура на матрицата.

Алтернатива на директните методи са итерационните методи, при които решението се получава като граница на редица от последователни приближения. На свой ред тези приближения удовлетворяват някакви нови алгебрични уравнения, които са сравнително лесни за решаване.

12.2 Основни итерационни схеми

Да разгледаме уравнението
Ах = b, A F^nxn, b Fⁿ, (1)
относно неизвестния вектор х Fⁿ, в което матрицата А е неособена. Всички итерационни схеми се основават на еквивалентно представяне на уравнението във вида
Мх = Nх + с, (2)
където М F^nxn е неособена матрица с проста структура, която позволява лесно решаване на уравнението Мх = d. Например М може да се избере като диагонална матрица.

Уравнението (2) се решава итеративно по схемата
Мx^(k+1) = Nх^(k) + c, k = 0, 1, 2, ..., (3)
при зададено начално приближение ж⁽⁰⁾.

Когато спектралният радиус rad(М^-1N) на матрицата М^-1N е по-малък от 1, последователните приближения х^(k) се схождат при k → ∞ към вектора
х = (М - N)^-1c
за всеки начален вектор x⁽⁰⁾. При това

където В := |М^-1N| и |R| := [|r_i,j|] е матричният модул на матрицата R := [r_i,j]. Ще напомним, че неравенството x y за векторите x, у с елементи x_i, у_i съответно се разбира поелементно, т.е., х₁ ≤ у₁,... ,х_n ≤ y_n.

В сила е и оценката
|x - x^(k+1)| B^k|x - x⁽⁰⁾|,
в която обаче векторът |х — х⁽⁰⁾| е предварително неизвестен.

За да бъдат уравненията (1) и (2) еквивалентни, т.е., за да имат те едно и също решение, е необходимо и достатъчно да бъде изпълнено условието
с=(М - N)А^-1b
(докажете!). В частност може да се избере А = М — N и с = b. Това представяне се нарича разцепване на матрицата А. Разцепването се прави с оглед на следните съображения:
- матрицата М трябва да е неособена и проста (в частност лесна за обръщане), така че уравнението Мх = d да се решава лесно за всяка дясна страна d;
- спектралният радиус на матрицата М^-1N трябва да е колкото се може по-малък от 1.

Подобни изисквания се налагат и в общия случай на еквивалентно представяне на уравнението (1) във вида (2).

Да представим матрицата А = [a_i,j] във вида
А = L + D + U,
където L, D и U са съответно строго долно триъгълната, диагоналната и строго горно тригъгълната части на А,

D := diag(a_1,1,.....,a_n,n)

Да предположим, че матрицата има ненулев диагонал. Тогава можем да изберем
М = D, N = L + U, с = b.
Това е най-простата итерационна схема, известна като итерация на Якоби. При тази схема последователните приближения х^(k)_i към елементите х_i на решението х се определят по формулите
(4)
Както е прието, при i = 1 първата сума от 1 до 0 се приема за празна, а при i = n за празна се приема втората сума от n + 1 до n.

Итерацията на Якоби допуска известно подобрение както следва. От (4) се вижда, че при определянето на х^(k+1)_i, i ≥ 1, се използват приближенията х^(k)₁,..... ,х^(k)_i-1 (вж. първата сума), докато вече са определени и приближенията x^(k+1)₁,...,х^(k+1)_i-1, за които се предполага, че са по-точни. Като заместим в първата сума тези по-точни приближения получаваме итерацията на Гаус-Зайдел
(5)

Формално итерацията на Гаус-Зайдел се получава като вземем
М = L + D, N = -U.

Може да се покаже, че
rad(D^-1(L + U)) < 1
и следователно итерацията на Якоби е сходяща за всяко начално приближение ж⁽⁰⁾, когато са изпълнени условията

(в този случай се казва, че матрицата А е с доминиращ диагонал).

В много приложения, например при числено решаване на елиптични частни диференциални уравнения, матрицата А е симетрична (U = L^T) и положително определена. В този случай
rad((L + D)^-1L^T) < 1
и итерацията на Гаус-Зайдел е сходяща.

Критичен момент при използване на итерационните схеми е изборът на матриците М и N оглед на минимизиране на спектралния радиус
ρ := rad(М^-1N) < 1
и съответно повишаване на скоростта на сходимост на итерационния процес.

Един от начините за минимизиране на ρ е както следва. Матриците М и N и векторът с се представят като като аффини изрази, зависещи от един скаларен параметър ω F, например
М = М(ω) := D + ωL,
N = N(ω) := (1 - ω)D - ωU,
с = с(ω) := ωb, 0 ≠ ω F
Параметърът ω се нарича релаксационен параметър, а съответната итерационна схема
М(ω)х^(k+1) = N(ω)х^(k) + с(ω), k = 0, 1, 2 ..., (6)
- метод на последователната релаксация. При ω = 1 получаваме класическата итерация на Гаус-Зайдел.

Параметърът ω се определя от условието за минимум на спектралния радиус
ρ(ω) := rad(М^-1(ω)N(ω)).
Това е трудна задача в общия случай. За щастие при някои важни за практиката случаи на симетрична положително определена матрица А е възможно релаксационният параметър да се намери с прости алгебрични операции.

Съществува и една друга група итерационни методи за решаване на уравнението Ах = b с неособена матрица А. Те се основават на минимизирането (или максимизирането) на някаква нелинейна функция x → f(x), която има глобален екстремум в точката х = А^-1b. Нека например F = R и матрицата А е симетрична положително определена. Да разгледаме израза
f(x) = 0.5x^TАх — b^Tх.
Минимумът на f е — b^TA^-1b/2 и се достига за х = А^-1b.

Един от най-простите начини за минимизиране на f е по метода на най-бързото спускане. Нека е дадена точка x⁽⁰⁾, такава че Ах⁽⁰⁾ ≠ b⁵. Функцията f намалява най-бързо по посока на вектора r⁽⁰⁾ := b — Ах⁽⁰⁾ ≠ 0, противоположен на градиента Ах⁽⁰⁾ - b на f в точката x⁽⁰⁾. Следователно съществува число μ > 0, такова че
f(x⁽⁰⁾ + μr⁽⁰⁾) < f(x⁽⁰⁾).

Ако определим μ от условието за минимум на израза f(х⁽⁰⁾ + μr⁽⁰⁾) получаваме
μ = μ₀ = (r⁽⁰⁾)^Tr⁽⁰⁾/(r⁽⁰⁾)^ТAr⁽⁰⁾.
С тази стойност на μ векторът x се актуализира до
x⁽¹⁾ := x⁽⁰⁾ + μ₀r⁽⁰⁾.
Аналогично се определят следващите приближения
х^(k+1) := х^(k) + μ_kr^(k), k = 1,2,....

Глобалната сходимост на метода се гарантира от неравенството, дадено в упражнение 12.6. Тази сходимост, обаче, може да не е достатъчно бърза в случаите, когато числото на обусловеност
соnd₂(А) = λ_max(A)/λ_min(A)
на матрицата А е много по-голямо от 1. Поради това се използват и други, по-усъвършенствани итерационни методи за минимизиране на функцията f, които тук не се разглеждат.

12.3 Упражнения

Упражнение 12.1 Покажете, че итерацията
Мх^(k+1) = Мх^(k) + b, А = М - N
с неособена матрица М може да се запише във вида
х^(k+1) = х^(k) + М^-1r^(k), r^(k) := b - Ах^(k)

Докажете, че ако при неособена матрица А тази итерация е сходяща за всяко начално приближение х⁽⁰⁾ то непременно е изпълнено неравенството rad!(М^-1N) < 1. Вярно ли е това твърдение когато матрицата А е особена?

Упражнение 12.2 Нека матрицата А = М — N е особена, а матрицата М - неособена. Възможно ли е да е изпълнено неравенството rad(М^-1N) < 1? Разсъждавайте чрез допускане на противното.

Упражнение 12.3 Нека спектралният радиус rad(В) на матрицата В F^nxn е по-малък от 1. Покажете, че за всяко положително число ε < 1 — rad(В) съществува векторна норма || • || : Fⁿ → R₊, такава че
||В|| < rad(B) + ε,
където ||В|| е съответната операторна норма на матрицата В.

Упражнение 12.4 Покажете, че при метода на най-бързото спускане е в сила оценката
F(x^(k+1)) ≤ [1 - (1/cond₂(A))].f(x^(k)) - b^TA^-1b/2cond₂(A).

Изпратете материали(програми), свързани с математика на: