Skip to content

Commit c90a5ae

Browse files
committed
Added a linear system solver
1 parent e6073f8 commit c90a5ae

File tree

2 files changed

+95
-0
lines changed

2 files changed

+95
-0
lines changed
Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
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 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+
* This variation uses @link <a href="https://en.wikipedia.org/wiki/Crout_matrix_decomposition">Crout Reduction</a>
17+
* and partial pivoting to decompose the matrix.
18+
* <b>This OVERWRITES the input matrix to save on memory</b>
19+
*
20+
* @param matrix - a square matrix of doubles
21+
* @param constants - an array of constant
22+
* @return solutions
23+
*/
24+
public static double[] solveSystem(double[][] matrix, double[] constants) {
25+
final double TOL = 0.00000001; // tolerance for round off
26+
for (int k = 0; k < matrix.length - 1; k++) {
27+
// find the largest value in column (to avoid zero pivots)
28+
double maxVal = Math.abs(matrix[k][k]);
29+
int maxIdx = k;
30+
for (int j = k + 1; j < matrix.length; j++) {
31+
if (Math.abs(matrix[j][k]) > maxVal) {
32+
maxVal = matrix[j][k];
33+
maxIdx = j;
34+
}
35+
}
36+
if (Math.abs(maxVal) < TOL) // hope the matrix works out
37+
continue;
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+
46+
for (int i = k + 1; i < matrix.length; i++) {
47+
// compute multipliers and save them in the column
48+
49+
matrix[i][k] /= matrix[k][k];
50+
for (int j = k + 1; j < matrix.length; j++) {
51+
matrix[i][j] -= matrix[i][k] * matrix[k][j];
52+
}
53+
constants[i] -= matrix[i][k] * constants[k];
54+
}
55+
}
56+
// back substitution
57+
double[] x = new double[constants.length];
58+
System.arraycopy(constants, 0, x, 0, constants.length);
59+
for (int i = matrix.length - 1; i >= 0; i--) {
60+
double sum = 0;
61+
for (int j = i + 1; j < matrix.length; j++) sum += matrix[i][j] * x[j];
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+
return x;
69+
}
70+
}
Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
package com.thealgorithms.matrix;
2+
3+
import org.junit.jupiter.params.ParameterizedTest;
4+
import org.junit.jupiter.params.provider.Arguments;
5+
import org.junit.jupiter.params.provider.MethodSource;
6+
7+
import java.util.stream.Stream;
8+
9+
import static org.junit.jupiter.api.Assertions.assertArrayEquals;
10+
11+
class SolveSystemTest {
12+
13+
@ParameterizedTest
14+
@MethodSource({"matrixGenerator"})
15+
void solveSystem(double[][] matrix, double[] constants, double[] solution) {
16+
17+
double[] expected = SolveSystem.solveSystem(matrix, constants);
18+
assertArrayEquals(expected, solution, 1.0E-10, "Solution does not match expected");
19+
}
20+
21+
private static Stream<Arguments> matrixGenerator() {
22+
23+
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}));
24+
}
25+
}

0 commit comments

Comments
 (0)