【数学】高斯-约旦消元法

给定 \(n\) 元一次方程组

\[\begin{cases} a_{1,1}x_1+a_{1,2}x_2+\cdots+a_{1,n}x_n=b_1\\ a_{2,1}x_1+a_{2,2}x_2+\cdots+a_{2,n}x_n=b_2\\ \cdots\\ a_{n,1}x_1+a_{n,2}x_2+\cdots+a_{n,n}x_n=b_n\\ \end{cases} \]

请求出方程组的解的情况:

  • 无解;

  • 无穷多解;

  • 唯一解。

对于这样的问题,我们可以使用 高斯消元法 进行求解,当然高斯消元法有一个回代的过程,代码略长,而且精度较低。

所以我们隆重推出 高斯-约旦消元法 !!!

回顾一下我们是怎么手算的,一般用的都是 加减消元法,普通高斯和高斯-约旦用的都是加减消元。

在此之前,我们需要了解一下矩阵初等变换。

在线性代数中,矩阵初等行变换 是指以下三种变换类型 :

  1. 交换矩阵的两行;

  2. 用一个非零数 \(k\) 乘矩阵的某一行所有元素;

  3. 把矩阵的某一行所有元素乘以一个数 \(k\) 后加到另一行对应的同一列的元素上;

类似地,把以上的 改为 便得到 矩阵初等列变换 的定义。

矩阵初等行变换与初等列变换合称为 矩阵初等变换

若矩阵 \(A\) 经过有限次的初等行变换变为矩阵 \(B\),则矩阵 \(A\) 与矩阵 \(B\) 行等价;若矩阵 \(A\) 经过有限次的初等列变换变为矩阵 \(B\),则矩阵 \(A\) 与矩阵 \(B\) 列等价;若矩阵 \(A\) 经过有限次的初等变换变为矩阵 \(B\),则矩阵 \(A\) 与矩阵 \(B\) 等价

当然列的用不着

首先有一个由系数构成的 \(n\times n\) 的矩阵

\[\begin{bmatrix} a_{1,1}&a_{1,2}&\cdots&a_{1,n}\\ a_{2,1}&a_{2,2}&\cdots&a_{2,n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{n,1}&a_{n,2}&\cdots&a_{n,n}\\ \end{bmatrix} \]

然后是一个由常数构成的 \(n\times 1\) 的列向量

\[\begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_n \end{bmatrix} \]

把它们放在一起构成一个 \(n\times(n+1)\) 的增广矩阵:

\[\begin{bmatrix} a_{1,1}&a_{1,2}&\cdots&a_{1,n}&\mid&b_1\\ a_{2,1}&a_{2,2}&\cdots&a_{2,n}&\mid&b_2\\ \vdots&\vdots&\ddots&\vdots&\mid&\vdots\\ a_{n,1}&a_{n,2}&\cdots&a_{n,n}&\mid&b_n\\ \end{bmatrix} \]

我们遍历每一列,对于第 \(i\) 列选出最大的 未处理过的行 上的数作为主元,意味着加减消元后除了主元这一行外其他行的第 \(i\) 列都为 \(0\)。

选最大的作为主元的原因是避免精度误差。

然后就是加减消元了。

举个例子

\[\begin{cases} 2x-y+z=1\\ 4x+y-z=5\\ x+y+z=0 \end{cases} \]

增广矩阵就是

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

先是第 \(1\) 列,选 \(4\) 为主元。

\(4\) 在第 \(2\) 行,将第 \(2\) 行与第 \(1\) 行交换。

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

对于第 \(2\) 行,第一列上的 \(2\) 是 \(4\) 的 \(\dfrac{1}{2}\)。

\[\begin{bmatrix} 4&1&-1&\mid&5\\ 2-4\times\dfrac{1}{2}&-1-1\times\dfrac{1}{2}&1-(-1)\times\dfrac{1}{2}&\mid&1-5\times\dfrac{1}{2}\\ 1&1&1&\mid&0 \end{bmatrix} \]

化简得

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

对第 \(2\) 行的处理同理

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

现在到了第 \(2\) 列,未处理过的是 \(2,3\) 行,选最大的 \(\dfrac{3}{2}\) 为主元。

将第 \(2\) 行与第 \(2\) 行交换(相当于没换)

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

消元得

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

第 \(3\) 列,未处理过的是第 \(3\) 行,选 \(\dfrac{1}{2}\) 作主元。

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

这样就把原来的矩阵通过初等变换,使得系数矩阵只有主对角线位置的元素非零,即一个对角矩阵。

上面那个矩阵的意思是

\[\begin{cases} 4x=4\\ \dfrac{3}{2}y=0\\ \dfrac{1}{2}z=-\dfrac{1}{2} \end{cases} \]

所以再把系数除过去就得到

\[\begin{cases} x=1\\ y=0\\ z=-1 \end{cases} \]

请自行检验。

时间复杂度为 \(\mathcal{O}(n^3)\)。

当然方程还可能是无解或无穷多解。

考虑一个一元一次方程 \(ax=b\) 的解的情况:

  1. 无解:\(a=0,b\ne0\);
  2. 无穷多解:\(a=b=0\);
  3. 唯一解:\(a\ne0\)。
  4. 所以当发现某一列的主元为 \(0\) 时,因为主元是最大的,所以意味着这一列全都是 \(0\),那么要么无解,要么有无穷多解。

P3389 【模板】高斯消元法

\(\text{Code}\)

#include <iostream>
#include <cstdio>
#include <cmath>
typedef double db;
using namespace std;

const int MAXN = 105;

int n;
db a[MAXN][MAXN];

bool Gauss_Jordan()
{
	for (int i = 1; i <= n; i++)
	{
		int mx = i;
		for (int j = i + 1; j <= n; j++)
		{
			if (fabs(a[j][i]) < fabs(a[mx][i]))
			{
				mx = j;
			}
		}
		if (mx != i)
		{
			swap(a[i], a[mx]);
		}
		if (!a[i][i])
		{
			return false;
		}
		for (int j = 1; j <= n; j++)
		{
			if (i != j)
			{
				db val = a[j][i] / a[i][i];
				for (int k = i + 1; k <= n + 1; k++)
				{
					a[j][k] -= a[i][k] * val;
				}
			}
		}
	}
	return true;
}

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_Jordan())
	{
		for (int i = 1; i <= n; i++)
		{
			printf("%.2lf\n", a[i][n + 1] / a[i][i]);
		}
	}
	else
	{
		puts("No Solution");
	}
	return 0;
}

关于判断无解和无穷多解

P2455 [SDOI2006]线性方程组

这下要具体到到底是无解还是无穷多解了。

这就是高斯-约旦的一个缺点:相比于普通高斯,它更难判断无解和无穷多解。

但是也是可以处理的。

具体地,我们用 \(r\) 来记录当前行,如果主元为 \(0\),那么 \(\operatorname{continue}\) 到下一列,但 \(r\) 不变;否则消元后令 \(r\gets r+1\)。

遍历完所有列后:

  • 若 \(r=n\),说明有唯一解;
  • 若 \(r<n\),说明第 \(r+1\sim n\) 行系数矩阵全都是 \(0\),若列向量全是 \(0\),说明有无穷多解,否则无解。

\(\text{Code}\)

#include <iostream>
#include <cstdio>
#include <cmath>
typedef double db;
using namespace std;

const int MAXN = 55;

int n;
db a[MAXN][MAXN], ans[MAXN];

int Gauss_Jordan()
{
	int r = 1;
	for (int i = 1; i <= n; i++)
	{
		int mx = r;
		for (int j = r + 1; j <= n; j++)
		{
			if (fabs(a[j][i]) > fabs(a[mx][i]))
			{
				mx = j;
			}
		}
		if (mx != r)
		{
			swap(a[r], a[mx]);
		}
		if (!a[r][i])
		{
			continue;
		}
		for (int j = 1; j <= n; j++)
		{
			if (j != r)
			{
				db val = a[j][i] / a[r][i];
				for (int k = i + 1; k <= n + 1; k++)
				{
					a[j][k] -= a[r][k] * val;
				}
			}
		}
		r++;
	}
	if (--r < n)
	{
		for (int i = r + 1; i <= n; i++)
		{
			if (a[i][n + 1])
			{
				return -1;
			}
		}
		return 0;
	}
	for (int i = 1; i <= n; i++)
	{
		ans[i] = a[i][n + 1] / a[i][i];
		if (!ans[i])
		{
			ans[i] = 0;
		}
	}
	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);
		}
	}
	int res = Gauss_Jordan();
	if (res != 1)
	{
		printf("%d\n", res);
	}
	else
	{
		for (int i = 1; i <= n; i++)
		{
			printf("x%d=%.2lf\n", i, ans[i]);
		}
	}
	return 0;
}
上一篇:【学军NOIP开放题2-A】塔纳(DP)


下一篇:【做题记录】图论杂题