高斯消元法
可以用于求解线性方程组,即n元1次方程组。利用矩阵,大致思路与普通解方程方法类似。只是更具一般性。将系数与右侧的常数存成一个矩阵,然后每次用第i行消去下面每行的第i个系数,最后就会得到一个一元方程,然后从后到前依次代回即可。
然后就是精度的问题,因为计算机中没有分数,所以只能用double来存,但是double的精度也是有限的,所以要想办法维护精度。
第一种优化办法就是在用第i行去消第i个系数时,从下面所有行中选择第i个系数的绝对值最大的一个与当前行交换,这样可以使每次消元时差距尽量大,然后就可以保持精度了。
第二种优化就是在进行消元时,不是用一个变量来储存两行首元素的商,而是每次直接计算。这样就必须防止第一行提前被修改,所以倒着进行修改。
代码:
1 #include<cstdio> 2 #include<cstring> 3 #include<iostream> 4 #include<cmath> 5 using namespace std; 6 const int N=110,M=1010; 7 int n,m; 8 double a[N][N],tmp[N][N],ans[N]; 9 inline void work() 10 { 11 for(int i=1;i<=n;++i) 12 { 13 int kk=i; 14 for(int j=i+1;j<=n;++j) 15 if(abs(a[kk][i])<abs(a[j][i])) 16 kk=j; 17 if(kk!=i) 18 for(int j=1;j<=n+1;++j) 19 swap(a[kk][j],a[i][j]); 20 for(int j=i+1;j<=n;++j) 21 { 22 for(int k=n+1;k>=i;--k) 23 a[j][k]=a[j][k]-a[j][i]/a[i][i]*a[i][k]; 24 } 25 } 26 27 for(int i=n;i>=1;--i) 28 { 29 for(int j=n;j>i;--j) 30 { 31 a[i][j]*=ans[j]; 32 a[i][n+1]-=a[i][j]; 33 } 34 ans[i]=a[i][n+1]/a[i][i]; 35 } 36 } 37 int main() 38 { 39 scanf("%d%d",&n,&m); 40 for(int i=1;i<=n;++i) 41 for(int j=1;j<=n;++j) 42 scanf("%lf",&tmp[i][j]); 43 for(int k=1;k<=m;++k) 44 { 45 memset(ans,0,sizeof(ans)); 46 memset(a,0,sizeof(a)); 47 for(int i=1;i<=n;++i) 48 for(int j=1;j<=n;++j) 49 a[i][j]=tmp[i][j]; 50 for(int j=1;j<=n;++j) 51 scanf("%lf",&a[j][n+1]); 52 work(); 53 for(int i=1;i<=n;++i) 54 printf("%.9llf ",ans[i]); 55 printf("\n"); 56 } 57 return 0; 58 }