Г.М. Кобельков - Курс лекций по численным методам (1160467), страница 10
Текст из файла (страница 10)
Рассмотрим матрицу B:она, ясное дело, тоже будет симметрической. Перепишем систему в терминах погрешностей:rn+1 − rn+ Arn = 0.(86)Bτ√ −1√Теперь сделаем замену переменной: положим yn := Brn , D :=B. Тогда получим следующее:yn+1 − ynt+DAD} yn = 0.| {zτ(87)CСовсем легко проверить, что матрица C тоже будет симметрической и положительно определённой.
Такимобразом, задача сведена к предыдущей.А чем это нам помогает? А вот чем. Имеет место равенствоλ(А)min = min(Dt ADy, y),(y, y)λ(А)max = max(Dt ADy, y).(y, y)(88)Сделаем замену z := Dy, тогда получим(Dt ADy, y)(Az, z)=.(y, y)(Bz, z)Определение. Отношение(Ax,x)(x,x)(89)называется отношением Релея.Поскольку матрицы A и B задают какие-то нормы, а в конечномерном пространстве все нормы эквивалентны, всегда найдутся такие константы γ1 и γ2 , чтоγ1 (By, y) 6 (Ay, y) 6 γ2 (By, y).(90)Отсюда следует, что λ(A)min > γ1 и λ(A)max 6 γ2 .
Отсюда следует, что если матрицу B выбрать «близкой» кматрице A, то после переобуславливания мы получим систему с хорошим спектром. Оценим скорость сходимости: имеем√√22kyn k = ( Brn , Brn ) = krn kB ,(91)кроме того, применяя те же методы оценки, что и раньше, получаемγ2 − γ1kyn+1 k2 6kyn k2 ,γ2 + γ1поэтому получаем оценкуγ2 − γ122krn+1 kB 6krn kB .γ2 + γ1(92)(93)Пример 3.1. Пусть Jk (λ) — матрица, у которой числами λ заполнена только k-я строка над диагональю(k = 0 — сама диагональ).
Положим Sk (λ) := 21 (Jk (λ) + Jkt (λ)). Рассмотрим матрицу A := S0 (3)+S1 (1)+S2 12 ++ S4 21 , у неё будет, понятное дело, 7 ненулевых диагоналей. Тогда положим B := S0 (3) + S1 (1) (она будеттрёхдиагональной). Несложно показать, что1(Bx, x) 6 (Ax, x) 6 2(Bx, x),22а число τ , стало быть, можно брать равным γ1 +γ= 2+2 1 = 45 .2236(94)4.3.4. Метод ЗейделяЭтот метод годится для решения систем Ax = b, у которых aii 6= 0. Представим матрицу A в виде A = B + C,где С — строго верхнетреугольная (на диагонали стоят нули), а B — нижнетреугольная. В силу условия надиагональные элементы, матрица B обратима.Процесс выглядит так:Bxn+1 + Cxn = b,(95)или, переходя на язык погрешностей,rn+1 = −B −1 Crn .Ясно, что сходимость такого метода будет тогда и только тогда, когда B −1 C < 1.Теорема 4.4.
Если A = At > 0, то метод сходится.(96)Замечание. Условие на диагональные элементы обеспечивается автоматически свойством A > 0. Мы покажем, что метод Зейделя и метод покоординатного спуска — это одно и то же. РассмотримфункционалF (x) = (Ax, x) − 2(b, x).(97)Как мы уже выясняли при изучении метода градиентного спуска, (единственный) минимум этого функционалаи есть искомое решение системы. На каждом шаге будем постепенно уточнять значения координат решения(xn1 , .
. . , xnm ), получая в конце шага новый вектор (xn+1, . . . , xn+1m ).1Пусть мы уже уточнили координаты с номерами 1 по k − 1. Уточним k-ю координату. Продифференцируемфункционал F по переменной xk и приравняем её к нулю, получимk−1mXX∂Fakj xnj − bk xk = 0,(98)= 2akj xn+1+ akk xk +j∂xkj=1j=k+1то есть это действительно покоординатный спуск. Теперь осталось понять, почему метод сходится. Метод работает таким образом, что на каждом шаге имеет место неравенство F (xn+1 ) < F (xn ). Пользуясь формулой (58),получаем, что его можно переписать в видеkxn+1 − xk2A2kxn − xkAПереходя к погрешностям, получаем2krn+1 kAkrn k2A=(99)< 1. −1B Crn 2Akrn k2A< 1.(100)Выражение в левой части не меняется при замене rn на crn , поэтому можно считать, что krn kA = 1. Такимобразом, получаем непрерывную функцию на единичной сфере, поэтому она достигает там своего максимумаM < 1. Значит, сходимость будет, но оценить её скорость невозможно.
4.3.5. Метод сопряжённых градиентовРассматриваем систему Ax = b с симметричной положительно определённой матрицей. Снова рассматриваемфункционалF (y) := (Ay, y) − 2(b, y).(101)ИмеемF (y) = (Ay, y) − (b, y) − (b, y) + (Ax, x) − (Ax, x) == (Ay, y) − (Ax, y) − (x, Ay) +(Ax, x) − (Ax, x) = (A(y − x), y − x) − (Ax, x) =| {z }(Ax,y)222= A−1 (Ay − b), Ay − b − kxkA = kAy − bkA−1 − kxkA . (102)Положим ξn := xn − x. Обычный итерационный процесс устроен примерно так: ξn = Pn (A)ξ0 , где Pn (A) —некоторый многочлен степени n, у которого Pn (0) = 1. Применим к ξn оператор A, получимAξn = APn (A)ξ0 = Pn (A)Aξ0 = Pn (A)r0 ,а с другой стороны Aξn = rn , поэтому Pn (A)r0 = rn .37(103)Пусть задано начальное приближение x0 . Хотим построить процесс, который при всяком n использует многочлен Pn (A), зависящий, вообще говоря, от x0 , такой что(104)krn kA−1 → min .PnПоложим vi := Ai r0 .
Минимизируем krn kA−1 (сейчас мы заново докажем, что перпендикуляр к подпространству — это минимальное расстояние до него). Пусть c0 , . . . , cn — коэффициенты многочлена Pn (здесь c0 = 1).Тогда имеемnXrn =ci vi .(105)i=0Продифференцируем по ck (k = 1, . . .
, n): ∂ −1∂A rn , rn + A−1 rn ,rn =∂ck∂ck−1−1= A vk , rn + A rn , vk = A−1 vk , rn + rn , A−1 vk = 2 A−1 vk , rn . (106)∂∂2krn kA−1 =(A−1 rn , rn ) =∂ck∂ckПриравнивая частные производные к нулю, находим, что rn ⊥A−1 vk = vk−1 = Ak−1 r0 при k = 1, . . . , n.Определение. Подпространство Крылова для оператора A = At > 0, порождённое вектором e ∈ V — этовот такое циклическое подпространство:e, Ae, A2 e, . . . .(107)Мы знаем, что у всякого симметрического оператора существует собственный базис. Пусть r0 как-то разложен по собственным векторам {ei } оператора A, отвечающим различным собственным значениям:r0 =qXr i ei ,(108)i=1причём все ri отличны от нуля.
Так всегда можно сделать, ибо линейная комбинация векторов из собственногоподпространства тоже является собственным вектором с тем же собственным значением, значит, если в разложении по собственному базису вектора r0 встретилось несколько векторов с одним собственным значением, тоих можно объединить в один.Покажем, что векторы r0 , Ar0 , A2 r0 , . . . , Aq−1 r0 линейно независимы. В самом деле, допустим, что имеетсянетривиальная линейная зависимость между ними:0=q−1Xcj Aj r0 =j=0q−1Xcjj=0qXri λji ei =i=1q Xq−1Xi=1j=0cj ri λji ei .(109)Но ei отвечают различным собственным значениям и потому линейно независимы.
Значит,q−1Xcj ri λji = 0,(110)i = 1, . . . , q.j=0Это линейная система q × q на cj с матрицей W (λ1 , . . . , λq )·diag(r1 , . . . , rq ), а поскольку λi различны, то матрицаневырождена. Значит, у неё есть только нулевое решение, поэтому все cj = 0.Отсюда следует корректность метода, поскольку доказанная линейная независимость означает, что rn (или,что то же самое, коэффициенты многочлена Pn ) находится из системы уравнений (rn , Ak−1 r0 ) однозначно.Рассмотрим подпространства Ln := hr0 , .
. . , An r0 i и Rn = hr0 , . . . , rn i. В силу того, что rn ⊥Ak r0 при k =0, . . . , n − 1, получаем, что векторы ri тоже линейно независимы, значит, они образуют базис пространства Rn .А поскольку rk = Pk (A)r0 , то они выражаются через базис пространства Ln . Отсюда следует, что Rk = Lk привсех k. Мы уже знаем, что rn ⊥Ln1 = Rn−1 , стало быть, rn ортогонален всем предыдущим ri .Покажем, что Rn = hr0 , . .
. , rn−1 , Arn−1 i. В самом деле,Arn−1 = APn−1 (A)r0 = cn−1 An r0 + wn−1 .| {z }| {z }6=0(111)∈Rn−1Теперь разложим rn по этому базису:rn =n−1Xγk rk + γn Arn−1 .k=038(112)Из взаимной ортогональности rj следует, что при j < n − 2 имеем γj = 0. В самом деле,0 = (rn , rj ) = γj (rj , rj ) + γn (Arn−1 , rj ).(113)Второе скалярное произведение равно нулю: (Arn−1 , rj ) = (rn−1 , Arj ) = 0, потому что Arj ∈ Rn−2 . Значит,γj = 0.Теперь вычислим три старших коэффициента, то есть γn−2 , γn−1 и γn . Имеем систему из двух уравнений стремя неизвестными:20 = (rn , rn−2 ) = γn−2 krn−2 k + γn (rn−1 , Arn−2 ),(114)20 = (rn , rn−1 ) = γn−1 krn−1 k + γn (rn−1 , Arn−1 ).Нужно ещё одно.
Мы уже получили, чтоrn = γn−2 rn−2 + γn−1 rn−1 + γn Arn−1 .(115)Приравняем коэффициенты при r0 в разложении этого вектора по базису r0 , . . . , An r0 , получим ещё одно уравнение. ИмеемPn (A)r0 = γn−2 Pn−2 (A) + γn−1 Pn−1 (A) + γn APn−1 (A)r0 .(116)Последнее слагаемое свободного члена не имеет. Поскольку при всех k имеем Pk (0) = 1, получаем ещё однонедостающее уравнение 1 = γn−2 + γn−1 .Решая полученную систему из трёх уравнений, находим γn−2 , γn−1 и γn .
А теперь остаётся только получитьформулу для решения (то есть от невязок перейти к векторам xn ):Axn − b = γn−2 (Axn−2 − b) + γn−1 (Axn−1 − b) + γn A(Axn−1 − b),(117)а в силу условия γn−2 + γn−1 = 1 это можно упростить:Axn = γn−2 Axn−2 + γn−1 Axn−1 + γn A(Axn−1 − b).(118)Осталось домножить слева на матрицу A−1 , и мы получим замечательную формулу:xn = γn−2 xn−2 + γn−1 xn−1 + γn (Axn−1 − b).(119)Это и есть рекуррентная формула для вычисления итераций решения.В силу конечномерности пространства, метод сойдётся не более, чем за rk A шагов (это и называется конечностью метода). На самом деле шагов будет ровно q (напомним, что через q было обозначено количествособственных векторов, участвующих в разложении r0 ).4.4.
Что делать, когда всё плохо?4.4.1. Метод регуляризации по ТихоновуПусть, как обычно, у нас есть система Ax = b, причём A = At > 0 и x — точное решение. Будем предполагать,что λmax ∼ 1, а λmin ≪ 1. Более точно, будем предполагать, что λ1 > λ2 > . . . > λm , причём первые n чиселимеют порядок O(1).Пусть e1 , . .
. , em — собственный ортонормированный базис оператора A, и пустьb=mXbk ek .(120)k=1Сделаем предположение относительно решения:x=nXbkek ,λk(121)k=1то есть на самом деле bn+1 = · · · = bm = 0.Поскольку в мире ничего идеального не бывает, правая часть системы на самом деле задана с погрешностью,то есть eb = b + δb (здесь δb — единый символ). Стало быть, на самом деле мы решаем систему Ax = eb, и у неёимеется точное решение xe.Будем оценивать погрешность, которая появляется от замены b на eb и пытаться её уменьшить.
Имеемxe=mXbk + δbkk=1λkek =nmXXbkδbkek +ek .λkλkk=1k=1| {z }x39(122)Во второй сумме в знаменателе стоят маленькие числа. Поэтому, даже если погрешность δb невелика, приделении почти на ноль она вырастет, и всё испортится.Поэтому применяется такая идея. Давайте «испортим» матрицу A. А именно, заменим матрицу A на матрицуAα := A + αE (как говорят, рассмотрим возмущение). Тогда, как несложно видеть,xα =mXbk + δbkλk + αk=1Вычислимxα − x =mXbk + δbkk=1λk + αek −(123)ek .nnmXXXbk11δbkek =bk−ek +ek .λkλk + α λkλk + αk=1k=1k=1|{z} |{z}=:ϕ1(124)=:ϕ2Заметим, что поскольку λ1 , . . . , λn имеют порядок O(1), разность11α−=−λk + α λk(λk + α)λkимеет порядок O(α).
Отсюда2kϕ1 k =nXb2kk=1α226 C kbk α2 ,+ α)2λ2k (λk(125)(126)где константа C1 как-то зависит от спектра матрицы A.Для оценки ϕ2 вообще забудем про то, что у нас в знаменателях стоят λk (мысленно положим их нулями),от этого сумма только возрастёт. Тогда получим оценку2kϕ2 k 6Отсюда, полагая C :=регуляризации:12kδbk .α2(127)√C1 , получаем оценку для разницы между точным решением и тем, что получится после1kδbk .αОсталось её минимизировать за счёт выбора α (обозначим это оптимальное решение α∗ ):kex − xk 6 C kbk α +C kbk −откудаα∗ =skδbk,C kbkkδbk= 0,α2∗pkex − xk 6 2 C kbk · kδbk.(128)(129)(130)4.4.2. Метод Поспелова для решения плохо обусловленных системЭтот метод применяется для решения несовместных систем.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
