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

数值分析实验 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
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
% 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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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
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
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
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
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
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
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
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
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, -1, 0, 1, 2, 0, 3, 1, -1, 2]
    • 方程组 2

      1
      [3271/27, -9948/71, 4790/161, -4331/72, 2357/216, -1447/54, 293/54, -109/54]
    • 方程组 3

      1
      [2, 1, -3, 0, 1, -2, 3, 0, 1, -1]
  • 系数矩阵的 LU 分解

    • 方程组 1

      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
      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

      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      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

      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
      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 消元法均使用向量化编程,故其运行效率较高