高斯消元与高斯-约旦消元

题目1

洛谷 P3389 【模板】高斯消元法
总的来说，就是求解一个 n n n 元一次方程组。

高斯消元

思路：

首先把所有系数看成一个矩阵：
[ a 1 , 1 a 1 , 2 a 1 , 3 a 1 , 4 a 1 , 5 a 2 , 1 a 2 , 2 a 2 , 3 a 2 , 4 a 2 , 5 a 3 , 1 a 3 , 2 a 3 , 3 a 3 , 4 a 3 , 5 a 4 , 1 a 4 , 2 a 4 , 3 a 4 , 4 a 4 , 5 a 5 , 1 a 5 , 2 a 5 , 3 a 5 , 4 a 5 , 5 ] \begin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4} & a_{1,5} \\ \\ a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & a_{2,5} \\ \\ a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & a_{3,5} \\ \\ a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4} & a_{4,5} \\ \\ a_{5,1} & a_{5,2} & a_{5,3} & a_{5,4} & a_{5,5} \\ \end{bmatrix} ​a1,1​a2,1​a3,1​a4,1​a5,1​​a1,2​a2,2​a3,2​a4,2​a5,2​​a1,3​a2,3​a3,3​a4,3​a5,3​​a1,4​a2,4​a3,4​a4,4​a5,4​​a1,5​a2,5​a3,5​a4,5​a5,5​​ ​
然后将所有方程等号右边的常数放到最后一列：
[ a 1 , 1 a 1 , 2 a 1 , 3 a 1 , 4 a 1 , 5 b 1 a 2 , 1 a 2 , 2 a 2 , 3 a 2 , 4 a 2 , 5 b 2 a 3 , 1 a 3 , 2 a 3 , 3 a 3 , 4 a 3 , 5 b 3 a 4 , 1 a 4 , 2 a 4 , 3 a 4 , 4 a 4 , 5 b 4 a 5 , 1 a 5 , 2 a 5 , 3 a 5 , 4 a 5 , 5 b 5 ] \begin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4} & a_{1,5} & b_1 \\ \\ a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & a_{2,5} & b_2 \\ \\ a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & a_{3,5} & b_3 \\ \\ a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4} & a_{4,5} & b_4 \\ \\ a_{5,1} & a_{5,2} & a_{5,3} & a_{5,4} & a_{5,5} & b_5 \\ \end{bmatrix} ​a1,1​a2,1​a3,1​a4,1​a5,1​​a1,2​a2,2​a3,2​a4,2​a5,2​​a1,3​a2,3​a3,3​a4,3​a5,3​​a1,4​a2,4​a3,4​a4,4​a5,4​​a1,5​a2,5​a3,5​a4,5​a5,5​​b1​b2​b3​b4​b5​​ ​
这样就构成了增广矩阵。

步骤：

1. 枚举每一行 i i i，并将这一行的第 i i i 个未知数作为这一行的主元 x i x_i xi​。 如果这个未知数的系数为 0 0 0，那就从这行开始，向下几行寻找一个该未知数系数不为 0 0 0 的方程，找到了就将这一行与当前行交换（找不到的情况就是 N o   S o l u t i o n No \space Solution No Solution，为什么后面会说）。
2. 下来再将该主元的系数化为 1 1 1，用这个主元去消下面几个方程的相同未知数的系数，即：把矩阵中与该主元在同一列且在该主元下面的系数消为 0 0 0。
3. 让后我们就可以从下往上带入求解。

完成 1 、 2 1、2 1、2 步后，我们就会得到一个除最后一个“常数列”外的上三角矩阵，即：除了最后一列，以这个矩阵的 “左上—右下” 对角线为界，其左下角（不包括对角线）都为 0 0 0。
如下图：
[ a 1 , 1 ′ a 1 , 2 ′ a 1 , 3 ′ a 1 , 4 ′ a 1 , 5 ′ b 1 ′ 0 a 2 , 2 ′ a 2 , 3 ′ a 2 , 4 ′ a 2 , 5 ′ b 2 ′ 0 0 a 3 , 3 ′ a 3 , 4 ′ a 3 , 5 ′ b 3 ′ 0 0 0 a 4 , 4 ′ a 4 , 5 ′ b 4 ′ 0 0 0 0 a 5 , 5 ′ b 5 ′ ] \begin{bmatrix} a_{1,1}' & a_{1,2}' & a_{1,3}' & a_{1,4}' & a_{1,5}' & b_1' \\ \\ 0 & a_{2,2}' & a_{2,3}' & a_{2,4}' & a_{2,5}' & b_2' \\ \\ 0 & 0 & a_{3,3}' & a_{3,4}' & a_{3,5}' & b_3' \\ \\ 0 & 0 & 0 & a_{4,4}' & a_{4,5}' & b_4' \\ \\ 0 & 0 & 0 & 0 & a_{5,5}' & b_5' \\ \end{bmatrix} ​a1,1′​0000​a1,2′​a2,2′​000​a1,3′​a2,3′​a3,3′​00​a1,4′​a2,4′​a3,4′​a4,4′​0​a1,5′​a2,5′​a3,5′​a4,5′​a5,5′​​b1′​b2′​b3′​b4′​b5′​​ ​
最后，从下往上带入求解每一元。

代码

#include <bits/stdc++.h>

#define mkpr make_pair
#define fir first
#define sec second

using namespace std;

typedef long long ll;
typedef pair<int, int> pii;

const int maxn = 1e2 + 7;
const double eps = 1e-9;

int n;
double a[maxn][maxn];
bool gauss() {
	for (int i = 1; i <= n; ++i) {
		int r = i;
		for (int j = i; j <= n; ++j)
		    if (abs(a[r][i]) < abs(a[j][i])) // 这里找该元系数最大的方程可以减小误差
		        r = j;
	    if (r != i) swap(a[r], a[i]);
	    if (abs(a[i][i]) < eps) return 0;
	    
	    for (int j = n + 1; j >= i; --j)
	        a[i][j] /= a[i][i];
	    for (int j = i + 1; j <= n; ++j)
	    	for (int k = n + 1; k >= i; --k)
	    		a[j][k] -= a[i][k] * a[j][i];
	}
	for (int i = n - 1; i >= 1; --i)
		for (int j = i + 1; j <= n; ++j)
			a[i][n + 1] -= a[i][j] * a[j][n + 1];
	return 1;
}
int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n + 1; ++j)
            scanf("%lf", &a[i][j]);
            
    if (gauss()) {
    	for (int i = 1; i <= n; ++i) 
		    printf("%.2lf\n", a[i][n + 1]);
	} else {
		printf("No Solution\n");
	}
	return 0;
}

高斯-约旦消元

思路

同上，将系数、常数视为一个矩阵。不同的是如何消元。

步骤：

1. 依旧是枚举 i i i 行，循环找到主元 x i x_i xi​ 的系数不为 0 0 0 的方程。
2. 再次拿 x i x_i xi​ 去消其他所有相同的元，“消去所有” 指的就是矩阵中与 a i , i a_{i,i} ai,i​（即此方程中 x i x_i xi​ 的系数） 在同一列的、除 a i , i a_{i,i} ai,i​ 自己的其他系数化为 0 0 0。

最后就得到了一条只有对角线的系数矩阵：
[ a 1 , 1 ′ 0 0 0 0 b 1 ′ 0 a 2 , 2 ′ 0 0 0 b 2 ′ 0 0 a 3 , 3 ′ 0 0 b 3 ′ 0 0 0 a 4 , 4 ′ 0 b 4 ′ 0 0 0 0 a 5 , 5 ′ b 5 ′ ] \begin{bmatrix} a_{1,1}' & 0 &0 & 0 & 0 & b_1' \\ \\ 0 & a_{2,2}' & 0 & 0 & 0 & b_2' \\ \\ 0 & 0 & a_{3,3}' & 0 & 0 & b_3' \\ \\ 0 & 0 & 0 & a_{4,4}' & 0 & b_4' \\ \\ 0 & 0 & 0 & 0 & a_{5,5}' & b_5' \\ \end{bmatrix} ​a1,1′​0000​0a2,2′​000​00a3,3′​00​000a4,4′​0​0000a5,5′​​b1′​b2′​b3′​b4′​b5′​​ ​
剩余系数以此化为 1 1 1 就能得到解。

代码

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int maxn = 1000 + 7;
const int inf  = 0x3f3f3f3f;
const double eps = 1e-4;

int n;
double a[maxn][maxn];
int guass_jordan() {
	for (int i = 1; i <= n; ++i) {
		int maxi = i;
		for (int j = i + 1; j <= n; ++j)
		    if (abs(a[maxi][i]) < abs(a[j][i])) 
			    maxi = j;
		if (maxi != i) swap(a[i], a[maxi]);
		if (abs(a[i][i]) <= eps) return 0;
		
	    for (int j = 1; j <= n; ++j) {
	    	if (j == i) continue;
	    	double mul = a[j][i] / a[i][i];
	    	for (int k = i; k <= n + 1; ++k)
	    	    a[j][k] -= mul * a[i][k];
		}
    }
    for (int i = 1; i <= n; ++i)
	    a[i][n + 1] /= a[i][i];
    return 1;
}
int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n + 1; ++j)
            scanf("%lf", &a[i][j]);
    
    if (guass_jordan()) {
    	for (int i = 1; i <= n; ++i)
    	    printf("%.2lf\n", a[i][n + 1]);
	}
    else printf("No Solution\n");
	return 0;
}

关于为什么系数为 0 0 0 就是无解或无限解

很显然我们可以在消元过程中遇到没有一个系数主元 x i x_i xi​ 就把它跳过，因为它不会影响剩下若干元的求解。
那么在高斯消元中，最后带入求解时，就会有一行的系数全是 0 0 0。此时再看等号右边，如果为 0 0 0，那就有无数解；如果不为 0 0 0，那就无解。

可以类比一元一次方程 k x = b kx=b kx=b，
k = 0 k=0 k=0 时，若 b = 0 b = 0 b=0 就有无数解，若 b ≠ 0 b \neq 0 b=0 就无解。

题目2

洛谷 P2455 [SDOI2006] 线性方程组

这题让具体判断是有无限解还是无解的情况。
根据上面的分析，当一个元的系数为 0 0 0 时，等号右边等于 0 0 0 就是有无数解，不等于 0 0 0 就是无解。

大体思路还是一样的，不过在寻找当前元 x i x_i xi​ 的系数非零方程时，不应该从第 i i i 个开始找，而应该是从第 1 1 1 个开始找。且如果找到的方程是在第 i i i 个方程前面（假如是 j j j），那么就要判断那个方程的主元 x j x_j xj​ 的系数是否为零。为零就用，不为零就不能用，因为如果不为零，主元 x j x_j xj​ 和 x i x_i xi​ 就会出现在同一个方程且系数都不为 0 0 0，此时就没有办法进行消元。
至于为什么从 1 1 1 开始，可以看看这篇题解。

注意： 无解的优先级比无数解高。即：在出现了若干行系数都为 0 0 0，而等号右边有的等于零，有有的不等于零的时，输出无解，而不是无数解。

代码

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int maxn = 50 + 7;
const int inf  = 0x3f3f3f3f;
const double eps = 1e-8;

int n;
double a[maxn][maxn];
int guass_jordan() {
	int newr = 1;  // 新的一行 
	for (int i = 1; i <= n; ++i) {  // 这里不再是枚举行，而是列 
		int maxi = i;
		for (int j = 1; j <= n; ++j) { // 这里枚举行 
		    if (abs(a[j][j]) > eps && j < i) continue;
		    if (abs(a[maxi][i]) < abs(a[j][i])) maxi = j;
	    }
		if (maxi != i) swap(a[i], a[maxi]);
		if (abs(a[i][i]) <= eps) continue;
	    for (int j = 1; j <= n; ++j) {
	    	if (j == i) continue;
	    	double mul = a[j][i] / a[i][i];
	    	for (int k = i; k <= n + 1; ++k)
	    	    a[j][k] -= mul * a[i][k];
		}
    }
	int res = 1;
	for (int i = 1; i <= n; ++i) {
		if (abs(a[i][i]) <= eps) {
			if (abs(a[i][n + 1]) > eps) res = -1;
		    else if (res != -1) res = 0;
		}
	}
    return res;
}
int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n + 1; ++j)
            scanf("%lf", &a[i][j]);
    int ans = guass_jordan();
    if (ans != 1) printf("%d\n", ans);
    else {
    	for (int i = 1; i <= n; ++i)
    	    printf("x%d=%.2lf\n", i, a[i][n + 1] / a[i][i]);
	}
	return 0;
}