2021-03-01-数值分析-Day04-Jacobi-GaussSeidel

3.1 - 6 迭代法

基本思想

对线性方程组 $A x = b$ ，当 $A$ 为低阶稠密矩阵时，选主元Gauss消去/Doolittle分解均为有效方法。
工程应用中 $A$ 常为高阶稀疏矩阵，适用迭代法求解：利好计算机存储与性能。

通法：

将 $A x = b$ 改写为 $x = B_{0} x + f$
任取初始值，如 $x^{(0)} = (0, 0, 0)^{T}$ ，将其带入得到方程组解 $x^{(1)} = B_{0} x^{(0)} + f$
依次得： $x^{(2)} = B_{0} x^{(1)} + f, x^{(2)} = B_{0} x^{(1)} + f, \dots, x^{(k)} = B_{0} x^{(k - 1)} + f, \dots$
即得向量序列 $x^{(0)}, x^{(1)}, \dots, x^{(k)}$ ，迭代公式 $x^{(k)} = B_{0} x^{(k - 1)} + f$
迭代次数较高时，向量序列 $x^{(k)}$ 有可能收敛至逼近精确解的值（不一定）。

根据 $x = B x + f$ 变形方式的不同，存在多种迭代算法。

定义1：
对于给定方程组 $x = B x + f$ ，用公式 $x^{k + 1} = B x^{(k)} + f, (k = 0, 1, 2, 3, \dots)$ 逐步代入求近似解的方法称为迭代法（或称为一阶定常迭代法，这里 $B$ 与 $k$ 无关）。
如果 $lim_{k \to \infty} x^{(k)} = x^{*}$ ，即向量序列收敛至精确解序列，则称此迭代法收敛，否则称发散。
由于需要研究 ${x^{(k)}}$ 的收敛性，引进误差向量 $ε^{(k + 1)} = x^{(k + 1)} - x^{*}$ 。易得到： $ε^{(k + 1)} = B ε^{(k)}$ ，递推得： $ε^{(k)} = B ε^{(k - 1)} = \dots = B^{k} ε^{(0)}$ 。

Jacobi迭代

通法细节：
将 $A$ 分裂成 $A = M - N$ ，其中 $M$ 为可选择的非奇异矩阵，使 $M x = d$ 易解，一般选择为 $A$ 的某种近似，称 $M$ 为分裂矩阵。
则原方程 $A x = b$ 转化为 $M x = N x + b$ ，即求解： $x = M^{- 1} N x + M^{- 1} b$ 。
可构造一阶定常迭代法： ${\begin{cases} x^{(0)} \\ x^{(k + 1)} = B x^{(k)} + f (k = 0, 1, \dots) \end{cases} B = M^{- 1} N = M^{- 1} (M - A) = I - M^{- 1} A, f = M^{- 1} b$ 称 $B = I - M^{- 1} A$ 为迭代法的迭代矩阵，选取 $M$ 阵就得到 $A x = b$ 的各种迭代法。

Jacobbi迭代：
设 $a_{i i} \neq 0 (i = 1, 2, \dots, n)$ ，并将 $A$ 写成三部分： $\begin{aligned} A & = (\begin{array}{c} a_{11} \\ a_{22} \\ ⋱ \\ a_{n n} \end{array}) \\ - (\begin{array}{c} 0 \\ - a_{21} & 0 \\ ⋮ & ⋮ & ⋱ \\ - a_{n - 1, 1} & - a_{n - 1, 2} & \dots & 0 \\ - a_{n, 1} & - a_{n, 2} & \dots & - a_{n, n - 1} & 0 \end{array}) \\ - (\begin{array}{c} 0 & - a_{12} & \dots & - a_{1, n - 1} & - a_{1, n} \\ 0 & \dots & - a_{2, n - 1} & - a_{2, n} \\ ⋱ & ⋮ & ⋮ \\ 0 & - a_{n - 1, n} \\ 0 \end{array}) \\ \equiv D - L - U \end{aligned}$ 由此，当 $a_{i i} \neq 0 (i = 1, 2, \dots, n)$ 时，选取 $M$ 为 $A$ 的对角元素部分，即选取 $M = D$ ， $A = D - N$ ，即得到 $A x = b$ 的Jacobi迭代法。 $\begin{aligned} {\begin{matrix} x^{(0)} \\ x^{(k + 1)} = B x^{(k)} + f (k = 0, 1, \dots) \end{matrix} \\ 其中 \begin{aligned} B = D^{- 1} (L + U) \equiv J \\ f = D^{- 1} b \end{aligned} \end{aligned}$ 称 $J$ 为Jacobi迭代法的迭代矩阵。

Jacobi迭代法的分量计算公式：记 $x^{(k)} = (x_{1}^{(k)}, x_{2}^{(k)}, \dots, x_{n}^{(k)})$ ，则有： $\begin{array}{l} \begin{aligned} D x^{k + 1} & = (L + U) x^{(k)} + b a_{i i} x_{i}^{k + 1} \\ = - \sum_{j = 1}^{i - 1} a_{i j} x_{j}^{(k)} - \sum_{j = i + 1}^{n} a_{i j} x_{j}^{(k)} + b_{i} (i = 1, 2, \dots, n) \end{aligned} \\ B = I - M^{- 1} A, f = M^{- 1} b \end{array}$ 因此，解 $A x = b$ 的Jacobi迭代法的计算公式为： ${\begin{cases} x^{(0)} = (x_{1}^{(0)}, \dots, x_{i}^{(0)}, \dots, x_{n}^{(0)})^{T} \\ x_{i}^{(k + 1)} = (b_{i} - \sum_{j = 1, j \neq i}^{n} a_{i j} x_{j}^{(k)}) / a_{i i} \\ i = 1, 2, \dots, n \\ k = 0, 1, \dots 表示迭代次数 \end{cases}$ 矩阵表示为： $\begin{array}{l} A x = b \Leftrightarrow (D - L - U) x = b \Leftrightarrow D x = (L + U) x + b \\ \begin{aligned} D x & = (\begin{array}{c} a_{11} \\ a_{22} \\ ⋱ \\ a_{n n} \end{array}) (\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \\ x_{n} \end{array}) \\ = (\begin{array}{c} 0 & - a_{12} & \dots & - a_{1, n - 1} & - a_{1, n} \\ - a_{21} & 0 & \dots & - a_{2, n - 1} & - a_{2, n} \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ - a_{n - 1, 1} & - a_{n - 1, 2} & \dots & 0 & - a_{n - 1, n} \\ - a_{n, 1} & - a_{n, 2} & \dots & - a_{n, n - 1} & 0 \end{array}) (\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \\ x_{n} \end{array}) + b \\ = (L + U) x + b \end{aligned} \end{array}$

Gauss-Seidel迭代法

改进Jacobi迭代法：选取分裂矩阵 $M$ 为 $A$ 的下三角部分，即 $M = D - L$ ， $A = M - U$ ，得到： $\begin{array}{l} {\begin{matrix} x^{(0)} 初始向量 \\ x^{(k + 1)} = B x^{(k)} + f (k = 0, 1, \dots) \end{matrix} \\ B = I - (D - L^{- 1}) A = (D - L)^{- 1} U \equiv G, f = (D - L)^{- 1} b \end{array}$ 称 $G$ 为Gauss-Seidel迭代法的迭代矩阵。

分量计算公式： $\begin{array}{l} D x^{(k + 1)} = D x^{(k)} - (L x^{(k + 1)} + U x^{k} - D x^{(k)} + b) \\ {\begin{matrix} x_{1}^{(k + 1)} = - \frac{a_{12}}{a_{11}} x_{2}^{(k)} - \frac{a_{13}}{a_{11}} x_{3}^{(k)} - \dots - \frac{a_{1 n}}{a_{11}} x_{n}^{(k)} + \frac{b_{2}}{a_{22}} \\ x_{2}^{(k + 1)} = - \frac{a_{21}}{a_{22}} x_{2}^{(k + 1)} - \frac{a_{23}}{a_{22}} x_{3}^{(k)} - \dots - \frac{a_{2 n}}{a_{22}} x_{n}^{(k)} + \frac{b_{2}}{a_{22}} \\ x_{3}^{(k + 1)} = - \frac{a_{31}}{a_{33}} x_{2}^{(k + 1)} - \frac{a_{32}}{a_{33}} x_{3}^{(k + 1)} - \dots - \frac{a_{3 n}}{a_{33}} x_{n}^{(k)} + \frac{b_{3}}{a_{33}} \\ \dots \\ x_{n}^{(k + 1)} = - \frac{a_{n 1}}{a_{n n}} x_{2}^{(k + 1)} - \frac{a_{n 2}}{a_{n n}} x_{3}^{(k + 1)} - \dots - \frac{a_{n, n - 1}}{a_{n n}} x_{n}^{(k + 1)} + \frac{b_{n}}{a_{n n}} \end{matrix} \\ 即 x_{i}^{k + 1} = \frac{b_{n} - \sum_{j = 1}^{i - 1} a_{i j} x_{j}^{(k + 1)} - \sum_{j = i + 1}^{n} a_{i j} x_{j}^{(k)}}{a_{i i}} \end{array}$

代码：Jacobi迭代

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

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_);

int main(){
    vector<vector<double>> A_{{10, 3, 1}, {2, -10, 3}, {1, 3, 10}};
    vector<double> b_{14, -5, 14};
    iterate(A_, b_);
}

void iterate(vector<vector<double>> A_, vector<double> b_) {
    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] * x[j];
            next_x[i] = 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;
        if (k > 100) break;
        ++k;
    }
}

代码：Gauss-Seidel迭代

/***区别***/
//Jacobi
for (int j = 0; j < size; ++j)
    if (j != i) tmp -= A[i][j] * x[j];
//Gauss-Seidel
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];
}

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

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_);

int main(){
    vector<vector<double>> A_{{10, 3, 1}, {2, -10, 3}, {1, 3, 10}};
    vector<double> b_{14, -5, 14};
    iterate(A_, b_);
}

void iterate(vector<vector<double>> A_, vector<double> b_) {
    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] = 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;
        if (k > 100) break;
        ++k;
    }
}

3.1 - 6 迭代法¶

基本思想¶

Jacobi迭代¶

Gauss-Seidel迭代法¶

代码：Jacobi迭代¶

代码：Gauss-Seidel迭代¶