五猴分桃-诺贝尓奖获得者出过的题目用北太天元也能轻松求解

把矩阵化成简化行阶梯型矩阵的函数rref 在北太天元中已经实现为内置函数了。下面是一个脚本实现的版本，可以供写脚本的同志们参考。

为了和内置的rref区分，这里的函数名改成了 rrefEx, 这个脚本参考了octave的rref.m的实现。

北太天元作简化行阶梯型矩阵的rrefEx.m 的内容如下

% 把矩阵 A 化成 简化行阶梯形矩阵

% [A, k] = rrefEx (A, tol)

% 右边的A 是输入的矩阵，可以是single或者double，但是必须是二维矩阵

% tol 是忍量，矩阵值比tol小的数会被置为0

% 左边的A 是输出的矩阵，是简化行阶梯形矩阵

% k 是主元所在的列号

% 例:

% a = [1];

% [r k] = rrefEx (a);

% a = [1 3; 4 5];

% [r k] = rrefEx (a);

% a = [1 3; 4 5; 7 9];

% [r k] = rrefEx (a);

% a = [1 2 3; 2 4 6; 7 2 0];

% [r k] = rrefEx (a);

% a = [1 2 1; 2 4 2.01; 2 4 2.1];

% tol = 0.02;

% [r k] = rrefEx (a, tol);

% tol = 0.2;

% [r k] = rrefEx (a, tol);

function [A, k] = rrefEx (A, tol)

   %如果输入的参数个数小于1，那么输出帮助信息

 if (nargin < 1)

      help rrefEx;

 end

 if (ndims (A) > 2)

   error ("rrefEx: A 必须是一个二维矩阵");

 end

 [rows, cols] = size (A);

 if (nargin < 2)

   if (isa (A, 'single'))

     tol = eps ('single') * max (rows, cols) * norm (A, inf ('single'));

   else

     tol = eps * max (rows, cols) * norm (A, inf);

   end

 end

 used = zeros (1, cols); %用来记录用过的列

 r = 1; 

 for c = 1:cols

   % 第c列的 从第r行到最后一行 的主元 （绝对值最大的元素)

   [m, pivot] = max (abs (A(r:rows,c))); % m 是主元的大小

   pivot = r + pivot - 1; % pivot 是主元所在的行号

   if (m <= tol)

     % 如果上面得到的主元m很小，那么跳过第c列，

       % 同时把这些小元素都强制修改成真正的0

     A(r:rows, c) = zeros (rows-r+1, 1);

   else

     % 记录主元所在的列号

     used(1, c) = 1;

     % 交换当前行和主元所在的行

     A([pivot, r], c:cols) = A([r, pivot], c:cols);

     % 现在主元所在的行已经放在了当前行r了，再用一个常数1/A(r,c)乘以这一行

         % 从而把主元 标准化为 1.

     A(r, c:cols) = A(r, c:cols) / A(r, c);

     % 把当前列的主元上下方的元素都消成0

     ridx = [1:r-1, r+1:rows];% 1:r-1 是主元上方的行号，r+1:rows是主元下方的行号

     A(ridx, c:cols) = A(ridx, c:cols) - A(ridx, c) * A(r, c:cols);

     %检查是否完成了所有行

     if (r++ == rows)

       break;

     end

   end

 end

 k = find (used);

end