高斯-约旦消元

高斯-约旦消元

介绍

之前看高斯消元感觉好麻烦，不知道怎么用程序实现，在网上搜了搜发现有个叫“高斯约旦消元的好东西”，实现起来非常方便不易错

高斯约旦消元直接把原矩阵转换为简化阶梯矩阵，大致流程如下：

从第一列开始依次枚举每一列未知数，每一次选择这一列上的系数不为 \(0\) 的一行（枚举到第 \(i\) 列则从第 \(i+1\) 行开始往下找，因为上面的行中有已经消元完毕的未知数，如果所有找到的行均为 \(0\) 则说明无解），利用初等行变换把这一行交换到与列编号相等的行号（第 \(i\) 列第 \(i\) 行），然后使用倍加行变换把矩阵中这一列上其他所有位置的系数都变为 \(0\)

最后矩阵中对角线上则是仅剩的未知数，用右侧的常数除以左侧的系数即可得到答案

Code

下面给出代码，中间有一些细节可以参考

#include<bits/stdc++.h>
using namespace std;
#define N 105
const double eps=1e-6;
int n;
double a[N][N],tmp;
signed main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;++i)
		for(int j=1;j<=n+1;++j)
			scanf("%lf",&a[i][j]);
	for(int i=1,maxn;i<=n;++i)//选择主元（列） 
	{
		maxn=i;
		for(int j=i+1;j<=n;++j)//从当前行开始向下枚举找绝对值最大值（第1~i行上有第1~i列的已经消元的部分，不能换下来） 
			if(fabs(a[maxn][i])<fabs(a[j][i]))//这里找绝对值最大值和找不为0的值没什么区别
				maxn=j;
		if(maxn!=i) swap(a[maxn],a[i]);//i行i列不是最大值就把最大绝对值那一行换到i行
		if(fabs(a[i][i])<eps) return puts("No Solution"),0;//最大绝对值都为0就说明无解 
		for(int j=1;j<=n;++j)//初等行变换消元 
			if(i!=j)
			{
				tmp=a[j][i]/a[i][i];
				a[j][i]=0;
				for(int k=i+1;k<=n+1;++k)//第i列左边都为0了，不用扫描
					a[j][k]-=a[i][k]*tmp;
			}
	}
	for(int i=1;i<=n;++i)//此时每一行都只有一个未知数系数不为0
		printf("%.2lf\n",a[i][n+1]/a[i][i]);//记得要除以系数 
    return 0;
}

---

该文为本人原创，转载请注明出处

博客园传送门