Java 实现线性方程组求解：LU 分解与部分主元法

Java 实现线性方程组求解：LU 分解与部分主元法

在数值计算中，求解线性方程组 $Ax=b$ 是基础且高频的任务。虽然可以直接调用库函数，但理解底层的 LU 分解（LUP）有助于处理更复杂的数值问题。这里分享一个基于 Java 的完整实现，重点展示了如何通过部分主元法（Partial Pivoting）来避免除零错误并提高精度。

核心思路

我们将系数矩阵 $A$ 分解为一个下三角矩阵 $L$ 和一个上三角矩阵 $U$，同时引入置换向量 $p$ 记录行交换情况。这样原方程转化为 $PAx = Pb$，即 $LUx = Pb$。求解分为两步：先向前代入求中间变量 $Y$，再向后回代求最终解 $X$。

代码实现

package com.data.struct;

public class LinearEquation {
    // 系数矩阵 A
    private double[][] A = new double[][]{
        {1, 2, 0},
        {3, 4, 4},
        {5, 6, 3}
    };
    // 常数向量 B
    private double[] B = new double[]{3, 7, 8};
    
    // 存储分解后的 L 和 U 矩阵
    private double[][] L = new double[3][3];
    private double[][] U = new double[3][3];
    // 行交换记录
    private int[] p = new int[3];
    // 解向量 X
    private double[] X = new double[3];

    public void solve() throws Exception {
        lupDecomposition();
        
        // 1. 前代过程：求解 Ly = Pb
        double[] Y = new double[3];
        for (int i = 0; i < 3; i++) {
            Y[i] = B[p[i]];
            double sub = 0;
            for (int j = 0; j <= i - 1; j++) {
                sub += L[i][j] * Y[j];
            }
            Y[i] = Y[i] - sub;
        }

        // 2. 回代过程：求解 Ux = Y
        for (int i = 2; i >= 0; i--) {
            double sub = 0;
            for (int j = i + 1; j < 3; j++) {
                sub += U[i][j] * X[j];
            }
            X[i] = (Y[i] - sub) / U[i][i];
        }
    }

    // LUP 分解核心逻辑
    private void lupDecomposition() throws Exception {
        // 初始化置换向量
        for (int i = 0; i < 3; i++) {
            p[i] = i;
        }

        for (int k = 0; k < 3; k++) {
            // 选主元：寻找当前列绝对值最大的行
            double s = 0;
            int k2 = k;
            for (int i = k; i < 3; i++) {
                if (Math.abs(A[i][k]) > s) {
                    s = Math.abs(A[i][k]);
                    k2 = i;
                }
            }
            
            // 如果主元为 0，说明矩阵奇异
            if (s == 0) {
                throw new Exception("error");
            }

            // 交换行记录
            int tmp = p[k];
            p[k] = p[k2];
            p[k2] = tmp;

            // 交换实际矩阵行
            for (int i = 0; i < 3; i++) {
                double tmp2 = A[k][i];
                A[k][i] = A[k2][i];
                A[k2][i] = tmp2;
            }

            // 高斯消元更新 A 为 L 和 U 的混合形式
            for (int i = k + 1; i < 3; i++) {
                A[i][k] = A[i][k] / A[k][k];
                for (int j = k + 1; j < 3; j++) {
                    A[i][j] = A[i][j] - A[i][k] * A[k][j];
                }
            }
        }

        // 从修改后的 A 中提取 L 和 U
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j <= i; j++) {
                if (i == j) {
                    L[i][j] = 1;
                } else {
                    L[i][j] = A[i][j];
                }
            }
            for (int j = i; j < 3; j++) {
                U[i][j] = A[i][j];
            }
        }
    }

    public void printL() {
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                System.out.print(L[i][j] + " ");
            }
            System.out.println();
        }
    }

    public void printU() {
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                System.out.print(U[i][j] + " ");
            }
            System.out.println();
        }
    }

    public void print() {
        System.out.println("Solution X:");
        for (int i = 0; i < 3; i++) {
            System.out.println("X[" + i + "] = " + X[i]);
        }
    }

    public static void main(String[] args) {
        try {
            LinearEquation solver = new LinearEquation();
            solver.solve();
            System.out.println("Matrix L:");
            solver.printL();
            System.out.println("Matrix U:");
            solver.printU();
            solver.print();
        } catch (Exception e) {
            e.printStackTrace();
        }
    }
}

关键点解析

1. 部分主元法：在每一列消元前，都会查找该列下方绝对值最大的元素作为主元。这能有效防止除以极小值导致的数值溢出或精度丢失。
2. 原地分解：为了节省空间，代码直接在原矩阵 A 上进行操作，利用下三角区域存 L（不含对角线），上三角区域存 U。
3. 异常处理：当检测到主元为 0 时抛出异常，这意味着矩阵不可逆，无法求解唯一解。

总结

这个实现虽然针对的是 3x3 矩阵，但逻辑适用于任意阶数。在实际工程中，如果涉及大规模计算，建议直接使用 Apache Commons Math 等成熟库，但在面试或算法学习场景下，手写一遍能帮你彻底搞懂数值线性代数的门道。