贡献者: addis
这里给出一个简单的程序演示高斯消元法的基本步骤。注意在实际应用中,我们解线性方程组一般使用 Matlab 的反斜杠算符:x=A \ y
,其中 $ \boldsymbol{\mathbf{A}} $ 是系数矩阵,$ \boldsymbol{\mathbf{y}} $ 是常数列,$ \boldsymbol{\mathbf{x}} $ 是未知量。
% 高斯消元法得到梯形系数矩阵
% 显示每个行变换的操作和结果
% 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 中的增广矩阵测试程序如下:
>> 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