Skip to content

Commit 5285a3d

Browse files
authored
Add a linear system solver (#6196)
1 parent 769e497 commit 5285a3d

File tree

2 files changed

+92
-0
lines changed

2 files changed

+92
-0
lines changed
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
package com.thealgorithms.matrix;
2+
3+
/**
4+
* This class implements an algorithm for solving a system of equations of the form Ax=b using gaussian elimination and back substitution.
5+
*
6+
* @link <a href="https://en.wikipedia.org/wiki/Gaussian_elimination">Gaussian Elimination Wiki</a>
7+
* @see InverseOfMatrix finds the full of inverse of a matrice, but is not required to solve a system.
8+
*/
9+
public final class SolveSystem {
10+
private SolveSystem() {
11+
}
12+
13+
/**
14+
* Problem: Given a matrix A and vector b, solve the linear system Ax = b for the vector x.\
15+
* <p>
16+
* <b>This OVERWRITES the input matrix to save on memory</b>
17+
*
18+
* @param matrix - a square matrix of doubles
19+
* @param constants - an array of constant
20+
* @return solutions
21+
*/
22+
public static double[] solveSystem(double[][] matrix, double[] constants) {
23+
final double tol = 0.00000001; // tolerance for round off
24+
for (int k = 0; k < matrix.length - 1; k++) {
25+
// find the largest value in column (to avoid zero pivots)
26+
double maxVal = Math.abs(matrix[k][k]);
27+
int maxIdx = k;
28+
for (int j = k + 1; j < matrix.length; j++) {
29+
if (Math.abs(matrix[j][k]) > maxVal) {
30+
maxVal = matrix[j][k];
31+
maxIdx = j;
32+
}
33+
}
34+
if (Math.abs(maxVal) < tol) {
35+
// hope the matrix works out
36+
continue;
37+
}
38+
// swap rows
39+
double[] temp = matrix[k];
40+
matrix[k] = matrix[maxIdx];
41+
matrix[maxIdx] = temp;
42+
double tempConst = constants[k];
43+
constants[k] = constants[maxIdx];
44+
constants[maxIdx] = tempConst;
45+
for (int i = k + 1; i < matrix.length; i++) {
46+
// compute multipliers and save them in the column
47+
matrix[i][k] /= matrix[k][k];
48+
for (int j = k + 1; j < matrix.length; j++) {
49+
matrix[i][j] -= matrix[i][k] * matrix[k][j];
50+
}
51+
constants[i] -= matrix[i][k] * constants[k];
52+
}
53+
}
54+
// back substitution
55+
double[] x = new double[constants.length];
56+
System.arraycopy(constants, 0, x, 0, constants.length);
57+
for (int i = matrix.length - 1; i >= 0; i--) {
58+
double sum = 0;
59+
for (int j = i + 1; j < matrix.length; j++) {
60+
sum += matrix[i][j] * x[j];
61+
}
62+
x[i] = constants[i] - sum;
63+
if (Math.abs(matrix[i][i]) > tol) {
64+
x[i] /= matrix[i][i];
65+
} else {
66+
throw new IllegalArgumentException("Matrix was found to be singular");
67+
}
68+
}
69+
return x;
70+
}
71+
}
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
package com.thealgorithms.matrix;
2+
3+
import static org.junit.jupiter.api.Assertions.assertArrayEquals;
4+
5+
import java.util.stream.Stream;
6+
import org.junit.jupiter.params.ParameterizedTest;
7+
import org.junit.jupiter.params.provider.Arguments;
8+
import org.junit.jupiter.params.provider.MethodSource;
9+
10+
class SolveSystemTest {
11+
12+
@ParameterizedTest
13+
@MethodSource({"matrixGenerator"})
14+
void solveSystem(double[][] matrix, double[] constants, double[] solution) {
15+
double[] expected = SolveSystem.solveSystem(matrix, constants);
16+
assertArrayEquals(expected, solution, 1.0E-10, "Solution does not match expected");
17+
}
18+
private static Stream<Arguments> matrixGenerator() {
19+
return Stream.of(Arguments.of(new double[][] {{-5, 8, -4}, {0, 6, 3}, {0, 0, -4}}, new double[] {38, -9, 20}, new double[] {-2, 1, -5}), Arguments.of(new double[][] {{-2, -1, -1}, {3, 4, 1}, {3, 6, 5}}, new double[] {-11, 19, 43}, new double[] {2, 2, 5}));
20+
}
21+
}

0 commit comments

Comments
 (0)