数值分析实验 - 解线性方程组的迭代法

发表于 更新于 分类于 ， ， ， 本文字数： 872 阅读时长 ≈ 3 分钟

数值分析实验 6 - 解线性方程组的迭代法

目的和意义

问题提出

对 数值分析实验 - 线性方程组的直接解法 所列目的和意义的线性方程组，试分别选用 Jacobi 迭代法，Gauss-Seidol 迭代法和 SOR 方法计算其解

要求

  1. 体会迭代法求解线性方程组，并能与消去法做以比较；
  2. 分别对不同精度要求，如 \(\epsilon=10^{-3},10^{-4},10^{-5}\), 由迭代次数体会该迭代法的收敛快慢；
  3. 对方程组 2, 3 使用 SOR 方法时，选取松弛因子 \(\omega=0.8, 0.9, 1, 1.1, 1.2\) 等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者；
  4. 给出各种算法的设计程序和计算结果

目的和意义

  1. 通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较；
  2. 运用所学的迭代法算法，解决各类线性方程组，编出算法程序；
  3. 体会上机计算时，终止步骤 \(\|x^{(k+1)}-x^{(k)}\|_{\infty}<\epsilon\) 或 \(k>\) (予给的迭代次数), 对迭代法敛散性的意义；
  4. 体会初始解 \(x^{(0)}\), 松弛因子的选取，对计算结果的影响

计算公式

令

\[ A=D-L-U \]

其中

  • \[ D=\begin{bmatrix} a_{11}\\ &a_{22}\\ \phantom{\vdots}&&\ddots\\ &&&a_{nn} \end{bmatrix} \]
  • \[ L=\begin{bmatrix} 0\\ -a_{21}&0\\ \vdots&\vdots&\ddots\\ -a_{n1}&-a_{n2}&\dots&0 \end{bmatrix} \]
  • \[ U=\begin{bmatrix} 0&-a_{12}&\dots&-a_{1n}\\ &0&\dots&-a_{2n}\\ &&\ddots&\vdots\\ &&&0 \end{bmatrix} \]

Jacobi 迭代法

令

\[ G=D^{-1}(L+U) \]

则

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

收敛条件为 \(\rho(G)<1\)

Gauss-Seidol 迭代法

令

\[ G=(D-L)^{-1}U \]

则

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

收敛条件为 \(\rho(G)<1\)

SOR 方法

令

\[ G=(D-\omega L)^{-1}((1-\omega)D+\omega U) \]

则

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

收敛条件为 \(\rho(G)<1\)

程序设计

主程序

Show code

main.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
% Exp.6

% @Author: Tifa
% @LastEditTime: 2021-06-22 23:03:26

% Data
A = {
[
    4 2 -3 -1 2 1 0 0 0 0;
    8 6 -5 -3 6 5 0 1 0 0;
    4 2 -2 -1 3 2 -1 0 3 1;
    0 -2 1 5 -1 3 -1 1 9 4;
    -4 2 6 -1 6 7 -3 3 2 3;
    8 6 -8 5 7 17 2 6 -3 5;
    0 2 -1 3 -4 2 5 3 0 1;
    16 10 -11 -9 17 34 2 -1 2 2;
    4 6 2 -7 13 9 2 0 12 4;
    0 0 -1 8 -3 -24 -8 6 3 -1;
    ], ...
    [
    4 2 -4 0 2 4 0 0;
    2 2 -1 -2 1 3 2 0;
    -4 -1 14 1 -8 -3 5 6;
    0 -2 1 6 -1 -4 -3 3;
    2 1 -8 -1 22 4 -10 -3;
    4 3 -3 -4 4 11 1 -4;
    0 2 5 -3 -10 1 14 2;
    0 0 6 3 -3 -4 2 19;
    ], ...
        diag(ones(10, 1) * 4) + diag(-ones(9, 1), 1) + diag(-ones(9, 1), -1)
    };

b = {
    [5 12 3 2 3 46 13 38 19 -21]'
    [0 -6 20 23 9 -22 -15 45]'
    [7 5 -13 2 6 -12 14 -4 5 -5]'
    };

now_equ = 2;
now_method = @jacobi;
now_epsi = 1e-4;

if isequal(now_method, @sor)
    omegas=0.8:0.05:1.9;
    [x, k] = sor(A{now_equ}, b{now_equ}, now_epsi, omegas);
    x
    plot(omegas, k)
    grid on
    xlabel('{\omega}')
    ylabel('Steps')
else
    [x, k] = now_method(A{now_equ}, b{now_equ}, now_epsi)
end

Jacobi 迭代法

Show code

jacobi.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function [x, k] = jacobi(A, b, epsi)

    % @Author: Tifa
    % @LastEditTime: 2021-06-22 23:03:26

    G = diag(diag(A)) \ -(tril(A, -1) + triu(A, 1));
    if max(abs(eig(G))) >= 1
        error('Not convergent!');
    end
    k = 1;
    x_pre = ones(length(b), 1);
    b = diag(diag(A)) \ b;
    x = G * x_pre + b;
    while norm(x - x_pre, inf) > epsi
        x_pre = x;
        x = G * x + b;
        k = k + 1;
    end
end

Gauss-Seidol 迭代法

Show code

gauss_seidol.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function [x, k] = gauss_seidol(A, b, epsi)

    % @Author: Tifa
    % @LastEditTime: 2021-06-22 23:03:26

    G = (diag(diag(A)) + tril(A, -1)) \ -triu(A, 1);
    if max(abs(eig(G))) >= 1
      error('Not convergent!');
    end
    k = 1;
    x_pre = ones(length(b), 1);
    b = (diag(diag(A)) + tril(A, -1)) \ b;
    x = G * x_pre + b;
    while norm(x - x_pre, inf) > epsi
        x_pre = x;
        x = G * x + b;
        k = k + 1;
    end
end

SOR 方法

Show code

sor.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
function [x, k] = sor(A, b, epsi, omegas)

    % @Author: Tifa
    % @LastEditTime: 2021-06-22 23:03:26

    x = zeros(length(b), length(omegas));
    k = zeros(1, length(omegas));
    for i = 1:length(omegas)
        G = (diag(diag(A)) + omegas(i) * tril(A, -1)) \ ((1 - omegas(i)) * diag(diag(A)) - omegas(i) * triu(A, 1));
        if max(abs(eig(G))) >= 1
            error('Not convergent!');
        end
        k(i) = 1;
        x_pre = ones(length(b), 1);
        b1 = omegas(i) * (diag(diag(A)) + omegas(i) * tril(A, -1)) \ b;
        xx = G * x_pre + b1;
        while norm(xx - x_pre, inf) > epsi
            x_pre = xx;
            xx = G * xx + b1;
            k(i) = k(i) + 1;
        end
        x(:, i) = xx;
    end
end

结果讨论和分析

结果

方程组的解经验证符合要求，为节省篇幅而略去

  • \(\epsilon=10^{-3}\) 时的迭代次数

    方程组 1方程组 2方程组 3
    Jacobi 迭代法不收敛不收敛15
    Gauss-Seidol 迭代法不收敛12567
    SOR 方法 (\(\omega=0.8\))不收敛189510
    SOR 方法 (\(\omega=1.2\))不收敛8429

  • \(\epsilon=10^{-4}\) 时的迭代次数

    方程组 1方程组 2方程组 3
    Jacobi 迭代法不收敛不收敛15
    Gauss-Seidol 迭代法不收敛16929
    SOR 方法 (\(\omega=0.8\))不收敛255113
    SOR 方法 (\(\omega=1.2\))不收敛113211

  • \(\epsilon=10^{-5}\) 时的迭代次数

    方程组 1方程组 2方程组 3
    Jacobi 迭代法不收敛不收敛15
    Gauss-Seidol 迭代法不收敛212911
    SOR 方法 (\(\omega=0.8\))不收敛320716
    SOR 方法 (\(\omega=1.2\))不收敛142313

  • SOR 方法中，方程组 2 的迭代次数与 \(\omega\) 的变化图像

  • SOR 方法中，方程组 3 的迭代次数与 \(\omega\) 的变化图像

结论

  • SOR 方法中
    • 对于方程组 2, \(\omega=1.9\) 时的迭代步数最小，且与 \(\omega=0.8\) 时的迭代步数差别较大
    • 对于方程组 3, \(\omega=1.0\) 时的迭代步数最小，不同 \(\omega\) 之间迭代步数差别不大