[旧]选主元消元法求矩阵的逆－程序示例

XuQiang

今天上机课上，信计2班的王静同学给出了《数值分析引论》P300 习题8的一个Matlab程序，见附件liezhuyuanGJ.m。请同学们参考
一下，并提一些意见。

另外，我们提一个关于此程序的问题。程序中用选主元消元法求解矩阵A的逆矩阵时，是使用每一个对角线A(k,k)元素所在的行去消所有行中的A(i,k),其中1<= i <=n。请问这种消元法有什么优点和缺点？
% liezhuyuanGJ.m
% 王静 07级信计二班 0237
% P300 习题8
%A为非奇异阵时，用列主元高斯-约当法求A的逆
%[A,B]用初等行变换将A变成E，对应的B变成A逆
clear; clc;
A = [2,1,-3,-1;3,1,0,7;-1,2,4,-2;1,0,-1,5];
n = length(A);
B = eye(n,n);
C = [A,B];
%d用于C的换行
d = zeros(1,2*n);
for k=1:n
max = 0;
maxrow = k;
%求从A(k,k)到A(n,k)中绝对值最大的元素
for i = k:n
%dangqian=|A(i,k)|
dangqian = A(i,k);
if A(i,k) < 0
dangqian = -1 * A(i,k);
end
%如果|A(i,k)|>max,将max改成该值，并记录该行
if dangqian > max
max = dangqian;
maxrow = i;
end
end
%从A(k,k)到A(n,k)元素都为0，A为奇异阵，退出
if max == 0
break;
end
%如果A(k,k)不是从A(k,k)到A(n,k)中绝对值最大的元素,将矩阵C的maxrow与k行互换
if maxrow ~= k
d = C(k,:);
C(k,:) = C(maxrow,:);
C(maxrow,:) = d;
end
%将C的C(k,k)应用行变换变为1
ckk = C(k,k);
for j = 1:2*n
C(k,j) = C(k,j)/ckk;
end
%将C第k列除C(k,k)应用行变换全变为0
for i = 1:n
if i == k
continue;
end
mik = C(i,k)/C(k,k);
%将第i行减去mik倍的第k行
for j = 1:2*n
C(i,j) = C(i,j) - mik * C(k,j);
end
end
end
%invA为A的逆
invA = C(1:n,5:8)

////////////////////////
以下内容为回复内容 /////
////////////////////////



friend726

斗胆在这说一下自己的想法吧。
首先，王静写的matlab程序，非常棒。巧妙的标记变量设置，程序层次清楚，注释恰当。非常欣赏。
二。提一下建议。
（1）最好将此程序设为函数，可使此程序的重用性提高。
例如：function [inv] inverse(input);
inv = zeros(size(input));
inv = C(1:n,n+1:2*n);
（2）下面语句中：我建议将“1”改为“k”。我认为会减少一些运算量。
%将C的C(k,k)应用行变换变为1
ckk = C(k,k);
for j = 1:2*n
C(k,j) = C(k,j)/ckk;
end
%将C第k列除C(k,k)应用行变换全变为0
（3）程序中最后一句，我认为改为 C(1:n,n+1:2*n)比较好。
其实我的方法，就是将原程序，转为函数，提高重用性。
三。至于xuqiang老师的问题。我认为举个例子最好说明问题。
例如：
（0.000000000000000000333333333333333，2；10，4）
10除于很小的数，误差是很大的。（或者很小的数做分母误差就会很大了。）
优点：选主元方法，就会很好的避免A（k，k）很小的情况，使数值精确和稳定性更好。
至于缺点，就是运算量更多一些。（多了n（n-1）/2步比较和至多n-1次行交换。我算的是这个。）

其实我在此说的都是多余的，xuqiang老师的问题，数值老师有讲。而我的对于程序的改进，只要作者想变为函数，自然就会解决。
所以，My name is nobody.



Nobody

补充一下，之所以高等代数没有选主元，而数值有选，主要是高代是理想的处理，而计算机确无法存储趋于无穷的小数或者无穷大的数。


XuQiang
选主元和变程序为函数的事情nobody说的很清楚了。

关于（2）下面语句中：我建议将“1”改为“k”。我认为会减少一些运算量。
%将C的C(k,k)应用行变换变为1
ckk = C(k,k);
for j = 1:2*n
C(k,j) = C(k,j)/ckk;
end
%将C第k列除C(k,k)应用行变换全变为0

我没有看明白是什么意思，这个消元求逆矩阵的过程是需要行变换为单位阵的，不知Nobody的理解是否有失偏颇。

最后，U r not nobody, u will be somebody someday.


friend726


至于老师提出的问题，我想举个例子吧。
例如：
A = [1, 2, 3; 1, 4, 2 ;1, 3, 5;] （只是举个例子，是不是非奇异矩阵暂且不管吧，只是说一下过程）
约化第一列会将A约化为 A = [1, 2, 3; 0, 2, -1; 0, 1, 3;]
约化第二列时，要将第二行单位化，但是A（2，1）是 0 呀！，根本不需要单位化的。而程序中j是从1开始的。如果改为k，就会避免0的单位化，减少一些运算。我想虽然是0，但是Computer还是要运算一次的吧。
谢谢老师的回复。有不当之处还望不吝指教。

XuQiang

为了节省计算量，是这样，没错。