Двухслойные итерационные схемы

Пусть решаем следующую СЛАУ:

\[ A x=f, \quad A \in \mathbb{R}^{n,n}, \quad x,f \in \mathbb{R}^{n,1} \]

До этого решали её прямыми методами - т.е. находили алгоритмами, которые за заранее определенное количество шагов находил решение (тот же метод Гаусса).

Сейчас же хотим решать систему итерационно - т.е. постепенно уточнять наше решение в зависимости от шага итерации. Останавливаемся при достижении определённой точности \(\varepsilon\).

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

В общем виде двухслойная итерационная схема выглядит следующим образом:

\[ R_k \frac{x^{(k+1)}-x^{(k)}}{\tau^{(k)}}+A x^{(k)}=f \]

где \(R_k \in \mathbb{R}^{n,n}, \tau^{(k)} \in \mathbb{R}\) - ненулевые параметры итерации. Матрицы \(R_k\) должны быть невырожденные и положительно определённые.

Очевидно, что если последовательность \(x^{(k)}\) является сходящейся, то она сходится к решению исходной системы. При этом на каждом шаге мы всё равно обязаны решить систему относительно \(x^{(k+1)}\) - это должно быть проще решения изначальной задачи.

Таким образом, коэффициенты итерации \(R_k , \tau^{(k)}\) надо подбирать исходя из:

1. Сходимости последовательности \(x^{(k)}\)

2. Простоты решения задачи \(R_k \cdot x^{(k+1)}=\tau^{(k)}f + \left(R_k \cdot x^{(k)} - A \tau^{(k)} \right) x^{(k)}\)

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

Приведём некоторые примеры двухслойных итераций.

Метод простых итераций

Рассмотрим систему \(Ax = f\). Пусть систему получилось привести к виду

\(x = Rx + b\).

В таком случае можно устроить итерационный процесс :

\(x^{(t+1)} = Rx^{(t)} + b\)

Если последовательность \(\{x^{(s)}\}\) сходится, то МПИ сходится.

Теорема : Если \(||R||\leq q < 1\), то МПИ сходится.

Приведение системы к МПИ:

\(Ax = f\)

\(\tau Ax = \tau f\)

\(x + \tau Ax = x + \tau f\)

\(x = (E - \tau A)x + \tau f\)

Параметр \(\tau\) выбирается из условия сходимости.

Оценка скорости сходимости

Пусть для решения на итерации \(k\) выполняется:

\(||x^{(k)} - x^*||\leq\varepsilon\)

Тогда:

\(||x^{(k)} - x^*||\leq q^k||\delta^0||\leq\varepsilon\)

\(q^k \leq \frac{\varepsilon}{||\delta^0||}\)

Откуда \(k\geq \frac{\ln\varepsilon - \ln||\delta^0||}{\ln q}\)

Теорема о сходимости МПИ.

Для сходимости МПИ необходимо и достаточно \(\max\limits_{i}|\lambda_i(R)|< 1\), где \(\rho(R) = \max\limits_{i}|\lambda_i(R)|\) - спектральный радиус матрицы.

Для доказательства будем пользоваться рядом Неймана.

Ряд Неймана - ряд вида \(\sum\limits_{n = 0}^{\infty} M^n\). Нас интересует его свойство:

• Если спектральный радиус меньше единицы, то ряд Неймана сходится к матрице \((E - A)^{-1} = \sum\limits_{n = 0}^{\infty} A^n\)

Докажем достаточность.

Пусть \(\rho(R) < 1\), тогда \(R^k \rightarrow 0\) при \(k \rightarrow \infty\) и ряд Неймана сходится. Применим последовательно формулу МПИ:

\(x^{(1)} = Rx^{(0)} + b\)

\(x^{(2)} = Rx^{(1)} + b = R^2x^{(0)} + (E + R)b\)

\(x^{(3)} = Rx^{(2)} + b = R^3x^{(0)} + (E + R + R^2)b\)

\(x^{(k+1)} = Rx^{(k)} + b = R^{k + 1}x^{(0)} + (E + R + R^2 + ... + R^{k})b\)

При \(k\rightarrow\infty\) матрица \(R^{k+1}\rightarrow 0 \), а ряд

\(E + R + R^2 + ... + R^{k} = (E - R)^{-1}\).

Исходную систему можно представить в виде:

\((E - R)x = b\)

Подставляя сюда \(x^{*}\):

\((E - R)(E - R)^{-1}b = b\)

Отсюда видно, что подобранный таким образом вектор является решением.

Метод Якоби

Пусть имеем систему:

\[ A \vec{x}=\vec{b} \]

Хотим привести его к итерационному виду:

\[ \vec{x}=B \vec{x}+\vec{g} \]

Метод Якоби основан на отделении диагонали матрицы \(A\) от других её элементов:

\[ A=(L+D+U) \]

где мы разделили нашу матрицу на сумму нижнетреугольной, диагональной и верхнетреугольной. Подставляя в изначальное уравнение:

\[ D\vec{x}=-(L+U) \vec{x}+\vec{b} \]

\[ \vec{x}=-D^{-1}(L+U) \vec{x}+D^{-1}\vec{b} \]

Получили итерационный вид. Осталось навесить номера итерации и

Тогда итерационная схема выглядит как:

\[ \vec{x}^{(k+1)}=B \vec{x}^{(k)}+\vec{g} \]

\[ B=-D^{-1}(L+U), \quad \vec{g} = D^{-1}\vec{b} \]

Или, в раскрытом виде:

\[ x_i^{(k+1)}=\frac{1}{a_{i i}}\left(b_i-\sum_{j \neq i} a_{i j} x_j^{(k)}\right), \quad i=1,2, \ldots, n \]

Теорема 1. Критерий сходимости метода Якоби

Метод Якоби сходится тогда и только тогда, когда \(\left|\lambda_i\right|<1\), где \(\lambda_i\) - корни следующего уравнения:

\[\begin{split} \operatorname{det}\left|\begin{array}{cccc} \lambda a_{11} & a_{12} & \ldots & a_{1 n} \\ a_{21} & \lambda a_{22} & \ldots & a_{2 n} \\ \ldots & \ldots & \ldots & \ldots \\ a_{n 1} & a_{n 2} & \ldots & \lambda a_{n n} \end{array}\right|=0 \end{split}\]

Теорема 2. Достаточное условие сходимости метода Якоби. Пусть \(\|B\|<1\). Тогда при любом выборе начального приближения \(\vec{x}^{(0)}\) :

1. Метод сходится

2. Скорость сходимости метода равна скорости сходимости геометрической прогрессии со знаменателем \(q=\|B\|\)

3. Верна оценка погрешности: \(\left\|\vec{x}^{(k)}-\vec{x}\right\|<q^k\left\|\vec{x}^{(0)}-\vec{x}\right\|\)

Теорема 3. Достаточное условие сходимости метода Якоби. Метод Якоби сходится, если для матрицы изначального уравнения выполнено условие диагонального преобладания:

\[ \forall i \to \left|a_{i i}\right| \geqslant \sum_{j \neq i}\left|a_{i j}\right| \]

Притом хотя бы одно из неравенств должно быть строгим.

Метод Зейделя

Пусть имеем систему:

\[ A \vec{x}=\vec{b} \]

Хотим привести его к итерационному виду:

\[ \vec{x}=B \vec{x}+\vec{g} \]

Опять разобъём матрицу \(A\) на сумму нижнетреугольной, диагональной и верхнетреугольной:

\[ A=(L+D+U) \]

И построим итерационную формулу из изначальной системы:

\[ (L+D) \vec{x}^{(k+1)}=-U \vec{x}^{(k)}+\vec{b} \]

\[ \vec{x}^{(k+1)}=-(L+D)^{-1}U \vec{x}^{(k)}+(L+D)^{-1}\vec{b} \]

Метод Гаусса — Зейделя можно рассматривать как модификацию метода Якоби. Основная идея модификации состоит в том, что новые значения \(\vec{x}^{(i)}\) используются здесь сразу же по мере получения, в то время как в методе Якоби они не используются до следующей итерации:

\[\begin{split} \left\{\begin{array}{l} x_1{ }^{(k+1)}=&c_{12} x_2^{(k)}&+&c_{13} x_3^{(k)}&+&\ldots&+&c_{1 n} x_n{ }^{(k)}&+&d_1 \\ x_2{ }^{(k+1)}=&c_{21} x_1^{(k+1)}&+&c_{23} x_3^{(k)}&+&\ldots&+&c_{2 n} x_n{ }^{(k)}&+&d_2 \\ \ldots \\ x_n{ }^{(k+1)}=&c_{n 1} x_1^{(k+1)}&+&c_{n 2} x_2^{(k+1)}&+&\ldots&+&c_{n(n-1)} x_{n-1}^{(k+1)}&+& d_n \end{array}\right. \end{split}\]

где

\[\begin{split} c_{i j}=\left\{\begin{array}{ll} -\frac{a_{i j}}{a_{i i}}, & j \neq i, \\ 0, & j=i . \end{array} \quad d_i=\frac{b_i}{a_{i i}}, \quad i=1, \ldots, n\right. \end{split}\]

Теорема 4. Критерий сходимости метода Зейделя

Метод Зейделя сходится тогда и только тогда, когда \(\left|\lambda_i\right|<1\), где \(\lambda_i\) - корни следующего уравнения:

\[\begin{split} \operatorname{det}\left|\begin{array}{cccc} \lambda a_{11} & a_{12} & \ldots & a_{1 n} \\ \lambda a_{21} & \lambda a_{22} & \ldots & a_{2 n} \\ \ldots & \ldots & \ldots & \ldots \\ \lambda a_{n 1} & \lambda a_{n 2} & \ldots & \lambda a_{n n} \end{array}\right|=0 \end{split}\]

Теорема 5. Достаточное условие сходимости метода Зейделя. Метод Зейделя сходится, если для матрицы изначального уравнения выполнено условие диагонального преобладания:

\[ \forall i \to \left|a_{i i}\right| \geqslant \sum_{j \neq i}\left|a_{i j}\right| \]

Притом хотя бы одно из неравенств должно быть строгим.

Методы неполной релаксации

В целом, метод Зейделя работает лучше метода Якоби. Это наталкивает на мысль о том, что можно ещё чутка улучшить метод, особо не меняя сути. Помните параметр \(\tau\) в общем определении двухслойной схемы? Сейчас он нам понадобится.

Пусть \(z_i^{(k + 1)}\) - \(i\)-я компонента \(k + 1\) приближения Зейделя. Пусть теперь следующую итерацию мы будем считать по формуле

\[x_i^{(k + 1)} = x_i^{(k)} + \omega (z_i^{(k + 1)} - x_i^{(k)})\]

При \(\omega = 1\) получается метод Зейделя в чистом виде. Метод при \(\omega \in(0, 1)\) называют методом нижней релаксации. Если \(\omega \in(1, 2)\) - методом верхней релаксации. Переписав выражение, видим, что по сути мы «чутка оставляем» прошлый слой при построении следующего:

\[x_i^{(k + 1)} = (1 - \omega)x_i^{(k)} + \omega z_i^{(k +1 )}\]

Матричная запись данного метода выглядит как:

\[(L + \omega D)\frac{x^{(k + 1)} - x^{(k)}}{\omega} + Ax^{(k)} = f\]

Откуда понятно, что параметр \(\tau=\omega\) имеет релаксационный смысл (т.е. своей рода «консервативности» метода).

Теорема 6.: Пусть матрица \(A = A^* > 0\), тогда метод релаксации сходится \(\forall\omega \in(0, 2)\).

При каких \(\omega\) метод сходится быстрее? В общем случае на этот вопрос нет универсального ответа, однако математиками было рассмотрено много частных случаев.

Рассмотрим один важный класс задач, для которого найдена такая \(\omega\).

Пусть \(P\) - матрица смены базиса такая, что : \(PAP^\top = \begin{bmatrix} D_1 & T_{12} \\ T_{21} & D_2 \end{bmatrix}\), где \(D_i\) - диагональные матрицы. Т.е. «разделили переменные на два кластера».

Для таких задач найден оптимальный параметр релаксации

\[\omega_{\text(опт)} = \frac{2}{1 + \sqrt{(1 - \rho^2(B))}}\]

Тут \(B\) - матрица перехода Якоби итерационного метода (\(\rho(B)\) - спектральный радиус, то есть наибольшее по модулю собственное число) в итерационном уравнении вида:

\[ \begin{gathered} \vec{x}=B \vec{x}+\vec{g} \end{gathered} \]

Сравнение двухслойных итерационных методов

Сравним скорости сходимости для конкретной системы:

\[\begin{split} \left\{\begin{array}{l} 2 x_{1}+x_{2}=1 \\ x_{1}+2 x_{2}=-1 \end{array}\right. \end{split}\]

Точное решение системы имеет вид:

\[\begin{split} \vec{x}^{*}=\left(\begin{array}{c} 1 \\ -1 \end{array}\right) \end{split}\]

Решим систему уравнений методами Якоби, Зейделя и ПВР. Тогда для метода Якоби получим:

\[\begin{split} \left\{\begin{array}{l} x_{1}^{(s+1)}=-0.5 x_{2}^{(s)}+0.5 \\ x_{2}^{(s+1)}=-0.5 x_{1}^{(s)}-0.5 \end{array}\right. \end{split}\]

Для метода Зейделя:

\[\begin{split} \left\{\begin{array}{l} x_{1}^{(s+1)}=-0.5 x_{2}^{(s)}+0.5 \\ x_{2}^{(s+1)}=-0.5 x_{1}^{(s+1)}-0.5 \end{array}\right. \end{split}\]

Для метода ПВР:

\[\begin{split} \left\{\begin{array}{l} x_{1}^{(s+1)}=(1-\omega) x_{1}^{(s)}+\omega\left(-0.5 x_{2}^{(s)}+0.5\right) \\ x_{2}^{(s+1)}=(1-\omega) x_{2}^{(s)}+\omega\left(-0.5 x_{1}^{(s+1)}-0.5\right) \end{array}\right. \end{split}\]

где значение параметра было выбрано \(\omega=1.1\).

Для каждого метода рассмотрим три итерации. Тогда получим следующие значения:

\[\begin{split} \begin{array}{|c|c|c|c|} \hline & x_1^{(3)} & x_2^{(3)} & \left\|\vec{x}^{(3)}-\vec{x}^*\right\| \\ \hline \text { Якоби } & 0.875 & -0.875 & 0.125 \\ \hline \text { Зейдель } & 0.969 & -0.984 & 0.031 \\ \hline \text { ПВР } & 1.0008 & -1.0009 & <0.0001 \\ \hline \end{array} \end{split}\]

Из таблицы видно, что метод ПВР сходится быстрее, чем методы Якоби и Зейделя. Вспомним, что для системы, в которой переменные можно разделить на два класса так, чтобы каждое уравнение связывало элемент одного класса с элементами другого, можно найти оптимальное значение параметра \(\omega_{\text {опт }}\). Для систем уравнений с двумя переменными такое разделение выполняется всегда. Таким образом, при условии сходимости метода Якоби для системы второго порядка всегда можно найти оптимальное значение параметра.

Рассмотрим, чему равен спектральный радиус матрицы перехода для метода Якоби.

\[\begin{split} \begin{gathered} \vec{x}=R \vec{x}+\vec{F} \\ R_{\text {Якоби }}=\left(\begin{array}{cc} 0 & -0.5 \\ -0.5 & 0 \end{array}\right) \\ \lambda\left(R_{\text {Якоби }}\right)=\pm \frac{1}{2}, \quad \Rightarrow \quad \rho\left(R_{\text {Якоби }}\right)=\frac{1}{2} \end{gathered} \end{split}\]

Тогда для метода ПВР получим:

\[ \omega_{\text {опт }}=\frac{2}{1+\sqrt{1-\left(\frac{1}{2}\right)^{2}}} \approx 1.07 \]

Именно поэтому значение \(\omega\) было выбрано равным \(1.1 .\)