Matlab erreur de précision de la fonction Rref() après le 12 colonne des matrices hilbert

Question

Ma question peut être un simple mais je ne pouvais pas penser à une explication logique à ma question:

Quand j'utilise

rref(hilb(8)), rref(hilb(9)), rref(hilb(10)), rref(hilb(11)) 

ça me donne le résultat auquel je m'attendais, une matrice unitaire.

Cependant, quand il vient à la

rref(hilb(12)) 

il ne donne pas une inversible matrice comme prévu. J'ai utilisé Wolfram et il donne la matrice unité pour le même cas donc je suis sûr qu'il aurait dû donner une matrice unitaire. Il peut y avoir une erreur d'arrondi ou quelque chose comme ça, mais 1/11 ou 1/7 ont aussi quelques décimales gênantes

alors pourquoi Matlab se comporte-t-il comme ça quand il s'agit de 12?

Answer 1

Cela ressemble en effet à une erreur de précision. Ceci est logique car le déterminant de la matrice de Hilbert d'ordre n tend vers 0 car n tend vers l'infini (see here). Cependant, vous pouvez utiliser rref with tol parameter:

[R,jb] = rref(A,tol) 

et prendre tol être très petit pour obtenir des résultats plus précis. Par exemple, rref(hilb(12),1e-20) vous donnera une matrice d'identité.

EDIT- plus de détails concernant le rôle du paramètre tol. Le code source rref est fourni au bas de la réponse. Le tol est utilisé lorsque nous recherchons un élément maximal en valeur absolue dans une certaine partie d'une colonne, pour trouver la ligne de pivot.

% 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) = zeros(m-i+1,1); 
     j = j + 1; 

Si tous les éléments sont plus petits que tol en valeur absolue, la partie correspondante de la colonne est remplie par des zéros. Cela semble être l'endroit où l'erreur de précision pour rref(hilb(12)) se produit. En réduisant le tol, nous évitons ce problème en rref(hilb(12),1e-20).

code source:

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. 
% 
% Roundoff errors may cause this algorithm to compute a different 
% value for the rank than RANK, ORTH and NULL. 
% 
% Class support for input A: 
%  float: double, single 
% 
% See also RANK, ORTH, NULL, QR, SVD. 

% Copyright 1984-2005 The MathWorks, Inc. 
% $Revision: 5.9.4.3 $ $Date: 2006/01/18 21:58:54 $ 

[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 = []; 
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) = zeros(m-i+1,1); 
     j = j + 1; 
    else 
     % Remember column index 
     jb = [jb j]; 
     % 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