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

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

概念

高斯消元是一种求解线性方程组的方法。线性方程组是由 \(M\) 个 \(N\) 元一次方程共同构成。线性方程组可以写成 \(M\) 行 \(N\) 列的矩阵，加上等号右边的常数，即可写成一个 \(M\) 行 \(N\) 列的增广矩阵。

如：

求解线性方程组可以进行对增广矩阵的三类操作：

1. 用一个非零的数乘某一行
2. 将其中一行的若干倍加到另一行上
3. 交换两行的位置

这三类操作称为矩阵的“初等行变换”，我们可以用若干次初等行变换求解方程组，如：

最后即可得到矩阵的“上三角矩阵：

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

再简化一番，可以得到一个”对角矩阵”：

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

通过初等行变换把增广矩阵简化为上三角矩阵或对角矩阵的过程就是高斯消元。高斯消元的思想是，对于每一个未知量 \(x_i\)，找到一个 \(x_i\) 系数非 0，\(x_{1..i-1}\) 系数全为0的方程，通过初等行变换消去其他方程的 \(x_i\) 项。而高斯消元法与高斯-约旦消元法不同的地方在于高斯消元法只消去后面方程组的 \(x_i\) 项，而高斯-约旦消元法将其他所有的 \(x_i\) 项都消去。所以高斯消元法最终得到的是上三角矩阵，而高斯-约旦消元法最终得到的则是对角矩阵。

实现

例题：P3389 P2455

高斯消元

首先找到一个第一列的数不为0的行（一般找第一列的数最大的行）（如果都为0就跳过当前步）

然后用它的第一列的数将下面行当前列的值化为0，变换过的初等矩阵与原矩阵等价，化为方程后依然成立

本矩阵第一列的数最大的为第三行就把第三行与第一行交换

然后下面行的当前列消去，如图

除了最后一列外，每一列都如此，最后得到上三角矩阵如图

这样我们很容易算出\(x3\)的值，再用\(x3\)的值算出\(x2\)，\(x1\)的值

Code

int gauss(int n) {
	int r, w = 0;
	for (int i = 0; w < n && i < n; i++, w++) { //枚举列（每一项）
		r = w;
		for (int j = w + 1; j < n; j++) 
			if (fabs(a[j][i]) > fabs(a[r][i])) //找到最大的系数
				r = j;
		if (fabs(a[r][i]) < eps) {      //是0跳过
			w--;
			continue;
		}
		if (w != r) {
			for (int j = 0; j <= n; j++)   //换到第一行
				swap(a[r][j], a[w][j]);
		}
		for (int k = w + 1; k < n; k++) {
			double p = a[k][i] / a[w][i];
			for (int j = w; j <= n; j++)     
				a[k][j] -= p * a[w][j];
		}	
	}
	return w; //有效方程数量+1
}

高斯-约旦消元

高斯-约旦消元与高斯消元大体差不多，但消元时上面计算过的行也要消去当前列，最后得到的是对角矩阵而不是上三角矩阵。

Code

int gauss_jordan(int n) {
	int r, w = 0;
	for (int i = 0; w < n && i < n; i++, w++) {
		r = w;
		for (int j = w + 1; j < n; j++) 
			if (fabs(a[j][i]) > fabs(a[r][i])) 
				r = j;
		if (fabs(a[r][i]) < eps) {
			w--;
			continue;
		}
		if (w != r) {
			for (int j = 0; j <= n; j++)
				swap(a[r][j], a[w][j]);
		}
		for (int k = 0; k < n; k++) {   //消去所有方程的xi项
			if (k == w) continue;
			double p = a[k][i] / a[w][i];
			for (int j = n; j >= w; j--) 
				a[k][j] -= p * a[w][j]; 
		}
	}
	return w;
}

判断无解与无穷解

注意到我们在函数最后返回了 \(w\)，而 \(w-1\) 也就是有效方程个数。（实际上一般增广矩阵从 0 存的话 \(w\) 才是有效方程个数，要从 0 开始记）

如果有效方程个数小于 \(n\)，那么就无解了！

然后判断无穷解和无解的情况：

• 如果剩下的方程左边=右边=0，那么就是无穷解答
• 如果剩下方程左边=0，右边!=0，那么就是无解

int d = gauss_jordan(n);
if(d < n){
	while(d < n)
		if(a[d++][n] != 0) return puts("-1") && 0;
	return putchar('0') && 0;
}