Баландин М.Ю., Шурина Э.П. - Методы решения СЛАУ большой размерности, страница 6

При этом коэффициенты ортогонализации в(4.22) и (4.23) одинаковы.Из сказанного следует, что можно записать аналог матричнойформулы (4.20)TWmAVm = Tm ,(4.26)где Tm — симметричная9 трехдиагональная матрица, элементы которой определяются соотношениями[Tm ]ii = αi (vi , wi ) ;[Tm ]i+1,i = [Tm ]i,i+1 = βi+1 (vi , wi ) ,(4.27)(4.28)легко выводимыми из (4.22)–(4.23).Замечание. Основным недостатком биортогонализации Ланцоша является возможность возникновения ситуации, когда (v i , wi ) == 0; при этом продолжение процесса становится невозможным из-занеопределенности коэффициента βi+1 .8Иначе способ определения коэффициентов β сделал бы невозможным построение всех следующих векторов.9Симметричность следует из того, что для построения vi+1 и wi+1 используютсяодни и те же коэффициенты ортогонализации.39Глава 5Методы крыловского типаВ §4.1 был построен класс проекционных методов, основой которого служит условие Петрова-Галеркина (4.1) и формула (4.5).

В §4.2были приведены примеры построения нескольких методов для достаточно простого случая одномерных подпространств K и L.Однако выбор m = 1, как правило, дает медленную сходимостьи соответствующие методы оказываются пригодными лишь дляСЛАУ небольшой размерности или специальных случаев. В настоящее время широкое распространение получили проекционныеметоды, которые в качестве K используют описанные в §4.4 подпространства Крылова размерности больше единицы.

Такие методыназываются методами крыловского типа 1 . Для упрощения вычислений (WTAV)−1 из формулы (4.5) в них используются соотношения(4.19)–(4.21).5.1.FOM: Метод полной ортогонализацииНаиболее характерным примером применения описанного в §4.1подхода является метод полной ортогонализации, известный также как метод Арнольди или FOM2 .Положим K = L и будем проектировать искомое решение на подпространство Крылова Km (v1 , A), где в качестве v1 выступает пронормированный вектор начальной невязкиv1 = r0 /β ,β = kr0 k2 .(5.1)Тогда V = W, для построения базиса можно воспользоваться ортогонализацией Арнольди (§4.5), а матрица W TAV в силу (4.20) будет12Англ.

Krylov sequence methods.Full Orthogonalization Method.40r0 := b − Ax0β := kr0 k2 ; v1 := r0 /βСоздать (m × m)-матрицу Hm ; положить Hm = 0Создать m-вектор f ; положить f = βe1Для j = 1, mwj := AvjДля i = 1, jhij := (wj , vi )wj := wj − hij viувеличить ihj+1,j := kwj k2Если hj+1,j = 0тоm := jПрервать цикл по jvj+1 := wj /hj+1,jувеличить jДля i = 2, mОбнулить hi,i−1 , выполнив гауссовское исключениедля i-й строки системы Hm y = fотносительно (i − 1)-йувеличить iНайти y из верхнетреугольной системы H m y = fxm := x0 +mPi=1yi viРис. 5.1. FOM41иметь верхнюю хессенбергову форму и следовательно, окажется легко обратимой.

Решение СЛАУ с такой матрицей сведется к (m − 1)гауссовскому исключению строк и последующему обратному ходупо верхнетреугольной матрице.Присутствующий в (4.5) сомножитель W Tr0 в силу (5.1) и (4.21)будет равен βe1 . Таким образом, формула (4.5) будет сведена кxm = x0 + βVH−1m e1 .(5.2)На основании (5.2) и описанного в §4.5 алгоритма ортогонализации Арнольди, необходимого для построения базиса V, метод FOMможет быть записан как на рис.5.1.В такой формулировке метод обладает серьезным недостатком:заранее неизвестно, какой должна быть размерность подпространства m, чтобы найденное решение xm имело достаточную точность.Однако существует теорема, позволяющая оценить уменьшение нормы невязки без явного нахождения решения.Теорема 5.1.1 Норма вектора невязки kr m k2 в методе FOM удовлетворяет соотношениюkrm k2 = hm+1,m |ym | .(5.3)Доказательство.

Найдем невязку rm , воспользовавшись свойством(4.19)rm = b − A(x0 + Vm y) = r0 − AVm y == βv1 − (Vm Hm + wm eTm )y .Подставим сюда выражения y = βH−1m e1 и wm = hm+1,m vm+1 :rm = βv1 − βVm e1 − hm+1,m (em , y)vm+1 == −hm+1,m ym vm+1 .С учетом того, что векторы vj нормированы, а коэффициенты hj+1,jвычисляются как нормы векторов и, следовательно, неотрицательны, отсюда и следует утверждение теоремы.2Проверка условия выхода krk2 < ε в соответствии с (5.3) можетбыть достаточно просто введена в алгоритм FOM, если гауссовскоеисключение строк H проводить сразу же в цикле ортогонализацииАрнольди (при этом коэффициенты исключения необходимо хранить в отдельном векторе, так как матрица H по мере нахождениябазиса расширяется).42Однако на практике чаще применяется другой подход, гораздоболее простой с точки зрения программной реализации. Его сутьзаключается в том, что заранее выбирается некоторая размерностьподпространства m, относительно небольшая по сравнению с порядком СЛАУ.

После того как решение x m найдено, проверяется,является ли оно достаточно точным. Если необходимая точностьеще не достигнута, то весь процесс повторяется с вектором x m вкачестве начального приближения:Выбрать m < n и начальное приближение x 0Начало(алгоритм рис.5.1 для выбранного m)x0 := xmПовторять, пока kb − Ax0 k2 > εТакая схема носит название перезапускаемого FOM 3 или FOM(m).Замечание. При программной реализации FOM и FOM(m) основные затраты памяти приходятся на хранение базисных векторовv1 , . .

. , vm подпространства Крылова; они составляют O(mn) вещественных чисел.Еще одним вариантом метода является неполная ортогонализация, когда вектор wj+1 ортогонализуется не ко всем ранее найденным v1 , . . . , vj , а только к k предыдущим векторам vj−m+1 , . . . , vj . Этопозволяет уменьшить расходы памяти при программной реализации метода; соответствующая схема носит название IOM 4 [13].5.2.Предобусловливание в схеме методаНа примере построенного в предыдущем параграфе FOM покажем,как в схему проекционного метода может быть введено предобусловливание (2.1) без необходимости явного вычисления матричногопроизведения.Перепишем алгоритм FOM в применении к системе (2.1) с матрицей предобусловливателя M. Очевидно, что в нем необходимо изменить только те фрагменты, где фигурирует матрица СЛАУ и векторправой части:3Англ.

Restarted FOM.Incomplete Orthogonalization Method. Встречается также название TruncatedFOM.443Построить матрицу предобусловливателя Mz := b − Ax0Найти r0 из системы Mr0 = zβ := kr0 k2 ; v1 := r0 /βСоздать (m × m)-матрицу Hm ; положить Hm = 0Создать m-вектор f ; положить f = βe1Для j = 1, mz := AvjНайти wj из системы Mwj = zДля i = 1, jhij := (wj , vi )wj := wj − hij viувеличить ihj+1,j := kwj k2Если hj+1,j = 0тоm := jПрервать цикл по jvj+1 := wj /hj+1,jувеличить jДля i = 2, mОбнулить hi,i−1 , выполнив гауссовское исключениедля i-й строки системы Hm y = fотносительно (i − 1)-йувеличить iНайти y из верхнетреугольной системы H m y = fxm := x0 +mPi=1yi viРис.

5.2. Предобусловленный FOM44r0 := M−1 (b − Ax0 )·········Для j = 1, mwj := M−1Avj·········увеличить j·········Pxm := x0 + mi=1 yi viТаким образом, вектор начальной невязки r 0 и векторы wj находятся из решения систем Mr0 = b−Ax0 и Mwj = Avj соответственно.Напомним, что для матрицы M в §2.1 было сформулировано требование легкой обратимости. Для ILU-предобусловливания решение таких систем сводится к решению двух треугольных систем, для SSORпредобусловливания — к решению двух треугольных систем и умножения на диагональную матрицу и т.д.Следовательно, алгоритм FOM с введенным в него предобусловливанием матрицей M можно записать как на рис.

5.2.Подобным же образом предобусловливание может быть введенои в схему любого другого проекционного метода, требующего лишьвозможности вычисления матрично-векторных произведений безпересчета элементов матрицы СЛАУ.5.3.GMRES: Метод обобщенныхминимальных невязокПусть теперь пространства K и L связаны соотношением L = AK,причем в качестве K по-прежнему используется подпространствоКрылова Km (v1 , A) с выбором v1 = r0 /β, β = kr0 k2 .

Для построенияортонормального базиса K снова воспользуемся ортогонализациейАрнольди.Однако вместо проектирования (4.5) будем рассматривать эквивалентную (см. теорему 4.3.2) задачу минимизации функционалаΦ2 (x) = krx k22 на пространстве K:y = arg min kb − A(x0 + Vm y)k22 .yВычислительная схема, построенная с использованием такихпредпосылок, называется методом обобщенных минимальных невязок (GMRES)5 .5General Minimal RESiduals.45Выбрать m < n и начальное приближение x 0p := mr0 := b − Ax0 ; β := kr0 k2Создать ((m + 1) × m)-матрицу HСоздать (m + 1)-вектор gНачалоH := 0; g := βe1v1 := r0 /βДля i = 1, mw := AviДля k = 1, ihki := (w, vk )w := w − hki vkувеличить khi+1,i := kwk2vi+1 := w/hi+1,iДля k = 1, i − 1(h1i , . .