2.2 Методы решения сеточных уравнений

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

Рассмотрим численные методы решения системы линейных уравнений вида:

Ax  = b,
(2.103)

где A - матрица коэффициентов, b - вектор правой части, x - вектор неизвестных.

2.2.1 Предобуславливание и LU-факторизация

Пусть M : det(M)0, т.е. M - невырожденная матрица. Домножим (2.103) на M-1, получим

  -1       -1
M   Ax = M   b
(2.104)

Введя обозначения = M-1A,ˆ
 b = M-1b, уравнение (2.103) запишется в виде

     ˆ
ˆAx = b
(2.105)

Хотя (2.103) эквивалентна уравнению (2.105), но спектральные характеристики A, отличаются, что ведет к изменению скорости сходимости численных методов.

Процесс перехода от (2.103) к (2.105) с целью улучшения характеристик матрицы для ускорения сходимости называется предобуславливанием, а матрица M - матрицей предобуславливания.

Но практике часто используют другой подход. Поскольку невязка ˆr системы (2.105) связана с невязкой r системы (2.103) соотношением

M ˆr = r,
(2.106)

которое справедливо и для z = Aq и = q:

    ˆ      -1
ˆz = Aq = M   Aq → M ˆz = Aq = z
Это позволяет вместо явного перехода от (2.103) к (2.105) вводить в схемы методов корректирующие шаги для учета влияния предобуславливающей матрицы.

LU-факторизация позволяет представить матрицу A в виде

A = LAUA,
(2.107)

где LA и UA – нижне- и верхнетреугольные матрицы соответственно.

Подобное представление позволяет легко решать СЛАУ вида Ax = b путем выполнения прямого хода для нижнетреугольной системы LAy = b и затем обратного хода для верхнетреугольной системы UAx = y. Однако алгоритм факторизации непригоден для СЛАУ с разреженными матрицами, так как ведет к заполнению портрета, т.е. появлению в матрицах LA и UA ненулевых элементов в тех позициях, для которых aij = 0, и как следствие, резкому увеличению объема памяти, требуемой для хранения матриц.

Вместо задачи нахождения факторизации (2.107) сформулируем другую задачу. Для заданной матрицы A потребуем представить ее в виде

A = LU + R,
(2.108)

где матрицы в правой части такие что PL PA и PU PA.

Тогда приближенное представление A LU называется неполной LU-факторизацией матрицы A или коротко ее ILU-разложением.

2.2.2 Классические итерационные методы и релаксация.

Исторически первые итерационные методы основывались на циклическом покомпонентном изменении вектора решения, осуществляемом таким образом, чтобы обнулить соответствующий коэффициент вектора невязки и тем самым уменьшить его норму. Подобная методика уточнения решения получила название релаксации.

Хотя в настоящее время такие методы в их классической формулировке уже практически не применяются, существуют определенные классы задач, для которых разработаны их модификации, хорошо себя зарекомендовавшие. Кроме того эти методы могут быть применены (и применяются) не в качестве самостоятельного средства решения СЛАУ, а для предобусловливания.

Рассмотрим численное решение системы (2.103). Представим матрицу A в виде разности трех матриц

A = D - E - F.
(2.109)

где матрица D содержит диагональные элементы матрицы A, E - только поддиагональные, а F - только наддиагональные элементы.

Запишем следующий итерационный процесс:

M xk+1 = N xk + b, A = M - N,
(2.110)

Тогда при MJAC = D получаем метод Якоби, а при MGS = D -E - метод Гаусса-Зейделя.

Для ускорения скорости сходимости введем некоторый вещественный параметр ω и вместо системы (2.103) рассмотрим масштабированную систему

ωAx  = ωb.
(2.111)

и вместо (2.109) воспользуемся представлением

ωA  = (D - ωE )- (ωF + (1- ω)D ),
(2.112)

где матрицы D, E и F имеют тот же смысл, что и в (2.109).

Тогда на основании (2.111) и (2.112) можно построить итерационную схему, похожую на метод Гаусса-Зейделя,

(D  - ωE )xk+1 = (ωF + (1- ω )D )xk + ωb.
(2.113)

Схема (2.113) называется методом последовательной верхней релаксации (SOR). Для нее

MSOR  = (D - ωE ),
NSOR = (ωF + (1- ω)D ).
(2.114)

2.2.3 Связь с предобуславливанием

Рассмотрим предложенные методы в качестве предобуславливателя. Запишем метод Якоби, Гаусса-Зейделя и SOR в следующей форме:

xk+1 = Gxk + f,
(2.115)

которое получается из (2.110) домножением левой и правой частей на M-1.

Уравнений (2.115) можно рассматривать как итерационную схему, построенную для решения СЛАУ

(I - G)x = f,
(2.116)

где

      -1      -1                -1
G = M   N = M   (M  - A ) = I - M A,
f = M -1b.

Таким образом, система (2.116) эквивалентна системе

M -1Ax = M -1b
(2.117)

т.е. предобусловленной СЛАУ (2.103) с матрицей предобусловливания M. Для рассмотренных методов эти матрицы равны:

MJA  = D,
MGS  = D - E,
        1-
MSOR =  ω(D - ωE ).

2.2.4 Методы подпространств Крылова

Рассмотрим численные методы решения системы Ax = b, где A в общем случае несимметричная, положительно определенная, сильно разреженная матрица. Наиболее эффективными и устойчивыми среди итерационных методов решения систем линейных уравнений с разреженными матрицами являются так называемые проекционные методы, и особенно тот их класс, который связан с проектированием на подпространства Крылова. Эти методы обладают целым рядом достоинств: они устойчивые, допускают эффективное распараллеливание, работу с различными строчными (столбцовыми) форматами и предобуславливателями разных типов.

Введем два подпространства K Rn,L Rn, и будем искать решение x K оптимальное относительно подпространства L, т.е. чтобы выполнялось условие:

∀l ∈ L : (Ax,l) = (b,l)
называемое условием Петрова-Галеркина. Сгруппировав обе части равенства по свойствам скалярного произведения, перепишем это условие в виде:
∀l ∈ L : (rx,l) = 0
т.е. rx = b - Ax L. Такая задача называется задачей проектирования решения x на подпространство K ортогонально подпространству L.

Запишем в более общей постановке. Пусть для исходной системы (2.103) известно некоторое приближение x0 к решению ˜x. Требуется уточнить его поправкой δ K таким образом, чтобы b - A(x0 + δ) L. Условие Петрова-Галеркина в этом случае можно записать в виде:

∀l ∈ L : (r0 - A δ,l) = 0.

Введем в подпространства K и L базисы V = (vi)i=1m и W = (wi)i=1m соответственно. Поскольку можно записать δ = V y, где y Rm - коэффициенты, тогда из условия ортогональности запишем

W T (r0 - AV y) = 0,
откуда WTAV y = WTr0 и
       T    -1  T
y = (W  AV )  W  r0.
Таким образом, решение должно уточняться в соответствии с формулой
˜x = x0 + V (W TAV )- 1W Tr0,
из которой сразу вытекает важное требование: в практических реализациях проекционных методов подпространства K и L и их базисы должны выбираться так, чтобы матрица WTAV либо была малой размерности, либо имела простую структуру, удобную для обращения.

Пусть имеется набор пар подпространств (Ki,Li)i=1q. Последовательное применение описанного ранее процесса ко всем таким парам приведет к решению xq, удовлетворяющему исходной системе Ax = b. Таким образом, в общем виде алгоритм любого проекционного класса может быть записан следующим образом:

Теорема. Пусть A - симметричная и положительно определенная (SPD, symmetric positive definite) и L = K тогда x˜ - решение проекционного метода на подпространстве K, для любого начального условия x0, если оно минимизирует A-норму ошибки

                                             1∕2
E(˜x) = minx ∈x0+KE (x), E(x) = (A(x* - x),x* - x )
(2.118)

Теорема. Пусть A - произвольная квадратная матрица и L = AK тогда ˜x - решение проекционного метода на подпространстве K ортогональное L, для любого начального условия x0, если оно минимизирует вторую норму невязки

R(˜x) = minx∈x0+KR  (x),  R(x) = ∥(b- Ax )∥2
(2.119)

При построении и реализации проекционных методов важную роль играют так называемые подпространства Крылова, часто выбираемые в качестве K. Подпространством Крылова размерности m, порожденным вектором v и матрицей A называется пространтсво

Km (v,A) = span(v,Av,A2v,...,Am -1v)
(2.120)

В качество v обычно выбирается невязка начального приближения r0 = b - Ax0. Тогда способ построения базисов подпространств и выбор подпространства L полностью определяет вычислительную схему метода.

Для построения базиса в подпространстве Крылова Km(v1,A) можно применить следующий подход. Найдем сначала векторы w1 = v1, w2 = Aw1, w3 = A2w1 = Aw2, ..., wm = Am-1w1 = Awm-1. Тогда

Km (v1,A ) = span(w1,w2,...wm ).
Перейдем теперь от (w1,...,wm) к (v1,...vm), применив процедуру ортогонализации
             ∑k
vk+1 = wk+1 -   αivi
             i=1
(2.121)

и затем пронормируем полученные вектора.

Пусть k векторов построены, тогда (2.121) можно записать как

             k
v    = Av - ∑  α (v ,v )
 k+1     k   i=1  i i j
(2.122)

Для выполнения условия ортогональности vk+1 ко всем предыдущим векторам умножим это равенство скалярно на vj(j k) и приравняем результат к нулю

         ∑k
(Avk,vj) -    αi(vi,vj) = 0
          i=1
(2.123)

Отсюда получаем выражение для коэффициентов αj

αj = (Avk,vj)
Описанный алгоритм называется ортогонализацией Арнольди. Этот метод для построения каждого нового вектора vk требует нахождения (k-1) скалярных произведений и столько же операций линейного комбинирования. Однако, как оказывается, при отказе от требования ортогональности базиса в пользу условия биортогонолизации можно построить процедуру меньшей сложности.

Пусть векторы v1 и w1 таковы, что (v1,w1) = 0 и системы векторов {vi}i=1m и {wi}i=1m определяются соотношениями:

vi+1 = Avi - αivi - βivi- 1, v0 = 0;
       T
wi+1 = A  wi - αiwi - βiwi-1, w0 = 0;
     (Av ,w )
αi = ---i--i;
      (vi,wi)
βi = --(vi,wi)--, β1 = 0.
     (vi-1,wi-1)
(2.124)

Тогда системы {vi}i=1m и {wi}i=1m являются биортогональными и каждая из систем {vi}i=1m и {wi}i=1m является линейно независимой и образует базис вKm(v1,A) и Km(w1,AT) соответственно.

Процедура построения векторов v и w согласно формулам (2.124) называется биортогонализацией Ланцоша.

Пусть пространства K и L связаны соотношением L = AK, причем в качестве K используется подпространство Крылова Km(v1,A) с выбором v1 = r0∕β,β = ||r0||2. Для построения ортонормального базиса K воспользуемся ортогонализацией Арнольди. Полученная вычислительная схема называется методом обобщенных минимальных невязок (GMRES).

А если пространства K и L свяжем соотношением L = K и применим процедуру биортогонализации Ланцоша получим метод бисопряженных градиентов (BiCG). В случае симметричной матрицы получаем метод сопряженных градиентов (CG). Метод бисопряженных градиентов исторически появился позднее, как обобщение CG на несимметричный случай, что и отразилось на его названии, так как вместо ортогонализации в нем используется биортогонализация.