高斯消元

在算法竞赛中，高斯消元法常用于求解 \(n\) 元一次方程组。

一个包含 \(n\) 个变量和 \(m\) 个方程的线性方程组可以表示为：

\[\begin{cases} a_{11}x_1 + a_{12}x_2 + \dots + a_{1n}x_n = b_1 \\ a_{21}x_1 + a_{22}x_2 + \dots + a_{2n}x_n = b_2 \\ \dots \\ a_{m1}x_1 + a_{m2}x_2 + \dots + a_{mn}x_n = b_m \end{cases} \]

其对应的增广矩阵为：

\[\begin{bmatrix} a_{11} & a_{12} & \dots & a_{1n} & b_1 \\ a_{21} & a_{22} & \dots & a_{2n} & b_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{m1} & a_{m2} & \dots & a_{mn} & b_m \end{bmatrix} \]

利用以下三种变换（初等行变换），可以将增广矩阵转化为行阶梯形矩阵或简化行阶梯形矩阵：

交换两行：相当于交换两个方程的位置，方程组的解显然不变。
将某一行乘以一个非零常数：相当于在等式两边同时乘以一个非零数。
将某一行的若干倍加到另一行上：相当于将一个方程的倍数加到另一个方程上，利用加减消元法消去变量。

通过这些变换，将复杂的方程组转化为一个上三角形（行阶梯形）或对角形（简化行阶梯形）的方程组。例如，将方程组转化为简化行阶梯形矩阵后：

\[\begin{bmatrix} 1 & 0 & \dots & 0 & c_1 \\ 0 & 1 & \dots & 0 & c_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & c_n \end{bmatrix} \]

这直接对应了方程组的解：

\[\begin{cases} x_1 = c_1 \\ x_2 = c_2 \\ \dots \\ x_n = c_n \end{cases} \]

因此，高斯消元的过程就是通过一系列保解的变换，将方程组简化为可以直接读出解的形式。

以下是 Gauss-Jordan 消元法的算法流程：

枚举列 \(c\)（从第 1 列到第 \(n\) 列）
- 在当前未选中的行中，找到第 \(c\) 列上绝对值最大的行 \(r\)（为了数值稳定性）。
- 如果该列最大值为 0，说明此变量可能是自由变量，跳过该列。
- 将行 \(r\) 交换到当前未处理的首行。
- 将当前行的第 \(c\) 列系数化为 1。
- 利用当前行，将其他所有行（包括上方和下方）的第 \(c\) 列系数消为 0。
判断解的情况
- 如果所有变量都有对应的唯一行，则有唯一解。
- 如果出现 \(0 = 0\) 的行，说明有无穷多解（存在自由变量）。
- 如果出现 \(0 = d \ (d \ne 0)\) 的行，说明无解。

下面给出一个具体的例子：

\[\begin{cases} x_1 + 2x_2 - x_3 = -6 \\ 2x_1 + x_2 - 3x_3 = -9 \\ -x_1 - x_2 + 2x_3 = 7 \end{cases} \]

第一步：写出增广矩阵

\[\begin{bmatrix} 1 & 2 & -1 & -6 \\ 2 & 1 & -3 & -9 \\ -1 & -1 & 2 & 7 \end{bmatrix} \]

第二步：处理第 1 列

寻找第 1 列绝对值最大的行：第 2 行的 \(2\) 最大。交换第 1 行和第 2 行：

\[\begin{bmatrix} 2 & 1 & -3 & -9 \\ 1 & 2 & -1 & -6 \\ -1 & -1 & 2 & 7 \end{bmatrix} \]

将第 1 行系数化为 1（全行除以 2）：

\[\begin{bmatrix} 1 & 0.5 & -1.5 & -4.5 \\ 1 & 2 & -1 & -6 \\ -1 & -1 & 2 & 7 \end{bmatrix} \]

消去其他行的第 1 列：
- \(R_2 = R_2 - 1 \times R_1\)
- \(R_3 = R_3 - (-1) \times R_1\)

\[\begin{bmatrix} 1 & 0.5 & -1.5 & -4.5 \\ 0 & 1.5 & 0.5 & -1.5 \\ 0 & -0.5 & 0.5 & 2.5 \end{bmatrix} \]

第三步：处理第 2 列

第 2 行的 \(1.5\) 已经是最大，无需交换。
将第 2 行系数化为 1（除以 1.5）：

\[\begin{bmatrix} 1 & 0.5 & -1.5 & -4.5 \\ 0 & 1 & 1/3 & -1 \\ 0 & -0.5 & 0.5 & 2.5 \end{bmatrix} \]

消去其他行的第 2 列：
- \(R_1 = R_1 - 0.5 \times R_2 \Rightarrow R_1 = [1, 0, -5/3, -4]\)
- \(R_3 = R_3 - (-0.5) \times R_2 \Rightarrow R_3 = [0, 0, 2/3, 2]\)

\[\begin{bmatrix} 1 & 0 & -5/3 & -4 \\ 0 & 1 & 1/3 & -1 \\ 0 & 0 & 2/3 & 2 \end{bmatrix} \]

第四步：处理第 3 列

将第 3 行系数化为 1（除以 2/3）：

\[\begin{bmatrix} 1 & 0 & -5/3 & -4 \\ 0 & 1 & 1/3 & -1 \\ 0 & 0 & 1 & 3 \end{bmatrix} \]

消去其他行的第 3 列：
- \(R_1 = R_1 - (-5/3) \times R_3 \Rightarrow R_1 = [1, 0, 0, 1]\)
- \(R_2 = R_2 - (1/3) \times R_3 \Rightarrow R_2 = [0, 1, 0, -2]\)

\[\begin{bmatrix} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & -2 \\ 0 & 0 & 1 & 3 \end{bmatrix} \]

最终结果：
\(x_1 = 1, x_2 = -2, x_3 = 3\)。

double a[N][N]; // 增广矩阵 n * (n + 1)
double ans[N]; // 存储解
int gauss(int n) { // n: 变量个数
    int r = 0; // 当前处理的行
    for (int c = 0; c < n; ++c) { // 处理每一列
        int t = r;
        for (int i = r; i < n; ++i) // 找当前列绝对值最大的行（选主元）
            if (fabs(a[i][c]) > fabs(a[t][c])) t = i;
        if (fabs(a[t][c]) < EPS) continue; // 该列全为0，跳过
        if (t != r) {
            for (int i = 0; i <= n; i++) swap(a[t][i], a[r][i]]);
        }
        // 将当前行首位化为1
        for (int i = n; i >= c; --i) a[r][i] /= a[r][c];
        // 用当前行消去其他所有行的第 c 列
        for (int i = 0; i < n; ++i) {
            if (i != r && fabs(a[i][c]) > EPS) {
                double f = a[i][c];
                for (int j = c; j <= n; ++j)
                    a[i][j] -= a[r][j] * f;
            }
        }
        r++;
    }
    if (r < n) {
        for (int i = r; i < n; ++i) {
            if (fabs(a[i][n]) > EPS) return 1; // 0 = d 形式
        }
        return 2; // 0 = 0 形式
    }
    for (int i = 0; i < n; ++i) ans[i] = a[i][n];
    return 0;
}

由于有三层循环（枚举列、枚举行、消元），时间复杂度为 \(O(n^3)\)。

例题：P3389 【模板】高斯消元法

给定包含 \(n\) 个方程和 \(n\) 个未知数的线性方程组，要求求解该方程组。如果存在唯一解，输出各未知数的值（保留两位小数）；否则输出 No Solution。

模板题。

参考代码

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const int N = 105;
const double EPS = 1e-6;
// a[i][j] 存储增广矩阵，第 i 行第 j 列系数
// a[i][n] 存储等号右侧的常数 b[i]
double a[N][N];
bool solve(int n) {
    int r = 0; // 当前处理的行索引
    for (int c = 0; c < n; c++) { // 枚举每一列
        int t = r;
        // 寻找主元：在当前列及以下寻找绝对值最大的行，以保证数值稳定性
        for (int i = r; i < n; i++) {
            if (fabs(a[i][c]) > fabs(a[t][c])) t = i;
        }
        // 如果最大主元接近 0，说明该列所有元素都接近 0，此列无法消元，跳过
        if (fabs(a[t][c]) < EPS) continue;
        // 交换当前行与主元所在的行
        if (t != r) {
            for (int i = 0; i <= n; i++) swap(a[t][i], a[r][i]);
        }
        // 将当前行的主元系数化为 1
        for (int i = n; i >= c; i--) a[r][i] /= a[r][c];
        // 消去其他所有行在该列的系数（包括上方和下方的行）
        for (int i = 0; i < n; i++) {
            if (i != r && fabs(a[i][c]) > EPS) {
                double f = a[i][c];
                for (int j = c; j <= n; j++) {
                    a[i][j] -= f * a[r][j];
                }
            }
        }
        r++;
    }
    if (r < n) return -1;
    else return 0;
}
int main()
{
    int n; scanf("%d", &n);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j <= n; j++) {
            scanf("%lf", &a[i][j]);
        }
    }
    if (solve(n) != 0) {    
        printf("No Solution\n");
    } else {
        // 输出解，保留两位小数
        for (int i = 0; i < n; i++) {
            printf("%.2f\n", a[i][n]);
        }
    }
    return 0;
}

例题：P2455 [SDOI2006] 线性方程组

模板题，和上一题相比需要细分三种情况。

参考代码

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const int N = 55;
const double EPS = 1e-6;
double a[N][N]; // a[i][j] 存储增广矩阵
// 求解线性方程组，返回值：
// 0: 唯一解；-1: 无解；1: 无穷多解
int solve(int n) {
    int r = 0; // 当前处理的行
    for (int c = 0; c < n; c++) { // 枚举每一列（变量）
        int t = r;
        // 选主元：在当前行及以下寻找该列绝对值最大的行
        for (int i = r; i < n; i++) {
            if (fabs(a[i][c]) > fabs(a[t][c])) t = i;
        }    
        if (fabs(a[t][c]) < EPS) continue;
        // 交换当前行与主元行
        if (t != r) {
            for (int i = 0; i <= n; i++) swap(a[t][i], a[r][i]);
        }
        // 归一化：将主元系数化为 1
        double d = a[r][c];
        for (int i = c; i <= n; i++) a[r][i] /= d;
        // 消元：利用当前行消去其他所有行的第 c 列
        for (int i = 0; i < n; i++) {
            if (i != r && fabs(a[i][c]) > EPS) {
                double f = a[i][c];
                for (int j = c; j <= n; j++) {
                    a[i][j] -= f * a[r][j];
                }
            }
        }
        r++; // 处理完一个有效行
    }
    // 检查是否有解
    if (r < n) {
        for (int i = r; i < n; i++) {
            if (fabs(a[i][n]) > EPS) return -1; // 无解
        }
        return 1; // 无穷多解
    }
    return 0;
}
int main()
{
    int n; scanf("%d", &n);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j <= n; j++) {
            scanf("%lf", &a[i][j]);
        }
    }
    int res = solve(n);
    if (res == -1) printf("-1\n");
    else if (res == 1) printf("0\n");
    else {
        // 唯一解：输出每个变量的结果
        for (int i = 0; i < n; i++) {
            printf("x%d=%.2f\n", i + 1, a[i][n]);
        }
    }
    return 0;
}

例题：P4035 [JSOI2008] 球形空间产生器

给定 \(n\) 维空间球面上的 \(n+1\) 个点，要求求出该球面的球心坐标。

设球心坐标为 \(X = (x_1, x_2, \dots, x_n)\)，球面上第 \(i\) 个点的坐标为 \(P_i = (p_{i,1}, p_{i,2}, \dots, p_{i,n})\)。

根据球心的定义，球心到球面上所有点的距离相等（等于半径 \(R\)）。因此，对于每一个点 \(i \in [1,n+1]\)，都有 \(\sum \limits_{j=1}^n (p_{i,j} - x_j)^2 = R^2\)。

这是一个包含 \(n+1\) 个方程和 \(n+1\) 个未知数（\(x_1, \dots x_n\) 以及 \(R\)）的二次方程组，为了求解，需要对其进行线性化。

通过将相邻的两个方程相减（例如第 \(i+1\) 个方程减去第 \(i\) 个方程），可以消去二次项 \(x_j^2\) 和 \(R^2\)。

\[\sum_{j=1}^n (p_{i+1,j} - x_j)^2 - \sum_{j=1}^n (p_{i,j} - x_j)^2 = 0 \]

展开并简化：

\[\sum_{j=1}^n (p_{i+1,j}^2 - 2p_{i+1,j}x_j + x_j^2 - (p_{i,j}^2 - 2p_{i,j}x_j + x_j^2)) = 0 \]

\[\sum_{j=1}^n (p_{i+1,j}^2 - p_{i,j}^2 - 2x_j(p_{i+1,j} - p_{i,j})) = 0 \]

整理得到线性方程组形式：

\[\sum_{j=1}^n 2(p_{i+1,j} - p_{i,j})x_j = \sum_{j=1}^n (p_{i+1,j}^2 - p_{i,j}^2) \]

这样就得到了一个包含 \(n\) 个方程和 \(n\) 个未知数 \(x_j\) 的线性方程组 \(Ax=b\)，其中 \(A_{i,j} = 2(p_{i+1,j} - p_{i,j})\)，\(B_i = \sum \limits_{j=1}^n (p_{i+1,j}^2 - p_{i,j}^2)\)。

于是便可通过高斯消元法求解。

参考代码

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
double p[15][15], m[15][15];
int main()
{
    int n; scanf("%d", &n);
    // 读取 n+1 个点的坐标
    for (int i = 1; i <= n + 1; i++) {
        for (int j = 1; j <= n; j++) {
            scanf("%lf", &p[i][j]);
        }
    }
    // 构建线性方程组 Ax = B
    // 每一行 i 是由第 i+1 个点和第 i 个点的距离公式相减得到的
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            m[i][j] = 2.0 * (p[i + 1][j] - p[i][j]);
            m[i][n + 1] += p[i + 1][j] * p[i + 1][j] - p[i][j] * p[i][j];
        }
    }
    // 高斯-约旦消元法
    for (int i = 1; i <= n; i++) {
        // 主元选择
        int pe = i;
        for (int j = i + 1; j <= n; j++) {
            if (fabs(m[j][i]) > fabs(m[pe][i])) {
                pe = j;
            }
        }
        // 交换行
        if (pe != i) {
            for (int j = 1; j <= n + 1; j++) {
                swap(m[i][j], m[pe][j]);
            }
        } 
        // 归一化当前行
        double d = m[i][i];
        for (int j = i; j <= n + 1; j++) {
            m[i][j] /= d;
        }
        // 消去其他行的当前列
        for (int j = 1; j <= n; j++) {
            if (i != j) {
                double f = m[j][i];
                for (int k = i; k <= n + 1; k++) {
                    m[j][k] -= f * m[i][k];
                }
            }
        }
    }
    // 输出结果，精确到小数点后 3 位
    for (int i = 1; i <= n; i++) {
        printf("%.3f ", m[i][n + 1]);
    }
    return 0;
}

posted @ 2026-03-14 10:23 RonChen 阅读(1) 评论(0) 收藏举报

刷新页面返回顶部

RonChen

高斯消元

例题：P3389 【模板】高斯消元法

例题：P2455 [SDOI2006] 线性方程组

例题：P4035 [JSOI2008] 球形空间产生器

公告