Matlab函数rref()不能求解线性方程组的精确解

最新推荐文章于 2022-12-29 17:10:56 发布

underratedtang

最新推荐文章于 2022-12-29 17:10:56 发布

阅读量1.3k

点赞数 3

文章标签： matlab 线性代数开发语言

本文链接：https://blog.csdn.net/xxxtlxxx/article/details/122056109

版权

一个同学分享给我一个有趣的问题，我们先来看他的代码：

%Codes for Matlab r2020b
clc;
clear;
A=[6 5 0 0 0 0 0 0 0 0 1 2
   1 6 5 0 0 0 0 0 0 0 1 0
   0 1 6 5 0 0 0 0 0 0 1 0
   0 0 1 6 5 0 0 0 0 0 1 1
   0 0 0 1 6 5 0 0 0 0 1 5
   0 0 0 0 1 6 5 0 0 0 1 0
   0 0 0 0 0 1 6 5 0 0 1 9
   0 0 0 0 0 0 1 6 5 0 1 0
   0 0 0 0 0 0 0 1 6 5 1 8
   0 0 0 0 0 0 0 0 1 6 1 0];
b=[1 2 3 4 5 6 7 8 9 10]';
B=[A,b];
rref(B)
b1=[0;0;0;0;0;0;0;0;0;0];
B1=[A b1];
rref(B1)

很明显，他想用化简增广矩阵至行最简型矩阵的方法得出一个齐次和一个非齐次线性方程组的解，这是我们在学习高等代数时求解类似问题的常用方法。可是在Matlab计算时出现了一个小问题。

他的代码跑出的结果：

他发现最终算出来的最简行矩阵中B1和B的第十一列第五行不一样，一个是0，一个是分数（小数）。按照正常情况，两个增广矩阵除了最后一列以外，其他位置应该完全一样。因为我们只会对矩阵做行变换而不会做列变换。

问题可能出在Matlab自带函数rref()的源代码上。在命令行窗口输入open rref即可打开rref函数源代码的.m文件。代码如下：

function [A, jb] = rref(A, tol)
%RREF   Reduced row echelon form.
%   R = RREF(A) produces the reduced row echelon form of A.
%
%   [R,jb] = RREF(A) also returns a vector, jb, so that:
%       r = length(jb) is this algorithm's idea of the rank of A,
%       x(jb) are the bound variables in a linear system, Ax = b,
%       A(:,jb) is a basis for the range of A,
%       R(1:r,jb) is the r-by-r identity matrix.
%
%   [R,jb] = RREF(A,TOL) uses the given tolerance in the rank tests.
%   [R,jb] = RREF(A,TOL) 在等级测试中使用给定的容差。
%
%   Roundoff errors may cause this algorithm to compute a different
%   value for the rank than RANK, ORTH and NULL.
%   舍入错误可能会导致此算法为排名计算与 RANK、ORTH 和 NULL 不同的值。
%
%   Class support for input A:
%      float: double, single
%
%   See also RANK, ORTH, NULL, QR, SVD.

%   Copyright 1984-2017 The MathWorks, Inc. 

[m,n] = size(A);

% Does it appear that elements of A are ratios of small integers?
[num, den] = rat(A);
rats = isequal(A, num./den);

% Compute the default tolerance if none was provided.如果没有提供，计算默认容差。
if (nargin < 2)
    tol = max(m,n)*eps(class(A))*norm(A,inf);
end

% Loop over the entire matrix.
i = 1;
j = 1;
jb = zeros(1,0);
while i <= m && j <= n
   % Find value and index of largest element in the remainder of column j.
   [p, k] = max(abs(A(i:m,j)));
   k = k+i-1;
   if p <= tol
      % The column is negligible, zero it out.该列可以忽略不计，将其归零。             
      A(i:m,j) = 0;
      j = j + 1;
   else
      % Remember column index
      jb = [jb j]; %#ok<AGROW>
      % Swap i-th and k-th rows.
      A([i k],j:n) = A([k i],j:n);
      % Divide the pivot row by the pivot element.
      A(i,j:n) = A(i,j:n)./A(i,j);
      % Subtract multiples of the pivot row from all the other rows.
      for k = [1:i-1 i+1:m]
         A(k,j:n) = A(k,j:n) - A(k,j).*A(i,j:n);
      end
      i = i + 1;
      j = j + 1;
   end
end

% Return "rational" numbers if appropriate.
if rats
    [num, den] = rat(A);
    A = num./den;
end

查看备注我们可以发现，rref函数在计算时进行了多次减少计算量的运算，这难免产生了误差。比如使用rat函数（此 MATLAB 函数返回 X 在默认容差 1e-6*norm(X(:),1) 内的有理分式近似值。近似值是一个包含截断的连分式展开项的字符向量。）。rref函数的help文档也提出了算法对于很小误差的忽略：

rref 通过部分主元消元法实现 Gauss-Jordan 消元法。(max(size(A))*eps*norm(A,inf)) 的默认容差检验可忽略不计的列元素，这些列元素会归零以减少舍入误差。

但不是所有结果都有误差，所以rref的结果作为参考是一个很好的选择。