高斯消元学习体会

高斯消元就像初中学的多元一次方程组（且当$n$元时有$n$个方程）的一种解法。 ## 一、解法 其中经过多次($n^3$)消元后可呈倒三角形式 $\begin{bmatrix} 1 & -1 & 3\\ 0 & 8 & 5\\ 0 & 0 & 4 \end{bmatrix}$ $\begin{bmatrix} 2\\ 3\\ 8 \end{bmatrix}$ 这样将$4x_3=8\Rightarrow x_3=2$带入②和①得 $\begin{bmatrix} 1 & -1 & 0\\ 0 & 8 & 0\\ 0 & 0 & 0 \end{bmatrix}$ $\begin{bmatrix} -4\\ -7\\ 0\end{bmatrix}$ 可知$8x_2=-7\Rightarrow x_2=-\frac{8}{7}$ 再次迭代可得 $\begin{bmatrix} 1 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{bmatrix}$ $\begin{bmatrix} -\frac{36}{7}\\ 0\\ 0\end{bmatrix}$ $x_1=-\frac{36}{7}$就解出来了。 为了解决问题方便，我们可以将矩阵坐标为$i,j$的位置化为1，即 $\begin{bmatrix} 1 & -1 & 3\\ 0 & 8 & 5\\ 0 & 0 & 4 \end{bmatrix}$ $\begin{bmatrix} 2\\ 3\\ 8 \end{bmatrix}\rightarrow\begin{bmatrix} 1 & -1 & 3\\ 0 & 1 & \frac{5}{8}\\ 0 & 0 & 1 \end{bmatrix}$ $\begin{bmatrix} 2\\ \frac{3}{8}\\ 2 \end{bmatrix}$ 就可以看得更清楚，让程序直接从$i+1$做计算 ## 二、精度 为了保证精度，我们可以让小的去扩大，而不是去对大的除，同时要保证被除数不能是0，防止RE（这里是因为当小数点后8位均为0时已无意义，甚至会爆double） ## 三、判断 判断是否有唯一解 如果做到某一列发现这一列的系数全部为0，说明所在这一列的$x$为“自由元”，则不为唯一解，这个$x$可以取任意值。 如果做完后发现系数有矛盾即消不掉的系数，则该矩阵无解。 ## 四、其他 如上这种做法，因为计算机不能存储分数（可以手动实现）精度丢失过大，但在一般的题中已经够用。（见代码高斯消元v1.0） 为防止~~变态题~~精度误差太大，做除法之前可以求出这一列的最小公倍数，互相消元时就不会出现小数了（还是要避免0的情况出现），但是这样就不能保证$i,j$位置上的数是1了，迭代回溯时要注意这一点。（代码待补充） ### 高斯消元v1.0 ```cpp #include<cstdio> #include<algorithm> double abs(double x) { if(x<0) return -x; return x; } double a[111][111]; int n; bool gs() { for(int i=1;i<=n;i++) { double maxx=a[i][i]; int tmp=i; for(int j=i;j<=n;j++) if(abs(a[j][i])>abs(maxx)) { maxx=a[j][i]; tmp=j; } if(abs(maxx)<=0.0000001) return false; if(i<n) std::swap(a[i],a[tmp]); double upd=a[i][i]; for(int j=i;j<=n+1;j++) a[i][j]/=upd; for(int k=i+1;k<=n;k++) { upd=a[k][i]; for(int j=i;j<=n+1;j++) a[k][j]-=upd*a[i][j]; } } for(int i=n;i>=1;i--) { double upd=a[i][n+1]; for(int j=i-1;j>=1;j--) { a[j][n+1]-=a[j][i]*upd; a[j][i]=0; } } return true; } int main() { scanf("%d",&n); for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) scanf("%lf",&a[i][j]); if(gs()) { for(int i=1;i<=n;i++) printf("%.2lf\n",a[i][n+1]); } else puts("No Solution"); return 0; } ```