2021-03-02-数值分析-Day05-SOR-迭代法收敛性

续：3.1 - 6 迭代法

逐次超松弛迭代法（SOR迭代法）

选取分裂矩阵 \(M\) 为带参数的下三角阵： \(M=\displaystyle\frac{1}{\omega}(D-\omega L),\ B=I-M^{-1}A,\ f=M^{-1}b\) ，其中 \(w>0\) 为可选择的松弛因子。

构造迭代法，迭代矩阵为： \(L_\omega=I-\omega(D-\omega L)^{-1}A=(D-\omega L)^{-1}((1-\omega)D+\omega U)\) 。

则解 \(Ax=b\) 的SOR方法即为： \(\begin{array}{l}\left\{\begin{array}{l}x^{(0)}\\x^{(k+1)}=L_\omega x^{(k)}+f\quad(k=0,1,\dots)\end{array}\right.\\其中\begin{array}{l}L_\omega=(D-\omega L)^{-1}((1-\omega)D+\omega U),\\f=\omega(D-\omega L)^{-1}b\end{array}\end{array}\) 。

推导得： \(\begin{array}{l}(D-\omega L)x^{(k+1)}=((1-\omega)D+\omega U)x^{(k)}+\omega b\\或Dx^{(k+1)}=Dx^{(k+1)}+\omega(b+Lx^{(k+1)}+Ux^{(k)}-Dx^{(k)})\end{array}\) 。

分量计算公式为： \(\displaystyle x_i^{(k+1)}=x_1^{(k)}+\omega(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum_{j=i}^{n}a_{ij}x_j^{(k)})/a_{ii}\) 。

可令 \(\Delta x_i=\omega(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum_{j=i}^{n}a_{ij}x_j^{(k)})/a_{ii}\) ，则 \(x_i^{(k+1)}=x_i^{(k)}+\Delta x_i\) 。

Gauss-Seidel迭代法是SOR迭代法的一个特殊形式（系数 \(\omega=1\) ）。

当 \(\omega<1\) 时，称为超松弛法；当 \(\omega>1\) 时，称为低松弛法。

计算机中，常用 \(\max\limits_{1\le i\le n}|\Delta x_i|=\max\limits_{1\le i\le n}\left|x_i^{(k+1)}-x_i^{(k)}\right|<\varepsilon\) 或者 \(\left\|r^{(k)}\right\|=\left\|b-Ax^{(k)}\right\|\) 作为迭代终止条件。

迭代法的收敛性

设 \(Ax=b\) ，其中 \(A\in R^{n\times n}\) 为非奇异矩阵，记 \(x^*\) 为原方程组精确解，且设有等价的方程组： \(Ax=b\Leftrightarrow x=Bx+f\) ，则 \(x^*=Bx^*+f\) 。设有一阶定常迭代法 \(x^{(k+1)}=Bx^{(k)}\) 。引进误差向量 \(\varepsilon^{(k)}=x^{(k)}\) ，得到误差向量递推公式 \(\varepsilon^{(k+1)}=B\varepsilon^{(k)}\ \Rightarrow\ \varepsilon^{(k)}=B^k\varepsilon^{(0)}\) 。则研究问题从 \(\varepsilon^{(k)}\to0\) 转换为 \(B^k\to0\) 。

定理1： \(\lim\limits_{x\to\infty}A_k=A\Leftrightarrow\lim\limits_{x\to\infty}\left\|A_k-A\right\|,\ \|\|为任一种算子范数\)

定理2： \(\lim\limits_{x\to\infty}A_k=A\Leftrightarrow\lim\limits_{x\to\infty}A_kx=Ax,\ x\in R^n\)

定理3： \(\lim\limits_{x\to\infty}A_k=0\Leftrightarrow\rho(A)<1,\ \rho\ 谱半径\)

迭代法收敛基本定理：对方程组 \(x=Bx+f\) ，及对应一阶定常迭代法 \(x^{(k+1)}=Bx^{(k)}+f\) ，迭代法收敛的充要条件为矩阵 \(B\) 的谱半径 \(\rho(B)<1\) 。

由于 \(\rho(B)\le\|B\|\) ，因此矩阵 \(B\) 的范数也可用于判别迭代法的收敛性：
迭代法收敛充分条件： \(\exists \|B\|=q<1\) ，迭代法收敛，且有 \(\begin{array}{l}\left\|x^*-x^{(k)}\right\|\le q^k\left\|x^*-x^{(0)}\right\|\\\left\|x^*-x^{(k)}\right\|\le\displaystyle\frac{q}{1-q}\left\|x^k-x^{(k-1)}\right\|\\\left\|x^*-x^{(k)}\right\|\le\displaystyle\frac{q}{1-q^k}\left\|x^1-x^0\right\|\end{array}\) 。根据第二式也可得到迭代法的一种终止条件参考量： \(\left\|x^k-x^{(k-1)}\right\|\) 。

几种特殊迭代法的收敛性

如果矩阵 \(A\) 的元素满足： \(|a_{ii}|>\displaystyle\sum_{j=1,j\neq i}^n|a_{ij}|,\ (i=1,2,\dots,n)\) ，则称 \(A\) 为严格对角占优阵。

若 \(|a_{ii}|\ge\displaystyle\sum_{j=1,j\neq i}^n|a_{ij}|,\ (i=1,2,\dots,n)\) ，且至少有一个不等式严格成立，则称 \(A\) 为弱对角占优阵。

对角占优定理：若矩阵 \(A\) 为严格对角占优矩阵，则 \(A\) 为非奇异矩阵。

Jacobi和Gauss迭代法收敛的充分条件：对方程 \(Ax=b\) ，如果 \(A\) 为严格对角占优阵，则解该方程的Jacobi迭代法和Gauss-Seidel迭代法均收敛。

SOR迭代法收敛的必要条件： \(0<\omega<2\) 。

SOR迭代法收敛的充分条件1：对方程 \(Ax=b\) ，如果 \(A\) 为对称正定矩阵， \(0<\omega<2\) ，则解该方程的SOR迭代法收敛。

SOR迭代法收敛的充分条件2：对方程 \(Ax=b\) ，如果 \(A\) 为严格对角占优矩阵（或弱对角占优不可约矩阵）， \(0<\omega\le1\) ，则解该方程的SOR迭代法收敛。

代码：SOR迭代与两种终止条件

//
// Created by xa on 2021-03-01.
//
#include <iostream>
#include <vector>
#include <algorithm>

using std::vector;

vector<vector<double>> A;
vector<double> b;
vector<double> x;
vector<double> next_x;

int size = 0;

void iterate(vector<vector<double>> A_, vector<double> b_);
void iterate(vector<vector<double>> A_, vector<double> b_, double omega);

int main(){
    vector<vector<double>> A_{{4, -2, -4}, {-2, 17, 10}, {-4, 10, 9}};
    vector<double> b_{10, 3, -7};
    iterate(A_, b_, 1.46);
}

void iterate(vector<vector<double>> A_, vector<double> b_) {iterate(A_,b_,1);}
void iterate(vector<vector<double>> A_, vector<double> b_, double omega) {
    if (omega <= 0 || omega >= 1) { std::cout << "SOR does not converge !" << std::endl; return; }

    A.assign(A_.begin(), A_.end());
    b.assign(b_.begin(), b_.end());
    size = A.size();

    vector<double> x(size);
    vector<double> next_x(size);

    int k = 0;
    while (true) {
        for (int i = 0; i < size; ++i) {
            double tmp = b[i];
            for (int j = 0; j < size; ++j) {
                if (j < i) tmp -= A[i][j] * next_x[j];
                if (j >= i) tmp -= A[i][j] * x[j];
            }
            next_x[i] = x[i] + omega * tmp / A[i][i];
        }

        for (int i = 0; i < size; ++i) {
            x[i] = next_x[i];
            std::cout << x[i] << ' ';
        }
        std::cout << " 第"<< k << "次" << std::endl;

        double r = 0;
        for (int i = 0; i < size; ++i) {
            double tmp = b[i];
            for (int j = 0; j < size; ++j)
                tmp -= A[i][j] * x[j];
            r += tmp * tmp;
        }
        if (r < 1e-10) break;
        ++k;
    }
}

void iterate(vector<vector<double>> A_, vector<double> b_, double omega) {
    A.assign(A_.begin(), A_.end());
    b.assign(b_.begin(), b_.end());
    size = A.size();

    vector<double> x(size);
    vector<double> next_x(size);

    int k = 0;
    while (true) {
        for (int i = 0; i < size; ++i) {
            double tmp = b[i];
            for (int j = 0; j < size; ++j) {
                if (j < i) tmp -= A[i][j] * next_x[j];
                if (j >= i) tmp -= A[i][j] * x[j];
            }
            next_x[i] = x[i] + omega * tmp / A[i][i];
        }

        double maxDelta = std::abs(*(std::max_element(std::begin(next_x), std::end(next_x))) -
                                   *(std::max_element(std::begin(x), std::end(x))));
        std::cout << maxDelta << " | ";

        for (int i = 0; i < size; ++i) {
            x[i] = next_x[i];
            std::cout << x[i] << ' ';
        }
        std::cout << " 第"<< k << "次" << std::endl;

        if (maxDelta < 1e-3) break;
        ++k;
    }
}