高斯消元法程序

预备知识　线性方程组高斯消元法

　　这里给出一个简单的程序演示高斯消元法的基本步骤．注意在实际应用中，我们解线性方程组一般使用 Matlab 的反斜杠算符：x=A \ y ，其中 $ \boldsymbol{\mathbf{A}} $ 是系数矩阵，$ \boldsymbol{\mathbf{y}} $ 是常数列，$ \boldsymbol{\mathbf{x}} $ 是未知量．

GaussEli.m

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
54
55
56
57
58
59
60
61
62
63

% 高斯消元法得到梯形系数矩阵
% 显示每个行变换的操作和结果
% A 为任意尺寸的矩阵
% 如果 A 是增广矩阵， 输入 augmented = true， 如果 A 是系数矩阵， augmented = false
% 输出的 A 是梯形矩阵， q 是一个行矢量， q(ii) 是第 ii 行第一个不为零的列标
% 如果第 ii 行全为零， 则 q(ii) = 0

function [A, q] = GaussEli(A, augmented)
% 系数矩阵的行数 Nr 和列数 Nc
[Nr, Nc] = size(A);
if (augmented)
	Nc = Nc - 1;
end
q = zeros(1,Nr); q(1) = 1;
disp(A);
for ii = 1:Nr-1
	% 处理第 ii 行
	if (ii > 1)
		q(ii) = q(ii - 1) + 1;
	end
	for jj = q(ii) : Nc
		% 检查第 jj 列
		if (max(abs(A(ii:end, q(ii)))) > 0)
			% 消元使 A(ii+1:end, q(ii)) 全为 0
            A = eli1(A, ii, q(ii));
            break;
		else
			% 不需要消元，检查下一列
            if (q(ii) < Nc)
                q(ii) = q(ii) + 1;
            else
                return;
            end
        end
    end
end
end

% 做行变换使 A(ii+1:end, q) 全为 0
% A(ii:end,q) 不能全为 0
function A = eli1(A, ii, q)
% 交换两行，使 A(ii,q) 在 A(ii:end, q) 中最大
[~, iimax] = max(abs(A(ii:end, q)));
iimax = ii + iimax - 1;
if (iimax > ii)
    disp(['交换: r_', num2str(ii), ' <-> r_', num2str(iimax)]);
    temp = A(ii,:);
    A(ii,:) = A(iimax,:);
    A(iimax,:) = temp;
    disp(A);
end

% 用第 ii 行消去所有 A(ii+1:end, q)
for jj = ii+1 : size(A,1)
    if (A(jj,q) == 0)
        continue;
    end
    k = -A(jj,q)/A(ii,q);
    disp(['消元: ', num2str(k), ' * r_', num2str(ii), ' + r_', num2str(jj)]);
    A(jj,:) = A(ii,:) * k + A(jj,:);
    disp(A);
end
end

　　与 “高斯消元法” 中的步骤略有不同的是，该程序在处理每一行 A(ii,:) 时都会试图做一个行交换使系数 A(ii,q(ii)) 的绝对值尽可能大．这是为了减小数值误差：试想如果 A(ii,q(ii)) 的解析值为 0，但由于数值误差，计算出来是一个很小的数（例如 1e-15），那么用其消元时就可能需要将第 ii 行乘以一个很大的系数（例如 1e15）再加到另一行上，导致程序不稳定．

　　用例 2 中的增广矩阵测试程序如下：

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
54
55
56
57
58

>> A = [1 1 -1 1; 2 2 -2 1; 1 1 0 2; 2 2 -1 5];
>> y = [3; 7; 3; 4];
>> [A1, q] =GaussEli([A,y],true); q
     1     1    -1     1     3
     2     2    -2     1     7
     1     1     0     2     3
     2     2    -1     5     4

交换: r_1 <-> r_2
     2     2    -2     1     7
     1     1    -1     1     3
     1     1     0     2     3
     2     2    -1     5     4

消元: -0.5 * r_1 + r_2
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0         0    0.5000   -0.5000
    1.0000    1.0000         0    2.0000    3.0000
    2.0000    2.0000   -1.0000    5.0000    4.0000

消元: -0.5 * r_1 + r_3
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0         0    0.5000   -0.5000
         0         0    1.0000    1.5000   -0.5000
    2.0000    2.0000   -1.0000    5.0000    4.0000

消元: -1 * r_1 + r_4
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0         0    0.5000   -0.5000
         0         0    1.0000    1.5000   -0.5000
         0         0    1.0000    4.0000   -3.0000

交换: r_2 <-> r_3
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0    1.0000    1.5000   -0.5000
         0         0         0    0.5000   -0.5000
         0         0    1.0000    4.0000   -3.0000

消元: -1 * r_2 + r_4
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0    1.0000    1.5000   -0.5000
         0         0         0    0.5000   -0.5000
         0         0         0    2.5000   -2.5000

交换: r_3 <-> r_4
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0    1.0000    1.5000   -0.5000
         0         0         0    2.5000   -2.5000
         0         0         0    0.5000   -0.5000

消元: -0.2 * r_3 + r_4
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0    1.0000    1.5000   -0.5000
         0         0         0    2.5000   -2.5000
         0         0         0         0         0
         
q =
     1     3     4     0

致读者：小时物理百科一直以来坚持所有内容免费且不做广告，这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化，例如大量广告，内容付费，会员制，甚至被收购。因此，我们鼓起勇气在此请求广大读者热心捐款，使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元，我们几天内就能脱离亏损状态，并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。