Skip to content

Commit 620ec4e

Browse files
authored
feat(solvers): implement Gaussian Elimination with Pivoting
- Implements forward elimination with partial pivoting for numerical stability. - Includes checks for singular or inconsistent systems. - Replaces SymPy back-substitution with a pure MATLAB vectorized loop. - Returns the solution vector x.
1 parent b42ad5e commit 620ec4e

File tree

1 file changed

+47
-0
lines changed

1 file changed

+47
-0
lines changed
Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,47 @@
1+
% +linearSolvers/gaussianElimination.m
2+
function x = gaussianElimination(A, b)
3+
% Solves Ax = b using Gaussian Elimination with partial pivoting
4+
% followed by back-substitution.
5+
6+
[n, ~] = size(A);
7+
Ab = [A, b(:)]; % Augmented matrix
8+
9+
% Forward Elimination with Partial Pivoting
10+
for i = 1:n-1
11+
% Find pivot row
12+
[~, maxRowIdx] = max(abs(Ab(i:n, i)));
13+
maxRowIdx = maxRowIdx + i - 1; % Adjust index relative to full matrix
14+
15+
% Swap rows if necessary
16+
if maxRowIdx ~= i
17+
Ab([i, maxRowIdx], :) = Ab([maxRowIdx, i], :);
18+
end
19+
20+
% Check for singular matrix (zero pivot after pivoting)
21+
if abs(Ab(i, i)) < eps % Use epsilon for floating point comparison
22+
error('Matrix is singular or nearly singular. Cannot solve using Gaussian Elimination.');
23+
end
24+
25+
% Eliminate column i in rows below i
26+
for j = i+1:n
27+
factor = Ab(j, i) / Ab(i, i);
28+
Ab(j, i:n+1) = Ab(j, i:n+1) - factor * Ab(i, i:n+1);
29+
end
30+
end
31+
32+
% Check for inconsistency in the last row
33+
if abs(Ab(n, n)) < eps && abs(Ab(n, n+1)) > eps
34+
error('System has no solution (inconsistent).');
35+
elseif abs(Ab(n, n)) < eps && abs(Ab(n, n+1)) < eps
36+
error('System has infinitely many solutions.'); % Or handle appropriately
37+
end
38+
39+
% Back Substitution (Pure MATLAB, no SymPy)
40+
x = zeros(n, 1);
41+
x(n) = Ab(n, n+1) / Ab(n, n);
42+
for i = n-1:-1:1
43+
sum_ax = Ab(i, i+1:n) * x(i+1:n);
44+
x(i) = (Ab(i, n+1) - sum_ax) / Ab(i, i);
45+
end
46+
47+
end

0 commit comments

Comments
 (0)