为了和教材上列主元高斯消元法的求解过程一致，编写了MATLAB程序，供大家参考，互相交流学习。
先写出函数
function x = gsxylzy(A,b)
% format rat;clc
% A = [-3 2 6;10 -7 0; 5 -1 5];
% b = [4; 7 ;6];
M = [A,b]
n = length(b);
%%% 消元过程
for i = 1:n-1
%交换列主元
for k = i:n-1
ak = max(abs(M(k:n,k))); index = find(M(:,k)==ak);
if isempty(index)
index = find(M(:,k)==-ak);
end
temp = M(index,:); M(index,:) = M(k,:); M(k,:) = temp;
end
fprintf(["第 ",num2str(i)," 次交换列主元n"])
M
% 消元
for k = i:n-1
M(k+1,:) = M(k+1,:)-M(i,:)*M(k+1,i)/M(i,i);
end
fprintf(["第 ",num2str(i)," 列消元n"])
M
end
%%% 回代过程
x = zeros(n,1);
x(n) = M(n,n+1)/M(n,n);
fprintf(["第 ",num2str(n)," 个变量 x(",num2str(n),")=%sn"],num2str(x(n)));
for k = (n-1):-1:1
x(k) = (M(k,n+1)-M(k,k+1:n)*x(k+1:n))/M(k,k);
fprintf(["第 ",num2str(k)," 个变量 x(",num2str(k),")=% sn"],num2str(x(k)));
end
主程序测试应用：
format rat;clc
A = [-3 2 6;10 -7 0; 5 -1 5];
b = [4; 7 ;6];
x = gsxylzy(A,b);
fprintf("方程组的解为n")
x
运行结果为
M =
-3 2 6 4
10 -7 0 7
5 -1 5 6
第 1 次交换列主元
M =
10 -7 0 7
-3 2 6 4
5 -1 5 6
第 1 列消元
M =
10 -7 0 7
0 -1/10 6 61/10
0 5/2 5 5/2
第 2 次交换列主元
M =
10 -7 0 7
0 5/2 5 5/2
0 -1/10 6 61/10
第 2 列消元
M =
10 -7 0 7
0 5/2 5 5/2
0 0 31/5 31/5
第 3 个变量 x(3)=1
第 2 个变量 x(2)=-1
第 1 个变量 x(1)=2.6645e-16
方程组的解为
x =
1/3752999689475414 %这里严格说是0
-1
1
>>
如果不想出现这种分数近似，也可以用符号形式，则不会出现近似值.
format rat;clc
A = sym([-3 2 6;10 -7 0; 5 -1 5]);
b = sym([4; 7 ;6]);
x = gsxylzy(A,b);
fprintf("方程组的解为n")
x
则运行结果为
M =
[ -3, 2, 6, 4]
[ 10, -7, 0, 7]
[ 5, -1, 5, 6]
第 1 次交换列主元
M =
[ 10, -7, 0, 7]
[ -3, 2, 6, 4]
[ 5, -1, 5, 6]
第 1 列消元
M =
[ 10, -7, 0, 7]
[ 0, -1/10, 6, 61/10]
[ 0, 5/2, 5, 5/2]
第 2 次交换列主元
M =
[ 10, -7, 0, 7]
[ 0, 5/2, 5, 5/2]
[ 0, -1/10, 6, 61/10]
第 2 列消元
M =
[ 10, -7, 0, 7]
[ 0, 5/2, 5, 5/2]
[ 0, 0, 31/5, 31/5]
第 3 个变量 x(3)=1
第 2 个变量 x(2)=-1
第 1 个变量 x(1)=0
方程组的解为
x =
0
-1
1
>>
这里使用了符号变量命令 sym，结果和教材完全一样。
