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
0 个评论
采纳的回答
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
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 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!