数值分析实验 - 线性方程组的直接解法

数值分析实验 5 - 线性方程组的直接解法

目的和意义

给出下列几个不同类型的线性方程组，请用适当算法计算其解

1. 线性方程组

   \[ \begin{bmatrix} 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\\ \end{bmatrix}\begin{bmatrix} x_1\\x_2\\x_3\\x_4\\x_5\\x_6\\x_7\\x_8\\x_9\\x_{10} \end{bmatrix}=\begin{bmatrix} 5\\12\\3\\2\\3\\46\\13\\38\\19\\-21 \end{bmatrix} \]

2. 对称正定阵系数阵线方程组

   \[ \begin{bmatrix} 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\\ \end{bmatrix}\begin{bmatrix} x_1\\x_2\\x_3\\x_4\\x_5\\x_6\\x_7\\x_8 \end{bmatrix}=\begin{bmatrix} 0\\-6\\20\\23\\9\\-22\\-15\\45 \end{bmatrix} \]

3. 三对角形线性方程组

   \[ \begin{bmatrix} 4&-1&&&&&&&&\\ -1&4&-1&&&&&&&\\ &-1&4&-1&&&&&&\\ &&-1&4&-1&&&&&\\ &&&-1&4&-1&&&&\\ &&&&-1&4&-1&&&\\ &&&&&-1&4&-1&&\\ &&&&&&-1&4&-1&\\ &&&&&&&-1&4&-1\\ &&&&&&&&-1&4\\ \end{bmatrix}\begin{bmatrix} x_1\\x_2\\x_3\\x_4\\x_5\\x_6\\x_7\\x_8\\x_9\\x_{10} \end{bmatrix}=\begin{bmatrix} 7\\5\\-13\\2\\6\\-12\\14\\-4\\5\\-5 \end{bmatrix} \]

要求

1. 对上述三个方程组分别利用 Gauss 顺序消去法与 Gauss 列主元消去法；平方根法与改进平方根法；追赶法求解 (选择其一)
2. 应用结构程序设计编出通用程序
3. 比较计算结果，分析数值解误差的原因
4. 尽可能利用相应模块输出系数矩阵的三角分解式

目的和意义

1. 通过该课题的实验，体会模块化结构程序设计方法的优点
2. 运用所学的计算方法，解决各类线性方程组的直接算法
3. 提高分析和解决问题的能力，做到学以致用
4. 通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点

计算公式

Gauss 消元法

任意情况均可使用

过程即将系数矩阵通过初等行变换化为上三角矩阵或下三角矩阵，之后回代即可

平方根法

当系数矩阵为对称正定矩阵时可以使用

设方程组为 \(Ax=b\)

对系数矩阵 \(A\) 做分解

\[ A=LL^T \]

其中 \(L=(l_{ij})\) 为下三角矩阵，

\[ l_{ij}=\begin{cases} \sqrt{a_{ii}-\sum_{k=1}^{i-1}l_{ik}^2},&i=j\\ \displaystyle{a_{ij}-\sum_{k=1}^{j-1}l_{ik}l_{jk}\over l_{jj}},&i\ne j \end{cases} \]

接下来计算

\[ \begin{cases} Ly=b\\ L^Tx=y \end{cases} \]

即可

改进平方根法

当系数矩阵为对称正定矩阵时可以使用

设方程组为 \(Ax=b\)

对系数矩阵 \(A\) 做分解

\[ A=LDL^T \]

其中 \(L=(l_{ij})\) 为下三角矩阵，\(D=(d_{i})\) 为对角矩阵，令 \(t_{ij}=l_{ij}d_j\), 有

\[ t_{ij}=a_{ij}-\sum_{k=1}^{j-1}t_{ik}l_{jk} \]

\[ d_{i}=a_{ii}-\sum_{k=1}^{i-1}t_{ik}l_{ik} \]

接下来计算

\[ \begin{cases} Ly=b\\ L^Tx=D^{-1}y \end{cases} \]

即可

追赶法

当系数矩阵为三对角矩阵时可以使用

设方程组为 \(Ax=d\)

对系数矩阵做分解

\[ A=\begin{bmatrix} b_1&c_1&&&\\ a_2&b_2&c_2&&\\ &\ddots&\ddots&\ddots&\\ &&a_{n-1}&b_{n-1}&c_{n-1}\\ &&&a_n&b_n \end{bmatrix}=\begin{bmatrix} \alpha_1&&&&\\ \gamma_2&\alpha_2&&&\\ &\ddots&\ddots&&\\ &&\gamma_{n-1}&\alpha_{n-1}&\\ &&&\gamma_n&\alpha_n \end{bmatrix}\begin{bmatrix} 1&\beta_1&&&\\ &1&\beta_2&&\\ &&\ddots&\ddots&\\ &&&1&\beta_{n-1}\\ &&&&1 \end{bmatrix} \]

有

\[ \alpha_i=b_i-a_i\beta_{i-1} \]

\[ \beta_i=\frac{c_i}{\alpha_i} \]

\[ \gamma_i=a_i \]

令

\[ L=\begin{bmatrix} \alpha_1&&&&\\ \gamma_2&\alpha_2&&&\\ &\ddots&\ddots&&\\ &&\gamma_{n-1}&\alpha_{n-1}&\\ &&&\gamma_n&\alpha_n \end{bmatrix},~U=\begin{bmatrix} 1&\beta_1&&&\\ &1&\beta_2&&\\ &&\ddots&\ddots&\\ &&&1&\beta_{n-1}\\ &&&&1 \end{bmatrix} \]

接下来计算

\[ \begin{cases} Ly=d\\ Ux=y \end{cases} \]

即可

程序设计

主程序

Show code

main.mview raw

% Exp.5

% @Author: Tifa
% @LastEditTime: 2021-06-09 00:22:44

% 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 = 1;
now_method = @gauss_advanced;

now_method(A{now_equ}, b{now_equ})

% LU matrix factorization
[L, U] = lu(sym(A{now_equ}))

输入数据检查

Show code

input_check.mview raw

function input_check(A, b)
    % @Author: Tifa
    % @LastEditTime: 2021-06-09 00:22:44

    if ~iscolumn(b)
        error('Error: b should be a column vector!')
    end
    if ~ismatrix(A)
        error('Error: A should be a matrix!')
    end
    if size(A, 1) ~= size(A, 2)
        error('Error: A should be a square matrix!')
    end
    if det(A) == 0
        error('Error: A should be a nonsingular matrix!')
    end
end

Gauss 顺序消去法

Show code

gauss_basic.mview raw

function x = gauss_basic(A, b)

    % @Author: Tifa
    % @LastEditTime: 2021-06-09 00:22:44

    input_check(A, b)

    A = [A b];
    len = length(b);

    for i = 1:len - 1
        A(i + 1:len, i + 1:len + 1) = A(i + 1:len, i + 1:len + 1) - A(i + 1:len, i) * A(i, i + 1:len + 1) / A(i, i);
    end

    b = A(:, end);
    x = zeros(len, 1);
    for i = len:-1:1
        x(i) = (b(i) - A(i, i + 1:len) * x(i + 1:len)) / A(i, i);
    end
end

Gauss 列主元消去法

Show code

gauss_advanced.mview raw

function x = gauss_advanced(A, b)

    % @Author: Tifa
    % @LastEditTime: 2021-06-09 00:22:44

    input_check(A, b)

    A = [A b];
    len = length(b);

    for i = 1:len - 1
        [~, idxA] = max(abs(A(i:len, i)));
        idxA = idxA + i - 1;
        tmp = A(idxA, :);
        A(idxA, :) = A(i, :);
        A(i, :) = tmp;
        A(i + 1:len, i + 1:len + 1) = A(i + 1:len, i + 1:len + 1) - A(i + 1:len, i) * A(i, i + 1:len + 1) / A(i, i);
    end

    b = A(:, end);
    x = zeros(len, 1);
    for i = len:-1:1
        x(i) = (b(i) - A(i, i + 1:len) * x(i + 1:len)) / A(i, i);
    end
end

平方根法

Show code

sqrt_basic.mview raw

function x = sqrt_basic(A, b)

    % @Author: Tifa
    % @LastEditTime: 2021-06-09 00:22:44

    input_check(A, b)

    len = length(b);
    for k = 1:len
        if det(A(1:k, 1:k)) <= 0
            error('A should be a positive definite matrix')
        end
    end
    if ~issymmetric(A)
        error('A should be a symmetric matrix')
    end

    for k = 1:len
        A(k, k) = sqrt(A(k, k));
        A(k + 1:len, k) = A(k + 1:len, k) / A(k, k);
        for j = k + 1:len
            A(j:len, j) = A(j:len, j) - A(j:len, j) * A(j, k);
        end
    end
    for j = 1:len - 1
        b(j) = b(j) / A(j, j);
        b(j + 1:len) = b(j + 1:len) - b(j) * A(j + 1:len, j);
    end
    for j = len:-1:2
        b(j) = b(j) / A(j, j);
        b(1:j - 1) = b(1:j - 1) - b(j) * A(1:j - 1, j);
    end
    x = b;
end

改进平方根法

Show code

sqrt_advanced.mview raw

function x = sqrt_advanced(A, b)

    % @Author: Tifa
    % @LastEditTime: 2021-06-09 00:22:44

    input_check(A, b)

    len = length(b);
    for k = 1:len
        if det(A(1:k, 1:k)) <= 0
            error('A should be a positive definite matrix')
        end
    end
    if ~issymmetric(A)
        error('A should be a symmetric matrix')
    end

    L = diag(ones(len, 1));
    D = zeros(len); D(1, 1) = A(1, 1);
    S = zeros(len);

    for i = 2:len
        for j = 1:i - 1
            S(i, j) = A(i, j) - S(i, 1:j - 1) * L(j, 1:j - 1)';
        end
        L(i, 1:i - 1) = S(i, 1:i - 1) / D(1:i - 1, 1:i - 1)';
        D(i, i) = A(i, i) - S(i, 1:i - 1) * L(i, 1:i - 1)';
    end

    x = zeros(len, 1);
    y = zeros(len, 1);
    for i = 1:len
        y(i) = (b(i) - L(i, 1:i - 1) * D(1:i - 1, 1:i - 1) * y(1:i - 1)) / D(i, i);
    end
    x(len) = y(len);
    for i = len - 1:-1:1
        x(i) = y(i) - L(i + 1:len, i)' * x(i + 1:len);
    end
end

追赶法

Show code

chase.mview raw

function x = chase(A, f)

    % @Author: Tifa
    % @LastEditTime: 2021-06-09 00:22:44

    input_check(A, f);

    if ~isempty(find(A ~= diag(diag(A, -1), -1) + diag(diag(A)) + diag(diag(A, 1), 1), 1))
        error('A should be a tri-diagonal matrix')
    end

    a = [0; diag(A, -1)];
    b = diag(A);
    c = [diag(A, 1); 0];

    len = length(f);
    u = zeros(len, 1);
    d = zeros(len, 1); d(1) = b(1);
    for i = 1:len - 1
        u(i) = c(i) / d(i);
        d(i + 1) = b(i + 1) - a(i + 1) * u(i);
    end

    y = zeros(len, 1); y(1) = f(1) / d(1);
    for i = 2:len
        y(i) = (f(i) - a(i) * y(i - 1)) / d(i);
    end
    x = zeros(len, 1); x(len) = y(len);
    for i = len - 1:-1:1
        x(i) = y(i) - u(i) * x(i + 1);
    end
end

结果讨论和分析

结果

• 解

  • 方程组 1

    [1, -1, 0, 1, 2, 0, 3, 1, -1, 2]

  • 方程组 2

    [3271/27, -9948/71, 4790/161, -4331/72, 2357/216, -1447/54, 293/54, -109/54]

  • 方程组 3

    [2, 1, -3, 0, 1, -2, 3, 0, 1, -1]

• 系数矩阵的 LU 分解

  • 方程组 1

    L =
    
    [  1,  0,  0,  0,  0,     0,   0,   0,        0, 0]
    [  2,  1,  0,  0,  0,     0,   0,   0,        0, 0]
    [  1,  0,  1,  0,  0,     0,   0,   0,        0, 0]
    [  0, -1,  2,  1,  0,     0,   0,   0,        0, 0]
    [ -1,  2,  1,  0,  1,     0,   0,   0,        0, 0]
    [  2,  1, -3,  2,  2,     1,   0,   0,        0, 0]
    [  0,  1, -2,  1, -1,  -2/5,   1,   0,        0, 0]
    [  4,  1,  0, -1,  2,  29/5,   3,   1,        0, 0]
    [  1,  2,  3, -1,  1,   2/5,  19,  -6,        1, 0]
    [  0,  0, -1,  2,  0, -31/5, -12, 3/2, -131/250, 1]
    
    U =
    
    [ 4, 2, -3, -1,  2, 1,   0,   0,    0,       0]
    [ 0, 2,  1, -1,  2, 3,   0,   1,    0,       0]
    [ 0, 0,  1,  0,  1, 1,  -1,   0,    3,       1]
    [ 0, 0,  0,  4, -1, 4,   1,   2,    3,       2]
    [ 0, 0,  0,  0,  3, 1,  -2,   1,   -1,       2]
    [ 0, 0,  0,  0,  0, 5,   1,  -1,    2,       0]
    [ 0, 0,  0,  0,  0, 0, 2/5, 3/5, 14/5,       3]
    [ 0, 0,  0,  0,  0, 0,   0,   2,  -13,      -9]
    [ 0, 0,  0,  0,  0, 0,   0,   0, -125,    -110]
    [ 0, 0,  0,  0,  0, 0,   0,   0,    0, -607/50]

  • 方程组 2

    L =
    
    [   1,  0,    0,  0,    0,  0,   0, 0]
    [ 1/2,  1,    0,  0,    0,  0,   0, 0]
    [  -1,  1,    1,  0,    0,  0,   0, 0]
    [   0, -2,  1/3,  1,    0,  0,   0, 0]
    [ 1/2,  0, -2/3,  1,    1,  0,   0, 0]
    [   1,  1,    0, -2,  1/4,  1,   0, 0]
    [   0,  2,  1/3,  0, -1/2,  1,   1, 0]
    [   0,  0,  2/3,  1,    0, -2, 1/2, 1]
    
    U =
    
    [ 4, 2, -4,  0,  2,  4,  0,  0]
    [ 0, 1,  1, -2,  0,  1,  2,  0]
    [ 0, 0,  9,  3, -6,  0,  3,  6]
    [ 0, 0,  0,  1,  1, -2,  0,  1]
    [ 0, 0,  0,  0, 16,  4, -8,  0]
    [ 0, 0,  0,  0,  0,  1,  1, -2]
    [ 0, 0,  0,  0,  0,  0,  4,  2]
    [ 0, 0,  0,  0,  0,  0,  0,  9]

  • 方程组 3

    L =
    
    [    1,     0,      0,       0,        0,         0,           0,            0,             0, 0]
    [ -1/4,     1,      0,       0,        0,         0,           0,            0,             0, 0]
    [    0, -4/15,      1,       0,        0,         0,           0,            0,             0, 0]
    [    0,     0, -15/56,       1,        0,         0,           0,            0,             0, 0]
    [    0,     0,      0, -56/209,        1,         0,           0,            0,             0, 0]
    [    0,     0,      0,       0, -209/780,         1,           0,            0,             0, 0]
    [    0,     0,      0,       0,        0, -780/2911,           1,            0,             0, 0]
    [    0,     0,      0,       0,        0,         0, -2911/10864,            1,             0, 0]
    [    0,     0,      0,       0,        0,         0,           0, -10864/40545,             1, 0]
    [    0,     0,      0,       0,        0,         0,           0,            0, -40545/151316, 1]
    
    U =
    
    [ 4,   -1,     0,      0,       0,        0,          0,           0,            0,             0]
    [ 0, 15/4,    -1,      0,       0,        0,          0,           0,            0,             0]
    [ 0,    0, 56/15,     -1,       0,        0,          0,           0,            0,             0]
    [ 0,    0,     0, 209/56,      -1,        0,          0,           0,            0,             0]
    [ 0,    0,     0,      0, 780/209,       -1,          0,           0,            0,             0]
    [ 0,    0,     0,      0,       0, 2911/780,         -1,           0,            0,             0]
    [ 0,    0,     0,      0,       0,        0, 10864/2911,          -1,            0,             0]
    [ 0,    0,     0,      0,       0,        0,          0, 40545/10864,           -1,             0]
    [ 0,    0,     0,      0,       0,        0,          0,           0, 151316/40545,            -1]
    [ 0,    0,     0,      0,       0,        0,          0,           0,            0, 564719/151316]

• 耗时 (s)

  	Gauss 顺序消去法	Gauss 列主元消去法	平方根法	改进平方根法	追赶法
  1	0.012	0.005	-	-	-
  2	0.004	0.006	0.015	0.013	-
  3	0.004	0.005	0.008	0.007	0.004

分析

• Gauss 消元法的适用范围最广
• 平方根法只适用于对称正定矩阵
• 追赶法只适用于三对角矩阵
• 由于两种 Gauss 消元法均使用向量化编程，故其运行效率较高