「学习笔记」高斯消元

简单说：高斯消元就是我们初中学的解方程组时用的加减消元法和代入消元法，只是高斯这个人最后总结了一下

过程#

给定方程组

$$\left \{ \begin{aligned} 3x + 2y + z = 10 \quad &(1)\\ 5x + y + 6z = 25 \quad &(2)\\ 2x + 3y + 4z = 20 \quad &(3)\\ \end{aligned} \right.$$

我们用 $(2)$ 式来消去 $x$

$$\begin{aligned} & 5x + y + 6z = 25\\ & 5x = 25 - y - 6z\\ & x = 5 - \frac{1}{5}y - \frac{6}{5}z \end{aligned}$$

代入得

$$\left \{ \begin{aligned} & 3(5 - \frac{1}{5}y - \frac{6}{5}z) + 2y + z = 10 \\ & 5x + y + 6z = 25 \quad (忽略)\\ & 2(5 - \frac{1}{5}y - \frac{6}{5}z) + 3y + 4z = 20 \end{aligned} \right.$$

即

$$\left \{ \begin{aligned} & \frac{7}{5}y - \frac{13}{5}z = -5 \quad &(1)\\ & \frac{13}{5}y + \frac{8}{5}z = 10 \quad &(2) \end{aligned} \right.$$

我们现在再用 $(2)$ 式来消去 $y$ 得

$$\begin{aligned} & \frac{13}{5}y + \frac{8}{5}z = 10\\ & 13y + 8z = 50\\ & 13y = 50 - 8z\\ & y = \frac{50}{13} - \frac{8}{13} \end{aligned}$$

代入得

$$\left \{ \begin{aligned} & \frac{7}{5}(\frac{50}{13} - \frac{8}{13}) - \frac{13}{5}z = -5 \\ & \frac{13}{5}y + \frac{8}{5}z = 10 (忽略) \end{aligned} \right.$$

即

$$\begin{aligned} & \frac{45}{13}z = \frac{135}{13}\\ & z = 3 \end{aligned}$$

我们再将 $z = 3$ 代会得 $y = 2$，再带回得 $x = 1$

代码#

这里有应对 double 和 int 两种数据类型的代码。

应对 double 数据的#

double a[N][N];
cin >> n;
for (int i = 1; i <= n; ++ i) {
    for (int j = 1; j <= n; ++ j) {
        cin >> a[i][j];
    }
    cin >> a[i][n + 1];
}
// a[1][1] * x1 + a[1][2] * x2 + ... + a[1][n] * xn = a[1][n + 1];
// a[2][1] * x1 + a[2][2] * x2 + ... + a[2][n] * xn = a[2][n + 1];
// ...
// a[n][1] * x1 + a[n][2] * x2 + ... + a[n][m] * xn = a[n][n + 1];

void gauss() {
    for (int i = 1; i <= n; ++ i) { // 现在要把 xi 从第 i + 1 个方程到第 n 个方程消掉
        for (int j = i; j <= n; ++ j) {
            if (fabs(a[j][i]) > fabs(a[i][i])) {
                for (int k = 1; k <= n + 1; ++ k) {
                    swap(a[i][k], a[j][k]);
                }
            }
        }
        for (int j = i + 1; j <= n; ++ j) { // 把 xi 从第 j 个方程消掉
            double ratio = a[j][i] / a[i][i];
            for (int k = 1; k <= n + 1; ++ k) {
                a[j][k] -= a[i][k] * ratio;
            }
        }
    }
    for (int i = n; i >= 1; -- i) { // 解 x[i] 的值
        for (int j = i + 1; j <= n; ++ j) {
            a[i][n + 1] -= a[i][j] * x[j];
        }
        x[i] = a[i][n + 1] / a[i][i];
    }
}

应对 int 数据的#

int a[N][N];
cin >> n;
for (int i = 1; i <= n; ++ i) {
    for (int j = 1; j <= n; ++ j) {
        cin >> a[i][j];
    }
    cin >> a[i][n + 1];
}
// a[1][1] * x1 + a[1][2] * x2 + ... + a[1][n] * xn = a[1][n + 1];
// a[2][1] * x1 + a[2][2] * x2 + ... + a[2][n] * xn = a[2][n + 1];
// ...
// a[n][1] * x1 + a[n][2] * x2 + ... + a[n][m] * xn = a[n][n + 1];

void gauss() {
    for (int i = 1; i <= n; ++ i) { // 现在要把 xi 从第 i + 1 个方程到第 n 个方程消掉
        for (int j = i; j <= n; ++ j) {
            if (a[j][i] != 0) {
                for (int k = 1; k <= n + 1; ++ k) {
                    swap(a[i][k], a[j][k]);
                }
                break;
            }
        }
        for (int j = i + 1; j <= n; ++ j) { // 把 xi 从第 j 个方程消掉
            if (a[j][i] == 0)    continue;
            int l = a[i][i] / gcd(abs(a[i][i]), abs(a[j][i])) * a[j][i];
            int ratioi = l / a[i][i];
            int ratioj = l / a[j][i];
            for (int k = 1; k <= n + 1; ++ k) {
                a[j][k] = a[j][k] * ratioj - a[i][k] * ratioi;
            }
        }
    }
    for (int i = n; i >= 1; -- i) { // 解 x[i] 的值
        for (int j = i + 1; j <= n; ++ j) {
            a[i][n + 1] -= a[i][j] * x[j];
        }
        x[i] = a[i][n + 1] / a[i][i];
    }
}