【数学】高斯-约旦消元法
给定 \(n\) 元一次方程组
请求出方程组的解的情况：
-
无解；
-
无穷多解；
-
唯一解。
对于这样的问题，我们可以使用 高斯消元法 进行求解，当然高斯消元法有一个回代的过程，代码略长，而且精度较低。
所以我们隆重推出 高斯-约旦消元法 ！！！
回顾一下我们是怎么手算的，一般用的都是 加减消元法，普通高斯和高斯-约旦用的都是加减消元。
在此之前，我们需要了解一下矩阵初等变换。
在线性代数中，矩阵初等行变换 是指以下三种变换类型 ：
-
交换矩阵的两行；
-
用一个非零数 \(k\) 乘矩阵的某一行所有元素；
-
把矩阵的某一行所有元素乘以一个数 \(k\) 后加到另一行对应的同一列的元素上；
类似地，把以上的 行 改为 列 便得到 矩阵初等列变换 的定义。
矩阵初等行变换与初等列变换合称为 矩阵初等变换。
若矩阵 \(A\) 经过有限次的初等行变换变为矩阵 \(B\)，则矩阵 \(A\) 与矩阵 \(B\) 行等价；若矩阵 \(A\) 经过有限次的初等列变换变为矩阵 \(B\)，则矩阵 \(A\) 与矩阵 \(B\) 列等价；若矩阵 \(A\) 经过有限次的初等变换变为矩阵 \(B\)，则矩阵 \(A\) 与矩阵 \(B\) 等价。
当然列的用不着
首先有一个由系数构成的 \(n\times n\) 的矩阵
然后是一个由常数构成的 \(n\times 1\) 的列向量
把它们放在一起构成一个 \(n\times(n+1)\) 的增广矩阵：
我们遍历每一列，对于第 \(i\) 列选出最大的 未处理过的行 上的数作为主元，意味着加减消元后除了主元这一行外其他行的第 \(i\) 列都为 \(0\)。
选最大的作为主元的原因是避免精度误差。
然后就是加减消元了。
举个例子
增广矩阵就是
先是第 \(1\) 列，选 \(4\) 为主元。
\(4\) 在第 \(2\) 行，将第 \(2\) 行与第 \(1\) 行交换。
对于第 \(2\) 行，第一列上的 \(2\) 是 \(4\) 的 \(\dfrac{1}{2}\)。
化简得
对第 \(3\) 行的处理同理
现在到了第 \(2\) 列，未处理过的是 \(2,3\) 行，选最大的 \(\dfrac{3}{4}\) 为主元。
将第 \(3\) 行与第 \(2\) 行交换
消元得
第 \(3\) 列，未处理过的是第 \(3\) 行，选 \(4\) 作主元。
这样就把原来的矩阵通过初等变换，使得系数矩阵只有主对角线位置的元素非零，即一个对角矩阵。
上面那个矩阵的意思是
所以再把系数除过去就得到
请自行检验。
时间复杂度为 \(\mathcal{O}(n^3)\)。
当然方程还可能是无解或无穷多解。
考虑一个一元一次方程 \(ax=b\) 的解的情况：
- 无解：\(a=0,b\ne0\)；
- 无穷多解：\(a=b=0\)；
- 唯一解：\(a\ne0\)。
- 所以当发现某一列的主元为 \(0\) 时，因为主元是最大的，所以意味着这一列全都是 \(0\)，那么要么无解，要么有无穷多解。
\(\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;
}
关于判断无解和无穷多解
这下要具体到到底是无解还是无穷多解了。
这就是高斯-约旦的一个缺点：相比于普通高斯，它更难判断无解和无穷多解。
但是也是可以处理的。
具体地，我们用 \(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;
}
