高斯消元详解
题目传送门：【模板】高斯消元法，[SDOI2006]线性方程组。
引入
给定一个 \(n\) 元一次方程组：
其中，\(a_{i,j},b_i\) 均为常数，求解 \(x_1,x_2,x_3,\cdots,x_n\)。
那么，我们便可以利用高斯消元法 \(\mathcal O\left(n^3\right)\) 求解。
高斯消元法
基本思想
首先将方程的增广矩阵利用行初等变换化为行最简形，然后以线性无关为准则对自由未知量赋值，最后列出表达方程组通解。
上面这句话看不懂忽略，下文有解释。
过程
其实就可以理解为加减消元法求解后回代。
例如对于方程组：
将 \(x\) 项系数化为 \(1\)：
上下相减得：
将 \(y\) 项系数化为 \(1\)：
上下相减：
将 \(z\) 项系数化为 \(1\)：
以上为加减消元，以下为回代。
将 \(z=3\) 带入 \(y+8z=26\) 得：
解得：
将 \(\begin{cases}y=2\\z=3\end{cases}\) 带入 \(x+2y+4z\) 得：
解得：
思想解释
首先将方程的增广矩阵利用行初等变换化为行最简形，然后以线性无关为准则对自由未知量赋值，最后列出表达方程组通解。
- “方程的增广矩阵”即系数矩阵后加了一列作为结果。
- “利用行初等变换化为行最简形”即如过程中的依次化当前项系数为 \(1\)。
- “以线性无关为准则对自由未知量赋值”即加减消元法。
- “列出表达方程组通解”即回代求解。
实现
首先，为了便于程序存储和分析，我们使用矩阵的形式来代替方程组。
同样是上面那个例子，得出矩阵：
高斯消元的加减消元的前提是主项系数都是 \(1\)。
因此我们得到：
相减：
化系数为 \(1\)：
再次相减：
系数化为 \(1\)：
至此，你可以发现此时的矩阵是一个“上三角“，即主对角线上系数均为 \(1\)，主对角线下均为 \(0\)。
这种时候，考虑如何求解 \(x_i\)。
很明显，当我们求得了 \(x_{i+1}\sim x_n\) 后，我们仅仅需要使 \(\large b_i\leftarrow b_i-\sum\limits_{j=i+1}^na_{i,j}\times x_j\)，然后 \(b_i\) 即答案。
因为 \(a_{i,1}\sim a_{i,i-1}\) 均为 \(0\)，自然不需要考虑。
就此，高斯消元实现。
无解、定解与无数解
定解
这应该是最好判断的，不是其余两种就行。
这不废话吗
无解
明显地，\(0=b_i\land b_i\ne0\) 是不成立的，即无解。
放到高斯消元中，就是若某行前 \(n\) 个数（所有系数）均为 \(0\)，但最终结果 \(b_i\) 却不为 \(0\)，则无解。
无数解
一般来讲，”无数解“都是至最终解出来了一个恒等式。那么在一个方程中，若等式左边的未知项等于一个常数，则未知项的值是为 \(0\) 的。那么未知项的系数也就是 \(0\)。即：若某行前 \(n\) 个数和最终结果 \(b_i\) 均为 \(0\)，则有无数解。
存储
一般来讲，为了方便运算，我们都不会单独放一个 \(b\) 数组用于存储答案，而是直接存储增广矩阵。
高斯消元明显是一个 \(n\times n\) 的系数矩阵与 \(b\) 数组形成的一个 \(n\times(n+1)\) 的增广矩阵。
开一个数组 \(a[N+1][N+1+1]\) 即可。
消元
一般来讲，”化系数为 \(1\)“时，我们都会找当前项系数绝对值最大的将其化为 \(1\)。
为什么？看到无解、无数解的判定。为了更快得出是否有解。
[SDOI2006]线性方程组 AC 代码
具体参见注释。
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
const double dx=1e-6;//dx:因为double精度问题，可近似视为0
const int N=100;
int n;
double a[N+1][N+1],x[N+1];
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
scanf("%d",&n);
for(int i=1;i<=n;i++){
for(int j=1;j<=n+1;j++)scanf("%lf",&a[i][j]);
}
int r,c;
//r:当前方程
//c:当前项
for(r=c=1;r<=n&&c<=n;r++,c++){
//找系数绝对值最大的方程
int u=r;
for(int i=r+1;i<=n;i++){
if(abs(a[i][c])>abs(a[u][c]))u=i;
}if(u!=r)swap(a[r],a[u]);//交换(便于处理)
//这一步可以使循环结束后运行r++,c++时,r并没有增加，即直接跳过这一项，先计算后面的项
if(abs(a[r][c])<dx){
r--;
continue;
}//这一步其实合并了化系数为1和加减消元,实在不能理解手推一下看看
for(int i=r+1;i<=n;i++){
if(abs(a[i][c])>dx){
double pl=a[i][c]/a[r][c];
for(int j=c;j<=n+1;j++)a[i][j]-=a[r][j]*pl;
a[i][c]=0;
}
}
}//若剩余方程存在且剩余方程内有非0系数,则无解
for(int i=r;i<=n;i++){
if(abs(a[i][c])>dx){
printf("-1\n");
return 0;
}
}//无数解
if(r<=n){
printf("0\n");
return 0;
}//回代还原
for(int i=n;i>=1;i--){
for(int j=i+1;j<=n;j++)a[i][n+1]-=a[i][j]*x[j];
x[i]=a[i][n+1]/a[i][i];
}//输出
for(int i=1;i<=n;i++)printf("x%d=%.2lf\n",i,x[i]);
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
【模板】高斯消元法 AC 代码
和上面几乎一摸一样......
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
const double dx=1e-6;
const int N=100;
int n;
double a[N+1][N+1],x[N+1];
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
scanf("%d",&n);
for(int i=1;i<=n;i++){
for(int j=1;j<=n+1;j++)scanf("%lf",&a[i][j]);
}
int r,c;
for(r=c=1;r<=n&&c<=n;r++,c++){
int u=r;
for(int i=r+1;i<=n;i++){
if(abs(a[i][c])>abs(a[u][c]))u=i;
}
if(u!=r)swap(a[r],a[u]);
if(abs(a[r][c])<dx){
r--;
continue;
}
for(int i=r+1;i<=n;i++){
if(abs(a[i][c])>dx){
double pl=a[i][c]/a[r][c];
for(int j=c;j<=n+1;j++)a[i][j]-=a[r][j]*pl;
a[i][c]=0;
}
}
}
for(int i=r;i<=n;i++){
if(abs(a[i][c])>dx){
printf("No Solution\n");//改动1
return 0;
}
}
if(r<=n){
printf("No Solution\n");//改动2
return 0;
}
for(int i=n;i>=1;i--){
for(int j=i+1;j<=n;j++)a[i][n+1]-=a[i][j]*x[j];
x[i]=a[i][n+1]/a[i][i];
}
for(int i=1;i<=n;i++)printf("%.2lf\n",x[i]);//改动3
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
浙公网安备 33010602011771号