# Treating ties in a Gaussian Elimination process

7 views (last 30 days)
Conor Pierce on 16 Apr 2024
Commented: John D'Errico on 16 Apr 2024
Hey there, I have an M-File written for performing Gaussian elimination (shown below) however I can't figure out how to deal with ties in the matrix for example if there are two elements that are both the max value being called on how do I get the function to just select the first one it comes across when finding the max. Any help is appreciated.
function [Ang_update] = Gauss_elim(f, Df, k)
% [Ang_update] = Gauss_elim(f, Df, k)
% Using Gaussian Elimination to calculate the next approximation of the
% Newton Raphson iteration
% Input f is the column vector of functions containing the unknown variables
% Input Df is the Jacobian matrix of the functions in question
% Output Ang_update is the list of angles at the most recent iteration
% Version 1, 08/04/2024 by Conor Pierce, Student Number 21302251
if ~isvector(f)
error('Input f is a column vector of the system of equation to be solved.')
end
if ~isvector(k)
error('k is a column vector (8 x 1) containing the users initial estimates as entered in the Newton Raphson function.')
end
k0 = k; % Vector of initial estimates
F = f; % Initialise f(kn)
% Finding the size of the matrix
s = size(F);
n = s(1);
% Set Df as the Jacobian and formulate succesive steps
J = Df;
A = J;
B = Df*k0 - F;
% Begin set up for Gaussian elimination
A1 = A;
B1 = B;
% Set up identity matrix
Identity = eye(n);
Q_total = Identity;
% Total Pivoting
for m = 1 : n-1
active = m:n; % Identify active rows and columns
[Y, I] = max(abs(A1(active, active))); % Y is max value in column, I is position of that value
[~, I1] = max(Y); % Column to pivot is I1
row_num = I(I1);
col_num = I1;
p = 1:n;
p([m, row_num + (m-1)]) = p([row_num + (m-1), m]); % Row swap
q = 1:n;
q([m, col_num + (m-1)]) = q([col_num + (m-1), m]); % Column swap
Q_total = Q_total(:, q);
A1 = A1(p, q);
B1 = B1(p, :);
% Subtraction
I_vec = zeros(1, n);
I_vec(1, m) = 1;
L1 = Identity - ([zeros(m, 1); A1(m+1: n,m)]*I_vec)/A1(m, m);
A1 = L1*A1;
B1 = L1*B1;
end
% Set up column vector for outputs
Ang_update = zeros(n, 1);
for m = 0:n-1
% Back substitution
Ang_update(n-m, 1) = (B1(n-m, 1) - (A1(n-m, :)*Ang_update))/A1(n-m, n-m);
end
Ang_update = Q_total*Ang_update; % Reorder
end

John D'Errico on 16 Apr 2024
Um, max does exactly that. Where is there a problem?
x = [1 2 3 4 3 4 1];
[xmax,ind] = max(x)
xmax = 4
ind = 4
Conor Pierce on 16 Apr 2024
Yes, but does this hold for an 8 x 8 matrix?
John D'Errico on 16 Apr 2024
Ok, you are using complete pivoting, not just row pivoting as is common. Max, by default, finds the maximum in each column, and you don't want that.
A = randi(10,5,5)
A = 5x5
1 4 3 3 6 2 8 6 2 6 7 8 8 3 9 3 5 2 8 4 3 8 6 9 1
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
[Amax,ind] = max(A)
Amax = 1x5
7 8 8 9 9
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
ind = 1x5
3 2 3 5 3
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
But if you do this:
[Amax,ind] = max(A,[],'all')
Amax = 9
ind = 20
It finds the largest overall element in the array. The index returned is a linear index, so the 9th element when the array is unrolled into a vector. You can recover the row and column indices of that by
[imax,jmax] = ind2sub(size(A),ind)
imax = 5
jmax = 4
Note that there were 2 instances of a 9. The first one was seen as A(5,4).