Problem with reordering using Gaussian Elimination.

Hi, I am writing code to solve a system of equations using Gaussian Elimination, and partial pivoting. My code is below and is returning me with the correct values for each co-efficient, however the reordering code I have tried to use is not working. Thanks!

Ra = size(A); n = Ra(1); Id = eye(n);

P_tot = Id; A1 = A; b1 = b;

% Checking size of inputs is correst

if ~isequal(size(A),[n n])

error("This function can only be used with square (n x n) matrices of co-efficients.")

end

if ~isequal(size(b),[n 1])

error("This function can only be used with n x 1 vectors of constants.")

end

for m = 1:n-1

% Selecting row containing pivotal co-efficeint.

row_index = m:n;

[~,Ir] = max(abs(A1(row_index, row_index)));

row_num = Ir(1);

% Partial Pivoting

p = [1:n]; p(m) = m-1+row_num; p(m-1+row_num) = m;

P = Id(p,:);

A1 = P*A1; b1 = P*b1; P_tot = P*P_tot;

% eliminate

I_vec = zeros(1,n); I_vec(1,m) = 1;

L1 = Id - ([zeros(m,1); A1(m+1:n,m)]*I_vec)/A1(m,m);

A1 = L1*A1; b1 = L1*b1;

end

% Back Substitiution

v = zeros(n,1);

for m = 0:n-1

v(n-m,1) = (b1(n-m,1) - (A1(n-m,:)*v))/A1(n-m,n-m);

end

v = P_tot*v; % reorder

end

[Alan Stevens]

What do you mean by "not working"? What do you expect the result to be (i.e. what should the new order be?)?

Answer 1

[Rohit Kulkarni]

Hi Sean,

In my understanding the issue in your MATLAB code for Gaussian elimination with partial pivoting lies in the calculation and application of row swaps during partial pivoting. Specifically, the way you're calculating 'row_num' and updating the permutation vector 'p' doesn't correctly reflect the intended row swaps in the context of the entire matrix. Here's a suggested correction:

Original Issue:

- Incorrect calculation of 'row_num' leading to improper row swaps.

Correction:

- Adjust 'row_num' to reflect its absolute position within the entire matrix.

- Correctly swap rows in the permutation vector 'p'.

Corrected Code Snippet:

% Correctly identifying the pivot row within the entire matrix
row_index = m:n;
[~,Ir] = max(abs(A1(row_index, m)));
row_num = row_index(Ir(1)); % Absolute row number for the pivot

% Partial Pivoting
p = 1:n;
p([m, row_num]) = p([row_num, m]); % Swap the intended rows
P = Id(p,:);
A1 = P*A1; b1 = P*b1; P_tot = P*P_tot;

Hope it helps!

Thanks