这里暂时并不会涉及到复杂的线性代数中的情况
化为上三角行列式之,有三种情况
- n个变元,n个方程:唯一解
- n个变元,小于n个方程,且其余方程不存在矛盾:无穷解
- 存在矛盾:无解
算法流程
1.从第一列开始,找到每一列绝对值最大的那一行
2.把上一步找到的那一行放到最上面
3.把最上面一行该列的系数变为1
4.把该列其他行系数变为0
5.如果有唯一解,从最后一个式子开始向上逐渐消元直至得到最简上三角行列式
代码实现
#include <iostream>
#include <cmath>
#include <algorithm>
#include <stdio.h>
using namespace std;
const int N = 110;
const double eps = 1e-6; // 这个数字的定法是个比较玄学的问题,因为不好说多小才算是0,0.01就不应该算作0,所以1e-6应该说是经验值
int n;
double a[N][N];
int guass()
{
// 完成高斯消元的模拟过程并返回解的情况
int c, r;
for (c = 0, r = 0; c < n; ++ c) // 列,每次循环确定一个未知量的解
{
// 第一步:找到c列绝对值最大的行数
int t = r; // c列系数最大的一行
for (int i = r; i < n; ++ i) // 行
if (fabs(a[i][c]) > fabs(a[t][c]))
t = i;
/**
* 实际目的为判断是否为0,但是由于精度问题,可能出现0.00000001这种类似的情况,所以采用这种方式,在浮点数二分时也采用过这种方式
* 这里的contineu不太好理解
* 我们需要从r的意义入手,我们希望每次r都能确定一个解,如果当前列的系数均为0,我们是没办法确定一个解的,所以去看下一列
* 所以最终如果是唯一解,那么r应该是等于方程个数的,因为每个方程确定一个解,如果不等于,就要看方程右边是不是存在某个常数不为0,如果存在,说明存在0 = 非0这种式子,所以是无解,否则是无穷多解
*/
if (fabs(a[t][c]) < eps) continue;
// 第二步:将t行换到最上方还未确定的一行
for (int i = 0; i < n + 1; ++ i) swap(a[r][i], a[t][i]);
// 第三步:将首位系数变为1
for (int i = n; i >= 0; -- i) a[r][i] /= a[r][c];
// 第四步:将下方c列系数变为0
for (int i = r + 1; i < n; ++ i)
for (int j = n; j >= 0; -- j) // 需要从后向前更新,以保留该行第一个系数
a[i][j] -= a[i][c] * a[r][j];
++ r;
}
/**
* 如果是唯一解,那么每一行确定一个一个解
* 假设有3个方程,3个未知量,那么r=0时处理第一个未知量,=1时第二个,=2时第三个,=3时结束了,也就是r最终停止的位置是不确定解的位置
* 所以如果是无解或者无穷解,r-1是最后一个确定方程解得位置,r是第一个全零行,所以判断是无解还是无穷解是从r开始的,而非r+1
*/
// 第五步:判断是否有解并处理可消去可消系数
if (r < n)
{
for (int i = r; i < n; ++ i)
if (fabs(a[r][n]) > eps)
return 2;
return 1;
}
// 用i行下面每一行来把i行对应的系数消为0,常数列对应改变
// 这里i必须从下向上走,因为消i行的某个系数时,必须保证j行其它系数为0了才能保证得到正确的结果,所以要先消下面的系数
for (int i = n - 1; i >= 0; -- i)
for (int j = i + 1; j < n; ++ j)
a[i][n] -= a[j][n] * a[i][j];
return 0;
}
int main()
{
cin >> n;
for (int i = 0; i < n; ++ i)
for (int j = 0; j < n + 1; ++ j)
cin >> a[i][j];
int t = guass();
if (!t)
{
// 此时有唯一解,最终方程的目标格式为每个式子只有一个x且该x前系数为1,所以最后的常数就是解
for (int i = 0; i < n; ++ i)
printf("%.2lf
", a[i][n]);
}
else if (t == 1) cout << "Infinite group solutions" << endl;
else cout << "No solution" << endl;
return 0;
}