高斯消元
高斯消元，一种十分实用的工具。
学高斯消元，首先需要知道什么是线性方程组：
\[\left\{\begin{matrix}
a_{11}x_1+a_{12}x_2+ \dots +a_{1m}x_m=b_1 \\
a_{21}x_1+a_{22}x_2+ \dots +a_{2m}x_m=b_2 \\
\vdots \\
a_{n1}x_1+a_{n2}x_2+ \dots +a_{nm}x_m=b_n \\
\end{matrix}\right.
\]
而高斯消元就可以帮我们解除这个方程所有的解。
为了方便，我们将他写成一个 \(n\) 行 \(m+1\) 列的 增广矩阵
\[\begin{bmatrix}
a_{11} & a_{12} & \dots & a_{1m} & b_1 \\
a_{21} & a_{22} & \dots & a_{2m} & b_2 \\
& & \vdots & & \\
a_{n1} & a_{n2} & \dots & a_{nm} & b_n
\end{bmatrix}
\]
接下来存在三中情况的初等行变换，使得方程的解不变：
- 交换两行
显然存在正确性
- 某行乘上非 \(0\) 的数
学过方程的都知道，这是等式的基本性质
- 把某行的 \(k\) 倍加到另一行
等式的基本性质中，等式左右同时加上相同的数，等式任然满足，而一个等式左右两边同时乘上 \(k\)，等式左右两边依然相等，也就是说让一个等式左右两边同时加上两个相同的数，等式仍成立。
因此我们就可以通过初等行变换，将一个增广矩阵变成一个梯形矩阵(以3*3为例)：
\[\begin{bmatrix}
a_{11} & a_{12} & a_{13} & b_1\\
0 & a_{22} & a_{23} & b_2\\
0 & 0 & a_{33} & b_3
\end{bmatrix}
\]
显然可以通过倒退求出方程的解。
这里讲一下无解情况
首先是矩阵的秩，将矩阵化简为梯形矩阵后，非零行的个数称之为矩阵的行秩(\(rank\))，同理还会有列秩，行秩与列秩，统称为矩阵的秩。
那么我们就是通过矩阵的秩来判断方程解的情况。
- 无解
系数矩阵的秩小于增广矩阵的秩
- 无数解
增广矩阵的秩 \(<n\)
- 唯一解
系数矩阵的秩 \(=\) 增广矩阵的秩 \(= 1\)。
高斯-约旦消元法
在高斯消元法的基础上进一步把梯形矩阵变成：
\[\begin{bmatrix}
1 & 0 & 0 & b_1\\
0 & 1 & 0 & b_2\\
0 & 0 & 1 & b_3
\end{bmatrix}
\]
代码
int guess(){
for(int i=1;i<=n;i++){
int maxh=i;
for(int j=1;j<=n;j++){
if(abs(a[j][i])>abs(a[maxh][i]) && j<i && abs(a[j][j])<EPS){
maxh=j;
}
if(abs(a[j][i])>abs(a[maxh][i]) && j>i){
maxh=j;
}
}
for(int j=1;j<=n+1;j++){
swap(a[i][j],a[maxh][j]);
}
if(abs(a[i][i])<EPS){
continue;
}
double d=a[i][i];
for(int j=1;j<=n+1;j++){
a[i][j]/=d;
}
for(int j=1;j<=n;j++){
if(i==j)continue;
double r=a[j][i];
for(int k=i;k<=n+1;k++){
a[j][k]-=r*a[i][k];
}
}
}
for(int i=1;i<=n;i++){
bool flag=1;
for(int j=1;j<=n;j++){
if(abs(a[i][j])>EPS){
flag=0;
}
}
if(flag && abs(a[i][n+1])>EPS){
return -1;//无解
}
}
for(int i=1;i<=n;i++){
if(a[i][i]<EPS){
return 0;//无限解
}
}
return 1;//有解
}
P5027 Barracuda
对于输入的每一次称重，都可以表示成一个方程。
之后我们可以枚举哪一个方程是错误的，然后根据题意判断是否会出现错误，最后就输出答案。
#include<bits/stdc++.h>
using namespace std;
const double EPS=0.000001;
int n;
double a[105][105];
void write(){
cout<<"\n";
for(int i=1;i<=n;i++){
for(int j=1;j<=n+1;j++)cout<<a[i][j]<<' ';
cout<<"\n";
}
return;
}
int guass(){
for(int i=1;i<=n;i++){
int maxh=i;
for(int j=1;j<=n;j++){
if(abs(a[j][i])>abs(a[maxh][i]) && j<i && abs(a[j][j])<EPS){
maxh=j;
}
if(abs(a[j][i])>abs(a[maxh][i]) && j>i){
maxh=j;
}
}
for(int j=1;j<=n+1;j++){
swap(a[i][j],a[maxh][j]);
}
if(abs(a[i][i])<EPS){
continue;
}
double d=a[i][i];
for(int j=1;j<=n+1;j++){
a[i][j]/=d;
}
for(int j=1;j<=n;j++){
if(i==j)continue;
double r=a[j][i];
for(int k=i;k<=n+1;k++){
a[j][k]-=r*a[i][k];
}
}
}
for(int i=1;i<=n;i++){
bool flag=1;
for(int j=1;j<=n;j++){
if(abs(a[i][j])>EPS){
flag=0;
}
}
if(flag && abs(a[i][n+1])>EPS){
return -1;
}
}
for(int i=1;i<=n;i++){
if(abs(a[i][i])<EPS || ceil(a[i][n+1])!=floor(a[i][n+1]) || a[i][n+1]<=EPS){
return 0;
}
}
return 1;
}
int m[105],w[105];
int f[105][105];
bool flag;
signed main(){
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>n;
for(int i=1;i<=n+1;i++){
cin>>m[i];
for(int j=1;j<=m[i];j++){
cin>>f[i][j];
}
cin>>w[i];
}
int id=0;
int ans=0;
for(int i=1;i<=n+1;i++){
memset(a,0,sizeof(a));
int cnt=0;
for(int j=1;j<=n+1;j++){
if(i==j){
cnt++;continue;
}
for(int k=1;k<=m[j];k++){
a[j-cnt][f[j][k]]=1;
}
a[j-cnt][n+1]=w[j];
}
int res=guass();
if(res!=1){
continue;
}
int maxn=0,maxid=0;
bool f=0;
for(int j=1;j<=n;j++){
if(a[j][n+1]>maxn){
maxn=a[j][n+1];
maxid=j;
}
}
for(int j=1;j<=n;j++){
if(maxid!=j && a[j][n+1]==maxn){
f=1;
break;
}
}
if(f){
continue;
}
if(flag){
cout<<"illegal";
return 0;
}
flag=1;
ans=maxn;
id=maxid;
}
if(!id){
cout<<"illegal";return 0;
}
cout<<id;
return 0;
}
P4035 [JSOI2008] 球形空间产生器
题目给出了 \(n\) 个未知数以及 \(n+1\) 个方程，则不适合使用高斯消元，也用不了高斯消元。我们考虑让两个式子相减。则发现他们从原本的：
\[(a_1-x_1)^2+(a_2-x_2)^2+\dots+(a_1-x_n)^2=r^2 \\
(b_1-x_1)^2+(b_2-x_2)^2+\dots+(b_n-x_n)^2=r^2
\]
变成了：
\[2\cdot(a_1-b_1)+2\cdot(a_2-b_2)+\dots+2\cdot(a_n-b_n)=a_1^2-b_1^2+a_2^2+b_2^2+\codt+a_n^2-b_n^2
\]
把所有相邻的方程相减，就得到了 \(n\) 个方程以及 \(n\) 个未知数，可以使用高斯消元求解。
#include<bits/stdc++.h>
using namespace std;
const double EPS=0.000001;
int n;
double a[105][105];
double b[105][105];
int guass(){
for(int i=1;i<=n;i++){
int maxh=i;
for(int j=1;j<=n;j++){
if(abs(a[j][i])>abs(a[maxh][i]) && j<i && abs(a[j][j])<EPS){
maxh=j;
}
if(abs(a[j][i])>abs(a[maxh][i]) && j>i){
maxh=j;
}
}
for(int j=1;j<=n+1;j++){
swap(a[i][j],a[maxh][j]);
}
if(abs(a[i][i])<EPS){
continue;
}
double d=a[i][i];
for(int j=1;j<=n+1;j++){
a[i][j]/=d;
}
for(int j=1;j<=n;j++){
if(i==j)continue;
double r=a[j][i];
for(int k=i;k<=n+1;k++){
a[j][k]-=r*a[i][k];
}
}
}
for(int i=1;i<=n;i++){
bool flag=1;
for(int j=1;j<=n;j++){
if(abs(a[i][j])>EPS){
flag=0;
}
}
if(flag && abs(a[i][n+1])>EPS){
return -1;
}
}
for(int i=1;i<=n;i++){
if(abs(a[i][i])<EPS || ceil(a[i][n+1])!=floor(a[i][n+1]) || a[i][n+1]<=EPS){
return 0;
}
}
return 1;
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>n;
for(int i=1;i<=n+1;i++){
for(int j=1;j<=n;j++){
cin>>b[i][j];
}
}
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
a[i][j]=2*(b[i][j]-b[i+1][j]);
a[i][n+1]+=b[i][j]*b[i][j]-b[i+1][j]*b[i+1][j];
}
}
guass();
for(int i=1;i<=n;i++){
if(abs(a[i][n+1])<EPS){
cout<<"0.000 ";
}else{
cout<<fixed<<setprecision(3)<<a[i][n+1]<<' ';
}
}
return 0;
}
浙公网安备 33010602011771号