概述
高斯消元是一种重要的求解线性方程组的方法,高斯消元其实没有名字那么高大上,本质上就是减消元(这不是初中知识么)。大致过程就是用矩阵存储各项的系数,把第i个方程的前i项通过加减前面的式子变为0,然后求解即可。
过程
存储
首先用矩阵存储各项系数,我们定义系数矩阵,矩阵的第i行第j列存储第i个方程的第j项的系数。增广矩阵在系数矩阵的基础上加
上一列,记录每个等式等号右边的常数。
for(int i=1;i<=n;i++){ for(int j=1;j<=n+1;j++){ scanf("%lf",&a[i][j]); } }
消元
消元的过程就是通过方程的加减把一些系数消去,最终的结果应该是第一个方程剩余n个未知数,第二个方程剩n-1个...最后一个方程只剩下一个未知数,我们设第i个方程剩下的系数是x[i],x[i+1]...x[n],这样,消元后的结果用矩阵表示出来就是一个阶梯形的矩阵,第i行的元素有a[i][i]~a[i][n]。
那么怎么进行操作呢?显然,我们可以枚举1~n,再枚举要消去第i项系数的行j(i+1~n),用第i行乘一个数和第j行进行加减使a[j][i]为0,那么第i行需要乘的数显然就是ratio=a[j][i]/a[i][i],然后将第j行每一个位置k减去a[i][k]*ratio即可。
for(int j=i+1;j<=n;j++){ double ratio=a[j][i]/a[i][i]; for(int k=1;k<=n+1;k++) a[j][k]=a[j][k]-ratio*a[i][k]; }
但是有一个问题,如果a[i][i]=0咋办?除法会炸掉。而且如果是一个极小的数,也可能会引起精度误差。
办法就是,找出第i项系数绝对值最大的行p,把这一行和第i行交换。需要注意的是,如果交换后a[i][j]仍然为0,说明有无数组解(xi可以取任何值)。
const double eps=1e-8;
int sign(double x) { if(fabs(x)<eps) return 0; if(x>0) return 1; else return -1; }//判断浮点数是否为零
int p=i; for(int j=i;j<=n;j++){ if(fabs(a[j][i])>fabs(a[p][i])) p=j; } for(int j=1;j<=n+1;j++) swap(a[i][j],a[p][j]); if(sign(a[i][i]==0)){ cout<<"No Solution"; flag=false; break; }
// 蜜汁乱码,可以直接看下面代码
到这里消元的过程就完成了。
求解
求解的思路就十分明了了,现在我们最后一项只有一个未知数,所以我们从后往前枚举,每次求出一个未知数后,把该值分别代入前面的方程,再解前面一个方程,就能求出所有未知数。
模板题:洛谷3389
代码:
#include<cstdio> #include<iostream> #include<algorithm> #include<cstring> #include<cmath> using namespace std; const double eps=1e-8; double a[110][110]; double x[110]; int sign(double x) { if(fabs(x)<eps) return 0; if(x>0) return 1; else return -1; }//判断浮点数是否为零 int main() { int n; bool flag=true; cin>>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;i<=n;i++){ int p=i; for(int j=i;j<=n;j++){ if(fabs(a[j][i])>fabs(a[p][i])) p=j; } for(int j=1;j<=n+1;j++) swap(a[i][j],a[p][j]); if(sign(a[i][i]==0)){ cout<<"No Solution"; flag=false; break; } for(int j=i+1;j<=n;j++){ double ratio=a[j][i]/a[i][i]; for(int k=1;k<=n+1;k++) a[j][k]=a[j][k]-ratio*a[i][k]; } } for(int i=n;i>=1;i--){ x[i]=a[i][n+1]/a[i][i]; for(int j=1;j<=i-1;j++) a[j][n+1]-=a[j][i]*x[i]; } if(flag) for(int i=1;i<=n;i++) printf("%.2f\n",x[i]); return 0; }