高斯消元&&luogu3389
高斯消元(Gauss)
高斯消元和我们做二元一次方程组差不多
流程:
1.把系数和右边的值就是用二维数组存下来->转化成矩阵
我们的目标是把这个矩阵装换成 上三角的形式
对角线系数全部为1,1下面都为0,为了下面的回带
1 4 2 3 0 1 7 9 0 0 1 2 0 0 0 1
2.利用 加减消元和等式两边除以一个数,一列一列的进行消元
顺便判断一下是否有解,对角线上系数不为0
3.求出上三角之后,我们倒着回代一下就可以求取解了
当选取主元的时候,由于是double类型,当对角线的系数太小时,此时用它做除数会带来误差扩散,使结果严重失真。所以我们在消元的过程中,如果出现主元相差较大,要选取最大数作为主元,并交换行列,(当然,消元完毕的上边不能考虑在内)
---参考数学一本通
代码
1 #include <iostream> 2 #include <cmath> 3 #include <cstdio> 4 using namespace std; 5 6 const double eps=1e-6; 7 int n; 8 double a[110][110]; 9 double ans[110]; 10 11 int main() 12 { 13 scanf("%d",&n); 14 for(int i=1; i<=n; ++i) 15 for(int j=1; j<=n+1; ++j) 16 scanf("%lf", &a[i][j]); 17 18 for(int i=1; i<=n; ++i) { 19 int pivot=i; 20 for(int j=i+1; j<=n; ++j)//选取较大主元 21 if(fabs(a[j][i]) > fabs(a[pivot][i])) pivot=j; 22 if(abs(a[pivot][i]) < eps) { //判断有无解,无穷解也当做无解 23 printf("No Solution"); 24 return 0; 25 } 26 if(pivot!=i) swap(a[i],a[pivot]);//直接交换 27 double tmp=a[i][i]; 28 for(int j=i; j<=n+1; ++j) { 29 a[i][j]/=tmp;//系数化为1 30 } 31 for(int j=i+1;j<=n;j++) {//下面的化为0 32 tmp=a[j][i]; 33 for(int k=i;k<=n+1;k++) { 34 a[j][k]-=a[i][k]*tmp; 35 } 36 } 37 } 38 ans[n]=a[n][n+1]; 39 for(int i=n-1; i>=1; i--) { 40 ans[i]=a[i][n+1]; 41 for(int j=i+1; j<=n; ++j) 42 ans[i]-=a[i][j]*ans[j]; 43 }//回带 44 for(int i=1;i<=n;++i) 45 printf("%.2lf\n",ans[i]); 46 }