This commit is contained in:
Sergiusz Warga 2021-03-06 20:58:15 +01:00
parent 87f957d903
commit 7809b5d87b
3 changed files with 33 additions and 16 deletions

View File

@ -1,6 +1,4 @@
function A = gaussian_elimination_with_complete_pivoting(A)
% det(A(mi, lm))
function [U, L] = gaussian_elimination_with_complete_pivoting(A)
[n, m] = size(A);
if n ~= m
@ -11,21 +9,35 @@ if det(A) == 0
error('Matrix is not nonsingular!')
end
A
% for k = 1 : n-1
for k = 1 : n-1
i = k:n;
j = k:n;
A(i, j);
maximum = max(abs(A(i, j)), [], 'all');
max_idx = find(abs(A==maximum));
[mi, lm] = ind2sub(size(A), max_idx(1));
A([k mi], 1:n) = deal(A([mi k], 1:n));
A(1:n, [k lm]) = deal(A(1:n, [lm k]));
p(k) = mi;
q(k) = lm;
A(i, j)
[max_val, rows_of_max_in_col] = max(abs(A(i, j)));
[max_val, max_col] = max(max_val);
max_row = rows_of_max_in_col(max_col);
% Assing value of mi and lambda in respect to the main A matrix
[mi, lm] = deal(max_row+k-1, max_col+k-1)
A([k mi], 1:n) = deal(A([mi k], 1:n))
A(1:n, [k lm]) = deal(A(1:n, [lm k]))
p(k) = mi
q(k) = lm
% Perform Gaussian elimination with the greatest pivot
if A(k, k) ~= 0
rows = k+1 : n;
A(rows, k) = A(rows, k)/A(k, k);
A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows);
end
end
U = triu(A);
L = tril(A, -1) + eye(n);
p
I = eye(n);
P = I(p, :)
q
Q = I(:, q)

View File

@ -1,10 +1,15 @@
clear all;
B = [2, -1, 0, 0; -1, 2, -1, 0;0, -1, 2, -1; 0, 0, -1, 2];
B = [2, -1, 0, 0;
-1, 2, -1, 0;
3, -1, 2, -1;
0, 4, -1, 2];
b = [0;0;0;5];
[U, Mk] = outer_product_gaussian_elimination(B);
[U, L] = outer_product_gaussian_elimination(B);
back_substitution(U, b);
A = gaussian_elimination_with_complete_pivoting(B);
[U, L] = gaussian_elimination_with_complete_pivoting(B)
L*U
% A = gauss_jordan_elimination(B)

View File

@ -1,5 +1,5 @@
% Algorithm 1: Outer Product Gaussian Elimination (Alg. 3.2.1)
function [U, Mk] = outer_product_gaussian_elimination(A)
function [U, L] = outer_product_gaussian_elimination(A)
[n, m] = size(A);
if n ~= m
@ -17,4 +17,4 @@ end
end
U = triu(A);
Mk = diag(A,-1); % Gauss vector
L = tril(A, -1) + eye(n);