4. Ітераційні методи для систем

4.1. Ітераційні методи розв’язання СЛАР

Систему

\begin{equation} \label{eq:4.1.1} A \vec x = \vec b \end{equation}

зводимо до вигляду

\begin{equation} \label{eq:4.1.2} \vec x = B \vec x + \vec f. \end{equation}

Будь яка система

\begin{equation} \label{eq:4.1.3} \vec x = \vec x - C \cdot (A \vec x - \vec b) \end{equation}

має вигляд \eqref{eq:4.1.2} і при $\det C \ne 0$ еквівалентна системі \eqref{eq:4.1.1}. Наприклад, для $C = \tau \cdot E$:

\begin{equation} \vec x = \vec x - \tau \cdot (A \vec x - \vec b). \tag{3’} \end{equation}

4.1.1. Метод простої ітерації

Цей метод застосовується до рівняння \eqref{eq:4.1.2}

\begin{equation} \label{eq:4.1.4} \vec x^{(k + 1)} = B \vec x^{(k)} + \vec f, \end{equation}

де $\vec x^{(0)}$ — початкове наближення, задано.

Теорема: Ітераційний процес збігається, тобто

\begin{equation} \left| \vec x^{(k)} - \vec x\right| \xrightarrow[k\to\infty]{} 0, \end{equation}

якщо

\begin{equation} \label{eq:4.1.5} |B| \le q < 1. \end{equation}

При цьому має місце оцінка

\begin{equation} \label{eq:4.1.6} \left|\vec x^{(n)} - \vec x\right| \le \dfrac{q^n}{1-q}\cdot\left|\vec x^{(1)} - \vec x^{(0)}\right|. \end{equation}

4.1.2. Метод Якобі

Припустимо $\forall i$: $a_{i,i} \ne 0$. Зведемо систему \eqref{eq:4.1.1} до вигляду

\begin{equation} x_i = -\Sum_{j = 1}^{i - 1} \dfrac{a_{i,j}}{a_{i,i}} \cdot x_j - \Sum_{j = i + 1}^n \dfrac{a_{i,j}}{a_{i,i}} \cdot x_j + \dfrac{b_i}{a_{i,i}}, \end{equation}

де $i = \overline{1, n}$.

Ітераційний процес запишемо у вигляді

\begin{equation} \label{eq:4.1.7} x_i^{(k+1)} = -\Sum_{j=1}^{i-1} \dfrac{a_{i,j}}{a_{i,i}} \cdot x_j^{(k)} - \Sum_{j=i+1}^n \dfrac{a_{i,j}}{a_{i,i}} \cdot x_j^{(k)} + \dfrac{b_i}{a_{i,i}}, \end{equation}

де $k = 0,1,\ldots$, а $i=\overline{1,n}$.

Теорема: Ітераційний процес збігається до розв’язку, якщо виконується умова

\begin{equation} \forall i: \Sum_{\substack{j = 1 \\ i \ne j}}^n |a_{i,j}| \le |a_{i,i}|. \end{equation}

Це умова діагональної переваги матриці $A$.

Теорема: Якщо ж

\begin{equation} \label{eq:4.1.8} \forall i: \Sum_{\substack{j = 1 \\ i \ne j}}^n |a_{i,j}| \le q\cdot|a_{i,i}|, \quad 0 \le q < 1. \end{equation}

то має місце оцінка точності:

\begin{equation} \left|\vec x^{(n)} - \vec x\right| \le \dfrac{q^n}{1-q} \cdot \left|\vec x^{(0)}-\vec x\right|. \end{equation}

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

В компонентному вигляді ітераційний метод Зейделя записується так:

\begin{equation} \label{eq:4.1.9} x_i^{(k+1)} = -\Sum_{j=1}^{i-1} \dfrac{a_{i,j}}{a_{i,i}} \cdot x_j^{(k+1)} - \Sum_{j=i+1}^n \dfrac{a_{i,j}}{a_{i,i}} \cdot x_j^{(k)} + \dfrac{b_i}{a_{i,i}}, \end{equation}

де $k = 0,1,\ldots$, а $i=\overline{1,n}$.

На відміну від методу Якобі на $k$-му-кроці попередні компоненти розв’язку беруться з $(k+1)$-ої ітерації.

Теорема: Достатня умова збіжності методу Зейделя — $A^\intercal = A > 0$.

4.1.4. Матрична інтерпретація методів Якобі і Зейделя

Подамо матрицю $A$ у вигляді

\begin{equation} A = A_1 + D + A_2, \end{equation}

де $A_1$ — нижній трикутник матриці $A$, $A_2$ — верхній трикутник матриці $A$, $D$ — її діагональ. Тоді систему \eqref{eq:4.1.1} запишемо у вигляді

\begin{equation} D \vec x = A_1 \vec x + A_2 \vec x + \vec b, \end{equation}

або

\begin{equation} \vec x = D^{-1} A_1 \vec x + D^{-1} A_2 \vec x + D^{-1} \vec b, \end{equation}

Матричний запис методу Якобі:

\begin{equation} \vec x^{(k+1)} = D^{-1} A_1 \vec x^{(k)} + D^{-1} A_2 \vec x^{(k)} + D^{-1} \vec b, \end{equation}

методу Зейделя:

\begin{equation} \vec x^{(k+1)} = D^{-1} A_1 \vec x^{(k+1)} + D^{-1} A_2 \vec x^{(k)} + D^{-1} \vec b, \end{equation}

Теорема: Необхідна і достатня умова збіжності методу Якобі: всі корені рівняння

\begin{equation} \det(D + \lambda(A_1 + A_2 )) = 0 \end{equation}

по модулю більше 1.

Теорема: Необхідна і достатня умова збіжності методу Зейделя: всі корені рівняння

\begin{equation} \det(A_1 + D + \lambda A_2) = 0 \end{equation}

по модулю більше 1.

4.1.5. Однокрокові (двошарові) ітераційні методи

Канонічною формою однокрокового ітераційного методу розв’язку СЛАР є його запис у вигляді

\begin{equation} \label{eq:4.1.10} B_k \dfrac{\vec x^{(k+1)} - \vec x^{(k)}}{\tau_{k+1}} + A \vec x^{(k)} = \vec b, \end{equation}

Тут $\{B_k\}$ — послідовність матриць (пере-обумовлюючі матриці), що задають ітераційний метод на кожному кроці; $\{\tau_{k+1}\}$ — ітераційні параметри.

Означення: Якщо $B_k = E$, то ітераційний процес називається явним

\begin{equation} \vec x^{(k+1)} = \vec x^{(k)} - \tau_{k+1} \left(A \vec x^{(k)} + \vec b\right). \end{equation}

Означення: Якщо $B_k \ne E$, то ітераційний процес називається неявним

\begin{equation} B_k \vec x^{(k+1)} = F^k. \end{equation}

У цьому випадку на кожній ітерації необхідно розв’язувати СЛАР.

Означення: Якщо $\tau_{k+1} \equiv \tau$, $B_k \equiv B$, то ітераційний процес називається стаціонарним; інакше — нестаціонарним.

Методам, що розглянуті вище відповідають:

• методу простої ітерації: $B_k = E$, $\tau_{k+1} = \tau$;

• методу Якобі: $B_k = D$, $\tau_{k+1} = 1$;

• методу Зейделя: $B_k = D + A_1$, $\tau_{k+1} = 1$.

4.1.6. Збіжності стаціонарних ітераційних процесів у випадку симетричних матриць

Розглянемо випадок симетричних матриць $A^\intercal=A$ і стаціонарний ітераційний процес $B_k \equiv E$, $\tau_{k+1} \equiv \tau$.

Нехай для $A$ справедливі нерівності

\begin{equation} \label{eq:4.1.11} \gamma_1 E \le A \le \gamma_2 E, \quad \gamma_1, \gamma_2 > 0. \end{equation}

Тоді при виборі $\tau = \tau_0 = \frac{2}{\gamma_1 + \gamma_2}$ ітераційний процес збігається. Найбільш точним значенням $\gamma_1$, $\gamma_2$ при яких виконуються обмеження \eqref{eq:4.1.11} є $\gamma_1 = \min \lambda_i(A)$, $\gamma_2 = \max \lambda_i(A)$. Тоді

\begin{equation} q = q_0 = \dfrac{\gamma_2 - \gamma_1}{\gamma_2 + \gamma_1} = \dfrac{1-\xi}{1+\xi}, \quad \xi = \dfrac{\gamma_1}{\gamma_2}. \end{equation}

і справедлива оцінка

\begin{equation} \left|\vec x^{(n)} - \vec x\right| \le \dfrac{q^n}{1-q} \cdot \left|\vec x^{(0)} - \vec x\right|. \end{equation}

Зауваження: аналогічно обчислюється $q$ і для методу релаксації розв’язання нелінійних рівнянь, де $\gamma_1 = m = \min \vert f'(x)\vert$, $\gamma_2 = M_1 = \max\vert f'(x)\vert$.

Явний метод з багатьма параметрами $\{\tau_k\}$:

\begin{equation} B \equiv E, \quad {\tau_k}: \Min_\tau q(\tau), \quad n=n(\varepsilon)\to\min, \end{equation}

які обчислюються за допомогою нулів багаточлена Чебишова, називаються ітераційним методом з чебишевським набором параметрів.

4.1.7. Метод верхньої релаксації

Узагальненням методу Зейделя є метод верхньої релаксації:

\begin{equation} (D + \omega A_1) \cdot\dfrac{\vec x^{(k+1)} + \vec x^{(k)}}{\omega} + A \vec x^{(k)} = \vec b, \end{equation}

де $D$ — діагональна матриця з елементами $a_{i,i}$ по діагоналі. $\omega > 0$ — заданий числовий параметр.

Тепер $B = D + \omega A_1$, $\tau = \omega$. Якщо $A^\intercal = A > 0$, то метод верхньої релаксації збігається при умові $0 < \omega < 2$. Параметр підбирається експериментально з умови мінімальної кількості ітерацій.

4.1.8. Методи варіаційного типу

До цих методів відносяться: метод мінімальних нев’язок, метод мінімальних поправок, метод найшвидшого спуску, метод спряжених градієнтів. Вони дозволяють обчислювати наближення без використання апріорної інформації про $\gamma_1$, $\gamma_2$ в \eqref{eq:4.1.11}.

Нехай $B = E$. Для методу мінімальних нев’язок параметри $\tau_{k+1}$ обчислюються з умови

\begin{equation} \left|\vec r^{(k+1)}\right|^2 = \left|\vec r^{(k)}\right|^2 - 2\tau_{k+1}\cdot\left\langle\vec r^{(k)}, A\vec r^{(k)}\right\rangle + \tau_{k+1}^2 \cdot\left|A\vec r^{(k)}\right|^2 \to \min. \end{equation}

Тому

\begin{equation} \tau_{k+1} = \dfrac{\left\langle A\vec r^{(k)}, \vec r^{(k)}\right\rangle }{\left|\vec r^{(k)}\right|^2}, \end{equation}

де $\vec r^{(k)} = A \vec x^{(k)} - \vec b$ — нев’язка.

Умова для завершення ітераційного процесу:

\begin{equation} \left|\vec r^{(n)}\right| < \varepsilon. \end{equation}

Швидкість збіжності цього методу співпадає із швидкістю методу простої ітерації з одним оптимальним параметром $\tau_0 = \frac{2}{\gamma_1+\gamma_2}$.

Аналогічно будуються методи з $B \ne E$. Матриця $B$ називається переобумовлювачем і дозволяє підвищити швидкість збіжності ітераційного процесу. Його вибирають з умов

• легко розв’язувати СЛАР $B \vec x^{(k)} = F_k$ (діагональний, трикутній, добуток трикутніх та інше);

• зменшення числа обумовленості матриці $B^{-1}A$ у порівнянні з $A$.

4.2. Методи розв’язання нелінійних систем

Розглянемо систему рівнянь

\begin{equation} \left\{ \begin{aligned} & f_1(x_1, \ldots, x_n) = 0, \newline & \ldots \newline & f_n(x_1,\ldots,x_n) = 0. \end{aligned} \right. \end{equation}

Перепишемо її у векторному вигляді:

\begin{equation} \label{eq:4.2.1} \vec f(\vec x) = 0. \end{equation}

4.2.1. Метод простої ітерації

В цьому методі рівняння \eqref{eq:4.2.1} зводиться до еквівалентного вигляду

\begin{equation} \label{eq:4.2.2} \vec x = \vec \Phi(\vec x). \end{equation}

Ітераційний процес представимо у вигляді:

\begin{equation} \label{eq:4.2.3} \vec x^{(k+1)} = \vec \Phi\left(\vec x^{(k)}\right). \end{equation}

початкове наближення $\vec x^{(0)}$ — задано.

Нехай оператор $\vec \Phi$ визначений на множині $H$. За теоремою про стискуючі відображення ітераційний процес \eqref{eq:4.2.3} сходиться, якщо виконується умова

\begin{equation} \label{eq:4.2.4} \left| \vec \Phi(\vec x) - \vec \Phi(\vec y) \right| \le q \cdot |\vec x - \vec y|, \quad 0 < q < 1, \end{equation}

або

\begin{equation} \label{eq:4.2.5} \left| \vec \Phi’(\vec x)\right| \le q < 1, \end{equation}

де $\vec x\in U_r$, $\vec \Phi'(\vec x) = \left(\dfrac{\partial \varphi_i}{\partial x_j}\right)\void_{i,j=1}^n$. Для похибки справедлива оцінка

\begin{equation} \left| \vec x^{(m)} - \vec x\right| \le \dfrac{q^n}{1 - q} \cdot\left|\vec x^{(0)} - \vec x\right|. \end{equation}

Частинним випадком методу простої ітерації є метод релаксації для рівняння \eqref{eq:4.2.1}:

\begin{equation} \vec x^{(k+1)} = \vec x^{(k)} - \tau \cdot \vec F\left(\vec x^{(k)}\right), \end{equation}

де $\tau < 2 / \left\|\vec F'(\vec x)\right\|$.

4.2.2. Метод Ньютона

Розглянемо рівняння

\begin{equation} \vec F(\vec x) = 0. \end{equation}

Представимо його у вигляді

\begin{equation} \label{eq:4.2.6} \vec F\left(\vec x^{(k)}\right) + \vec F’\left(\vec \xi^{(k)}\right)\cdot\left(\vec x - \vec x^{(k)}\right) = 0, \end{equation}

де

\begin{equation} \vec \xi^{(k)} = \vec x^{(k)} + \theta_k \cdot \left(\vec x^{(k)} - \vec x\right), \end{equation}

де $0 < \theta_k < 1$. Тут $\vec F'(\vec x) = \left(\dfrac{\partial f_i}{\partial x_j}\right)\void_{i,j=1}^n$ — матриця Якобі для $\vec F(\vec x)$. Можемо наближено вважати $\vec \xi^{(k)} \approx \vec x^{(k)}$. Тоді з \eqref{eq:4.2.6} матимемо

\begin{equation} \label{eq:4.2.7} \vec F\left(\vec x^{(k)}\right) + \vec F’\left(\vec x^{(k)}\right) \cdot\left(\vec x^{(k+1)} - \vec x^{(k)}\right) = 0. \end{equation}

Ітераційний процес представимо у вигляді:

\begin{equation} \label{eq:4.2.8} \vec x^{(k+1)} = \vec x^{(k)} - \vec F’\left(\vec x^{(k)}\right)^{-1} \cdot \vec F\left(\vec x^{(k)}\right). \end{equation}

Для реалізації методу Ньютона потрібно, щоб існувала обернена матриця

\begin{equation} \vec F’\left(\vec x^{(k)}\right)^{-1}. \end{equation}

Можна не шукати обернену матрицю, а розв’язувати на кожній ітерації СЛАР

\begin{align} A_k \vec z^{(k)} &= \vec F\left(\vec x^{(k)}\right), \label{eq:4.2.9} \newline \vec x^{(k + 1)} &= \vec x^{(k)} - \vec z^{(k)}, \nonumber \end{align}

де $k=0,1,2,\ldots$, і $\vec x^{(0)}$ — задано, а матриця $A_k = \vec F'\left(\vec x^{(k)}\right)$.

Метод має квадратичну збіжність, якщо добре вибрано початкове наближення. Складність методу (при умові використання методу Гаусса розв’язання СЛАР \eqref{eq:4.2.9} на кожній ітерації $Q_n = \frac23 n^3+O(n^2)$, де $n$ — розмірність системи \eqref{eq:4.2.1}.

4.2.3. Модифікований метод Ньютона

Ітераційний процес має вигляд:

\begin{equation} \vec x^{(k+1)} = \vec x^{(k)} - \vec F’\left(\vec x^{(0)}\right)^{-1} \cdot\vec F\left(\vec x^{(k)}\right). \end{equation}

Тепер обернена матриця обчислюється тільки на нульовій ітерації. На інших — обчислення нового наближення зводиться до множення матриці $A_0 = \vec F'\left(\vec x^{(0)}\right)^{-1}$ на вектор $\vec F\left(\vec x^{(k)}\right)$ та додавання до $\vec x^{(k)}$.

Запишемо метод у вигляді системи лінійних рівнянь (аналог \eqref{eq:4.2.9})

\begin{align} A_0 \vec z^{(k)} &= \vec F\left(\vec x^{(k)}\right), \label{eq:4.2.10} \newline \vec x^{(k + 1)} &= \vec x^{(k)} - \vec z^{(k)}, \nonumber \end{align}

де $k=0,1,2,\ldots$.

Оскільки матиця $A_0$ розкладається на трикутні (або обертається) один раз, то складність цього методу на одній ітерації (окрім нульової) $Q_n = O(n^2)$. Але цей метод має лінійну швидкість збіжності.

Можливе циклічне застосування модифікованого методу Ньютона, тобто коли обернену матрицю похідних шукаємо та обертаємо через певне число кроків ітераційного процесу.

Задача 9: Побудувати аналог методу січних для систем нелінійних рівнянь.