diff --git a/src/main/java/com/thealgorithms/maths/LUDecomposition.java b/src/main/java/com/thealgorithms/maths/LUDecomposition.java new file mode 100644 index 000000000000..3f34350720af --- /dev/null +++ b/src/main/java/com/thealgorithms/maths/LUDecomposition.java @@ -0,0 +1,129 @@ +package com.thealgorithms.maths; + +/** + * @brief Implementation of LU Decomposition using the Doolittle algorithm + * @details Decomposes a square matrix A into a lower triangular matrix L and + * an upper triangular matrix U such that A = L * U. The diagonal of L contains + * all ones (Doolittle convention). This decomposition is useful for solving + * systems of linear equations, computing determinants, and finding inverses. + * @see LU Decomposition + */ +public final class LUDecomposition { + + private LUDecomposition() { + } + + /** + * @brief Performs LU decomposition on a square matrix using the Doolittle algorithm + * @param matrix a square matrix + * @return a 2D array where the lower triangle (excluding diagonal) contains L + * elements (with implicit 1s on the diagonal) and the upper triangle + * (including diagonal) contains U elements + * @throws IllegalArgumentException if the matrix is not square + * @throws ArithmeticException if a zero pivot is encountered + */ + public static double[][] decompose(double[][] matrix) { + int n = matrix.length; + for (double[] row : matrix) { + if (row.length != n) { + throw new IllegalArgumentException("Matrix must be square."); + } + } + + double[][] lu = new double[n][n]; + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + lu[i][j] = matrix[i][j]; + } + } + + for (int k = 0; k < n; k++) { + for (int j = k; j < n; j++) { + double sum = 0; + for (int s = 0; s < k; s++) { + sum += lu[k][s] * lu[s][j]; + } + lu[k][j] -= sum; + } + + if (lu[k][k] == 0) { + throw new ArithmeticException("Zero pivot encountered. Matrix may be singular."); + } + + for (int i = k + 1; i < n; i++) { + double sum = 0; + for (int s = 0; s < k; s++) { + sum += lu[i][s] * lu[s][k]; + } + lu[i][k] = (lu[i][k] - sum) / lu[k][k]; + } + } + + return lu; + } + + /** + * @brief Extracts the lower triangular matrix L from the combined LU matrix + * @param lu the combined LU matrix from decompose() + * @return the lower triangular matrix L with 1s on the diagonal + */ + public static double[][] getLowerMatrix(double[][] lu) { + int n = lu.length; + double[][] lower = new double[n][n]; + for (int i = 0; i < n; i++) { + lower[i][i] = 1.0; + for (int j = 0; j < i; j++) { + lower[i][j] = lu[i][j]; + } + } + return lower; + } + + /** + * @brief Extracts the upper triangular matrix U from the combined LU matrix + * @param lu the combined LU matrix from decompose() + * @return the upper triangular matrix U + */ + public static double[][] getUpperMatrix(double[][] lu) { + int n = lu.length; + double[][] upper = new double[n][n]; + for (int i = 0; i < n; i++) { + for (int j = i; j < n; j++) { + upper[i][j] = lu[i][j]; + } + } + return upper; + } + + /** + * @brief Solves a system of linear equations Ax = b using LU decomposition + * @param lu the combined LU matrix from decompose() + * @param b the right-hand side vector + * @return the solution vector x + */ + public static double[] solve(double[][] lu, double[] b) { + int n = lu.length; + double[] y = new double[n]; + double[] x = new double[n]; + + // Forward substitution: solve Ly = b + for (int i = 0; i < n; i++) { + double sum = 0; + for (int j = 0; j < i; j++) { + sum += lu[i][j] * y[j]; + } + y[i] = b[i] - sum; + } + + // Back substitution: solve Ux = y + for (int i = n - 1; i >= 0; i--) { + double sum = 0; + for (int j = i + 1; j < n; j++) { + sum += lu[i][j] * x[j]; + } + x[i] = (y[i] - sum) / lu[i][i]; + } + + return x; + } +} diff --git a/src/test/java/com/thealgorithms/maths/LUDecompositionTest.java b/src/test/java/com/thealgorithms/maths/LUDecompositionTest.java new file mode 100644 index 000000000000..daa7ebcea23e --- /dev/null +++ b/src/test/java/com/thealgorithms/maths/LUDecompositionTest.java @@ -0,0 +1,106 @@ +package com.thealgorithms.maths; + +import static org.junit.jupiter.api.Assertions.assertArrayEquals; +import static org.junit.jupiter.api.Assertions.assertThrows; + +import org.junit.jupiter.api.Test; + +public class LUDecompositionTest { + + private static final double DELTA = 1e-9; + + @Test + public void testDecomposeSimpleMatrix() { + double[][] matrix = {{2, 1, 1}, {4, 3, 3}, {8, 7, 9}}; + double[][] lu = LUDecomposition.decompose(matrix); + + double[][] lower = LUDecomposition.getLowerMatrix(lu); + double[][] upper = LUDecomposition.getUpperMatrix(lu); + + assertArrayEquals(new double[] {1, 1, 1}, new double[] {lower[0][0], lower[1][1], lower[2][2]}, DELTA); + + double[][] product = multiply(lower, upper); + assertArrayEquals(new double[] {2, 1, 1, 4, 3, 3, 8, 7, 9}, flatten(product), DELTA); + } + + @Test + public void testDecomposeTwoByTwo() { + double[][] matrix = {{1, 2}, {3, 4}}; + double[][] lu = LUDecomposition.decompose(matrix); + + double[][] lower = LUDecomposition.getLowerMatrix(lu); + double[][] upper = LUDecomposition.getUpperMatrix(lu); + + double[][] product = multiply(lower, upper); + assertArrayEquals(new double[] {1, 2, 3, 4}, flatten(product), DELTA); + } + + @Test + public void testDecomposeIdentityMatrix() { + double[][] matrix = {{1, 0}, {0, 1}}; + double[][] lu = LUDecomposition.decompose(matrix); + + double[][] lower = LUDecomposition.getLowerMatrix(lu); + double[][] upper = LUDecomposition.getUpperMatrix(lu); + + assertArrayEquals(new double[] {1, 0, 0, 1}, flatten(lower), DELTA); + assertArrayEquals(new double[] {1, 0, 0, 1}, flatten(upper), DELTA); + } + + @Test + public void testDecomposeNonSquareMatrixThrows() { + double[][] matrix = {{1, 2, 3}, {4, 5, 6}}; + assertThrows(IllegalArgumentException.class, () -> LUDecomposition.decompose(matrix)); + } + + @Test + public void testDecomposeSingularMatrixThrows() { + double[][] matrix = {{0, 1}, {1, 0}}; + assertThrows(ArithmeticException.class, () -> LUDecomposition.decompose(matrix)); + } + + @Test + public void testSolveLinearSystem() { + double[][] matrix = {{2, 1, 1}, {4, 3, 3}, {8, 7, 9}}; + double[] b = {8, 20, 46}; + double[][] lu = LUDecomposition.decompose(matrix); + double[] solution = LUDecomposition.solve(lu, b); + + assertArrayEquals(new double[] {1, 3, 3}, solution, DELTA); + } + + @Test + public void testSolveTwoByTwoSystem() { + double[][] matrix = {{2, 1}, {1, 3}}; + double[] b = {5, 7}; + double[][] lu = LUDecomposition.decompose(matrix); + double[] solution = LUDecomposition.solve(lu, b); + + assertArrayEquals(new double[] {1.6, 1.8}, solution, DELTA); + } + + private static double[][] multiply(double[][] a, double[][] b) { + int n = a.length; + double[][] result = new double[n][n]; + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + for (int k = 0; k < n; k++) { + result[i][j] += a[i][k] * b[k][j]; + } + } + } + return result; + } + + private static double[] flatten(double[][] matrix) { + int n = matrix.length; + double[] result = new double[n * n]; + int idx = 0; + for (double[] row : matrix) { + for (double val : row) { + result[idx++] = val; + } + } + return result; + } +}