c++高斯消元求解矩阵

// Gaussian elimination to solve Ax = b
vector<double> gaussian_elimination0(vector<vector<double>> A, vector<double> b) {
    int n = A.size();
    for (int i = 0; i < n; ++i) {
        A[i].push_back(b[i]);
    }

    for (int p = 0; p < n; ++p) {
        int max_row = p;
        for (int i = p + 1; i < n; ++i) {
            if (abs(A[i][p]) > abs(A[max_row][p])) {
                max_row = i;
            }
        }
        swap(A[p], A[max_row]);

        if (abs(A[p][p]) < 1e-10) {
            throw runtime_error("Matrix is singular");
        }

        for (int i = p + 1; i < n; ++i) {
            double alpha = A[i][p] / A[p][p];
            for (int j = p; j <= n; ++j) {
                A[i][j] -= alpha * A[p][j];
            }
        }
    }

    vector<double> x(n);
    for (int i = n - 1; i >= 0; --i) {
        double sum = 0.0;
        for (int j = i + 1; j < n; ++j) {
            sum += A[i][j] * x[j];
        }
        x[i] = (A[i][n] - sum) / A[i][i];
    }
    return x;
}

这个函数 gaussian_elimination0 使用高斯消元法（Gaussian Elimination）来解线性方程组 ( Ax = b )，其中 ( A ) 是一个 ( n \times n ) 的方阵，( b ) 是一个长度为 ( n ) 的向量，( x ) 是待求的解向量。

函数功能

• 输入：
  • A：一个 ( n \times n ) 的二维向量，表示系数矩阵。
  • b：一个长度为 ( n ) 的一维向量，表示方程组的右端项。
• 输出：
  • 一个长度为 ( n ) 的向量 ( x )，表示线性方程组 ( Ax = b ) 的解。
• 过程：
  1. 将向量 ( b ) 附加到矩阵 ( A ) 的最后一列，形成增广矩阵 ( [A | b] )。
  2. 通过高斯消元法将矩阵化为上三角形式：
     • 寻找主元（pivot），通过行交换确保主对角线元素绝对值最大。
     • 如果主元接近 0（小于 ( 1e-10 )），抛出异常（矩阵奇异）。
     • 对主元以下的行进行消元，消除主元列的非零元素。
  3. 使用回代（back substitution）从上三角矩阵求解 ( x )。
• 异常：
  • 如果矩阵 ( A ) 是奇异的（即不可逆），函数会抛出 runtime_error 异常。

代码逻辑

1. 增广矩阵：将 ( b ) 的元素追加到 ( A ) 的每一行，形成 ( [A | b] )。
2. 高斯消元：
   • 遍历每一列 ( p )，找到绝对值最大的主元行 ( max_row )。
   • 交换当前行 ( p ) 和主元行 ( max_row )。
   • 检查主元是否接近 0，若是则抛出异常。
   • 对 ( p ) 以下的行进行消元，使主元列的 ( p ) 以下元素为 0。
3. 回代求解：
   • 从最后一行开始，计算 ( x[i] ) 的值，利用已经求出的 ( x[j] )（( j > i )）。
   • 公式为 ( x[i] = (A[i][n] - \sum_{j=i+1}^{n-1} A[i][j] \cdot x[j]) / A[i][i] )。
4. 返回解向量 ( x )。

用途

• 解线性方程组，例如： [ \begin{cases} 2x_1 + x_2 = 5 \ x_1 + 3x_2 = 7 \end{cases} ] 其中 ( A = \begin{bmatrix} 2 & 1 \ 1 & 3 \end{bmatrix} ), ( b = \begin{bmatrix} 5 \ 7 \end{bmatrix} )，函数返回 ( x = \begin{bmatrix} x_1 \ x_2 \end{bmatrix} )。

注意事项

• 输入矩阵 ( A ) 必须是方阵，且不可奇异。
• 算法对数值稳定性敏感，主元选择（pivot selection）提高了稳定性，但仍可能因浮点误差导致问题。
• 如果矩阵规模较大或稀疏，可能需要其他方法（如 LU 分解或迭代法）来提高效率。