Reduced row echelon form technique

25 次查看(过去 30 天)
I am trying to use a code to calculate the reduced row echelon form of a matrix without the function rref. I make a random matrix A and and then make a matrix new_A = (A-lambda*I). The intent is to eventually find the nullspace of new_A without the null function. When comparing the code to the rref command the results match a majority of the time, however on certain occasions there is a discrepancy in the final column. My understanding is that rref must be unique for a given matrix. Any ideas on what could be the source?
clear all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAKES RANDOM MATRIX A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m = 3;
casenum = randi(3,1)
if casenum == 1
A = randn(m);
end
if casenum == 2
A = i*randn(m);
end
if casenum == 3
pownum = randi(2,m);
A = randn(m) + randn(m).*(i.^pownum);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAKES MATRIX new_A FROM (A-lambda*I) & CALCULATES MATLAB NULL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I = eye(m,m);
mat_eigs = eig(A);
new_A = A - mat_eigs(1)*I % currently only uses 1st eigenvalue to test
mat_null = null(new_A);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FINDS RREF OF new_A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mat_rref = rref(new_A)
tol = 1e-12;
i = 1;
j = 1;
while (i <= m) && (j <= m)
% Find value and index of largest element in the remainder of column j
[p,k] = max(abs(new_A(i:m,j)));
k = k+i-1;
if (p <= tol)
% The column is negligible, zero it out
new_A(i:m,j) = 0;
j = j + 1;
else
% Swap i-th and k-th rows
new_A([i k],j:m) = new_A([k i],j:m);
% Divide the pivot row by the pivot element
Ai = new_A(i,j:m) / new_A(i,j);
% Subtract multiples of the pivot row from all the other rows
new_A(:,j:m) = new_A(:,j:m) - new_A(:,j)*Ai;
new_A(i,j:m) = Ai;
i = i + 1;
j = j + 1;
end
end
code_rref = new_A

采纳的回答

Richard Brown
Richard Brown 2013-8-7
It'll be a difference between yours and MATLAB's choice of tolerance. Your matrix should be rank 2, however I notice that MATLAB is often computing the RREF to be the identity.
If you read the docs for rref they specify what they use as tolerance:
max(size(A))*eps *norm(A,inf))
  1 个评论
Lingling Fan
Lingling Fan 2019-6-13
Hi,
Thanks for your answer!
How could we change the tolerance of rref, as I did not find in docs.
Is there something like "options" like fmincon setup for tolerance?
Thanks in advance!

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by