Общий подход к построению проекционных методов

Проекционные методы и подпространства Крылова

Ускорение сходимости релаксационных методов

Скорость сходимости методов, основанных на расщеплении, непосредственно связана со спектральным радиусом матрицы K ^–1 R; с другой стороны, выбор K ограничен требованием легкой обратимости. Одним из распространенных способов улучшения сходимости является введение параметра. Пусть w– некоторое вещественное число. Рассмотрим вместо системы (1.1) масштабированную систему

wAx=wb,

(2.10)

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

wA=(D–wE)–(wF+(1–w)D),

(2.11)

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

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

(D–wE)x k +1=(wF+(1–w)D)x k +wb.

(2.12)

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

K_SOR(w)=D–wE R_SOR(w)=wF+(1–w)D.

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

Выражение (2.11) остается тождеством, если в нем поменять местами матрицы E и F. Такая перестановка дает обратный метод последовательной верхней релаксации:

(D–wF)x _{k +1}=[wE+(1–w)D]x _k +wb.

Последовательное применение прямого и обратного методов SOR дает симметричный метод последовательной верхней релаксации (SSOR)

(D–wE)x _{k +1/2}=[wF+(1–w)D]x _k +wb; (D–wF)x _{k +1}=[wE+(1–w)D]x _{k +1/2}+wb.

Рассмотрим систему (1.1) и сформируем для нее следующую задачу. Пусть заданы некоторые два подпространства K Ì R^N и L Ì R^N. Требуется найти такой вектор x Î K, который обеспечивал бы решение исходной системы, «оптимальное относительно подпространства L», т.е. чтобы выполнялось условие

l Î L: (Ax, l)=(b, l),

называемое условием Петрова-Галеркина. Сгруппировав обе части равенства по свойствам скалярного произведения и заметив, что b – Ax = r _x, это условие можно переписать в виде

l Î L: (rx, l)=0,

(2.13)

т.е. r _x= b – Ax ^ L. Такая задача называется задачей проектирования x на подпространство K ортогонально к подпространству L.

В более общей постановке задача выглядит следующим образом. Для исходной системы (1.1) известно некоторое приближение x ₀ к решению x _*. Требуется уточнить его поправку d _xÎ K таким образом, чтобы b – A (x ₀+ d _x) ^ L. Условие Петрова-Галеркина в этом случае можно записать в виде

l Î L: (rx₀+dx, l)=((b–Ax₀)–Ad_x, l)=(r₀–Ad_x, l)=0.

Пусть dim K =dim L = m (где dim – размер подпространств). Введем в подпространство K и L базисы и соответственно. Нетрудно видеть, что (2.19) выполняется тогда и только тогда, когда

j (1£ j £ m): (r0–Adx,w j)=0.

(2.14)

При введении для базисов матричных обозначений V =| n ₁| n ₂|¼ | n _m | и W =| w ₁| w ₂|¼ | w _m | можно записать d _x= Vy где yÌ R^m – вектор коэффициентов. Тогда (2.14) может быть записано в виде

W ^T (r0–AVy)=0,

(2.15)

откуда W ^T AVy = W ^T r ₀ и

y=(W ^T AV)^–1W ^T r0.

(2.16)

Таким образом, решение должно уточняться в соответствии с формулой

x1=x0+V(W ^T AV)^–1W ^T r0,

(2.17)

из которой сразу вытекает важное требование: в практических реализациях проекционных методов подпространства и их базисы должны выбираться, так чтобы W ^T AV либо была малой размерности, либо имела простую структуру, удобную для обращения.

Из (2.15) также вытекает соотношение

Vy=A^–1r₀=A^–1(b–Ax₀)=x_*–x₀,

т.е. Vy = d _x представляет собой проекцию на подпространство K разности между точным решением и начальным приближением.

Пусть имеется набор пар подпространств таких, что K ₁Å K ₂Å… K_q Å= R^N и L ₁Å L ₂Å… L_q Å= R^N (запись K Å L обозначает симметричную разность множеств K и L, т.е. элемент принадлежит K Å L, если он находится либо в K, либо в L, но не в пересечении K и L, т.е. K Å L = K È L – K Ç L). Тогда очевидно, что последовательное применение описанного ранее процесса ко всем таким парам приведет к решению x _q, удовлетворяющему исходной системе Ax = b. Соответственно, в общем, виде алгоритм любого метода проекционного класса может быть записан следующим образом:

Начало

Выбор пары подпространств K и L

Построение для K и L базисов V=|n₁|n₂|¼ |n _m | и W=|w₁|w₂|¼ |w _m |

r=b–Ax

y=(W ^T AV)^–1W ^T r

x=x+Vy

Продолжать до достижения сходимости

2.2.2 Случай одномерных подпространств K и L

Наиболее простой ситуацией является случай, когда пространства K и L одномерны. Пусть K и L подпространства являющиеся множествами векторов представимых в виде линейной комбинации векторов { n ₁, n ₂,…} и { w ₁, w ₂,…} т.е. K =span{ n } и L =span{ w }. Тогда (2.17) принимает вид

x k +1=x k +g k n k,

(2.18)

причем коэффициент g _k легко находится из условия ортогональности r _k – A (g _k n _k) ^ w _k:

(r _k –Ag _k n _k, w _k)=(r _k, w _k)–g _k (An _k, w _k)=0,

откуда g _k =(r _k, w _k)/(An _k, w _k).

Если выбрать n _k = w _k = r _k, тогда (2.18) принимает вид

x k +1=x k +r k (r k, r k)/(Ar k, r k).

(2.19)

Поскольку выражение в знаменателе представляет собой квадратичную форму r _k^T Ar _k, сходимость процесса гарантирована, если матрица A является симметричной и положительно определенной. Данный метод называется методом наискорейшего спуска. В практических задачах метод наискорейшего спуска обладает достаточно медленной сходимостью.

При выборе n _k = A ^T r _k и w _k = An _k формула (2.18) примет вид

x k =1=x k +A ^T r k (r k, AA ^T r k)/(AA ^T r k, AA ^T r k)= =x k +A ^T r k (A ^T r k, A ^T r k)/(A ^T AA ^T r k, A ^T r k).

(2.20)

Данный метод называется методом наискорейшего уменьшения невязки; условием его сходимости является невырожденность матрицы A. Сравнивая (2.19) и (2.20), нетрудно убедится, что метод наискорейшего уменьшения невязки совпадает с методом наискорейшего спуска, примененным к системе A ^T Ax = A ^T b.

Если положить K = L и на различных итерациях в качестве вектора n _k циклически с повторением выбирать единичные орты e ₁, e ₂,…, e _N, e ₁,…то получится рассмотренный ранее метод Гаусса-Зейделя (обратный порядок выбора соответствует обратному методу Гаусса-Зейделя).

Выбор на k -й итерации n _k = w _k = A ^T e _k дает, так называемый, ABS-класс методов. Чтобы для различных итераций выполнилось условие K_i ^ K_j, возможно либо хранение всех n _k с их ортогонализацией по мере нахождения (схема Хуанга), либо пересчет матрицы A (схема Абрамова). Первый вариант ведет к увеличению расходов памяти для реализации алгоритма, второй – к изменению заполненности матрицы. Следовательно, данные методы пригодны лишь для решения небольших систем; с другой стороны, их единственным условием сходимости является существование решения как такового, а потому данный класс пригоден для СЛАУ с вырожденными и неквадратными матрицами. Непосредственно из определения векторов n _k следует, что в данных методах на k -й итерации поправка к решению вычисляется из условия обращения k -го уравнения в тождество. Если предполагать, что матрица A квадратная и невырожденная, вычислительная схема имеет следующий вид (вариант с пересчетом матрицы, метод ABR1ORT):

Алгоритм метода ABR1ORT

Выполнять для i =1,…, N

Выполнять для j = i +1,…, N

a _j* =a _j* –a×a _i*

b_j = b_j –a× b_i

увеличить j

увеличить i

Пример использования метода ABR1ORT в системе TALGAT:

- действительный случай: SET "solve" ABR1ORT_r real_m right_r;

- комплексный случай: SET "solve" ABR1ORT complex_m right_c;

Дата добавления: 2014-12-27