高斯消元详解
简述
x+y+3z=6 - ①
2x+4y+3z=8 - ②
2x+3y+12z=4 - ③
想必所有人都解过上面的方程,那我们在这里模仿一下解题步骤
(1)先利用①式消去②式和③式的x项
=> ②-2①:2y-6z=-4 - ④
③-2①:y+6z=-8 - ⑤
(2)利用④式消去⑤式的y项
=> ⑤-1/2④:3z=6
∴ z=2
(3)将z=2带入④式
=> 2y-12=-4
∴ y=4
(4)将y=4,z=2带入①式
=> x+4+6=6
∴ x=-4
(5)得出解:x=-4,y=4,z=2
我们为了算法容易实现,将所有式子的系数和它们右面的数写在一个矩阵A之中
然后我们考虑之前是如何解出答案的,回顾前面的例子不难发现我们无非是一个一个将未知数消去,这样我们就可以得到了解题思路
即:每一次用第i行的第i个数aii作为标准,对大于i的每一行k的大于i的第j个数均减去aki/aii*aki,这样我们就可以成功的将之后的每一个式子消去一个未知数了。
但这实际是远远不够的,因为有数据可以卡掉这种做法,比如:
我们不难发现如果我们的矩阵长这个样子,那我们就无法将第三个式子的第二项消掉,所以我们应该想一种可以不受某项为0这种情况干扰的做法。要做到这一点我们要知道对于一个方程组,交换任意方程的位置解不变,所以我们对于第i行只要找到大于等于i的行中第一个第i个位置不为0的行与之交换就行了。
具体实现请参照下面的代码。
代码(p3389)
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#include<cmath>
#include<cstdlib>
#include<queue>
#include<ctime>
#include<vector>
#include<set>
#include<map>
#include<stack>
using namespace std;
#define ri register int
const double eps=1e-8;
double a[1100][1100];
int n;
inline int read(){
int x=0,f=1;char s=getchar();
while(s<'0'||s>'9'){if(s=='-')f=-1;s=getchar();}
while(s>='0'&&s<='9'){x=(x<<3)+(x<<1)+(s-'0');s=getchar();}
return x*f;
}
int main()
{ n=read();
for(ri i=0;i<n;i++){
for(ri j=0;j<n;j++)
a[i][j]=read();
a[i][n]=read();
}
for(ri i=0;i<n;i++){
int p=i;
for(ri j=i;j<n;j++)
if(fabs(a[j][i]-a[p][i])<=eps)
p=j;
for(ri j=0;j<=n;j++)swap(a[i][j],a[p][j]);
if(fabs(a[i][i])<=eps){
printf("No Solution\n");
return 0;
}
for(ri j=i+1;j<=n;j++)a[i][j]/=a[i][i];
for(ri j=0;j<n;j++)
if(j!=i)
for(ri k=i+1;k<=n;k++)
a[j][k]-=a[j][i]*a[i][k];
}
for(ri i=0;i<n;i++)
printf("%0.2lf\n",a[i][n]);
return 0;
}