高斯消元
高斯消元
高斯消元可真是复杂啊。。。。。
在 OI 中一般用于两点:求解线性方程组(不常见) & 求线性基(常见)
高斯消元求解n元一次线性方程组的板子题:
P3389 【模板】高斯消元法
举个栗子:
2x + y - z = 8
-3x - y + 2z = -11
-2x + y + 2z = -3
先将它存到矩阵中:
②+①* (2/3)
③+①
接着对①变换
得到x,y,z;
但是我们想到,如果它有在原方程中就有两个或多个方程本质上是一样的,那么这就是无解
比如:
最后得出:
0=9,矛盾 这就是无解
又比如:
元的数量 > 方程式的数量 这就是无解
这里我们引入一个定理:
一个矩阵的行列式如果不为0,方程组有唯一解,否则无解或者无穷多解
貌似就可以用行列式求解了
行列式:
上代码(from lh)(注释from lz):
#include<iostream> #include<algorithm> #include<cstdio> #include<cstring> #include<cmath> #include<cstdlib> using namespace std; typedef long long ll; typedef long double ld; typedef pair<int,int> pr; const double pi=acos(-1); #define rep(i,a,n) for(int i=a;i<=n;i++) #define per(i,n,a) for(int i=n;i>=a;i--) #define Rep(i,u) for(int i=head[u];i;i=Next[i]) #define clr(a) memset(a,0,sizeof a) #define pb push_back #define mp make_pair #define fi first #define sc second ld eps=1e-9; ll pp=1000000007; ll mo(ll a,ll pp){if(a>=0 && a<pp)return a;a%=pp;if(a<0)a+=pp;return a;} ll powmod(ll a,ll b,ll pp){ll ans=1;for(;b;b>>=1,a=mo(a*a,pp))if(b&1)ans=mo(ans*a,pp);return ans;} ll read(){ ll ans=0; char last=' ',ch=getchar(); while(ch<'0' || ch>'9')last=ch,ch=getchar(); while(ch>='0' && ch<='9')ans=ans*10+ch-'0',ch=getchar(); if(last=='-')ans=-ans; return ans; } //head(lh的代码头) int n,m; double a[100][100]; bool check(int k){ if(fabs(a[k][n+1])<eps)return 1;// 如果小于0,返回1 for(int i=1;i<=n;i++) if(fabs(a[k][i])>eps)return 1;//第k行的每个数都大于0 return 0; } int main(){ n=read();m=read();//m*n的矩阵, // a_i,1 a_i,2 ... a_i,n a_i,n+1 for(int i=1;i<=m;i++) for(int j=1;j<=n+1;j++)a[i][j]=read();//输入矩阵元素(因为除了这个m*n的矩阵,还要有一行答案qwq) (看洛谷题吧qwq) for(int j=1;j<=m;j++){ for(int k=1;k<=n+1;k++)cout<<a[j][k]<<" "; puts(""); }//先把矩阵输出辽一遍 int flag=0; for(int i=1;i<=n;i++){ int t=i; while(a[t][i]==0 && t<=n)t+=1; if(t==n+1){ flag=1; continue; }//if判断是否有一行全为0,如果全为0 那么我们就少了至少1个方程,那么这个方程显然无唯一解 for(int j=1;j<=n+1;j++)swap(a[i][j],a[t][j]);//交换第i行和第t行的值,目的在于交换第i行和第t行后可以使a[1][1]!=0 double kk=a[i][i];//找到对角线qwq for(int j=1;j<=n+1;j++)a[i][j]/=kk;//把这一行的每一项都除以对角线,使得对角线为1 for(int j=1;j<=m;j++) //这里真的要详细地讲一讲咯 if(i!=j){//首先i!=j这样就不再消对角线上的数了 double kk=a[j][i];//定义kk为第j行第i列的数 for(int k=1;k<=n+1;k++)//把每一项消成0qwq(先把第一列除了对角线全消为0然后第二第三……) // 关于处理对角线这里以i=1,j=3做例子:当k=3时,a[3][3]-=a[3][1]*a[1][3]而这个时候a[3][1]已经为0,故对对角线无影响 a[j][k]-=kk*a[i][k]; } puts("------------");//画一个杠杠来分割每一次消元结果 for(int j=1;j<=m;j++){ for(int k=1;k<=n+1;k++)cout<<a[j][k]<<" ";//输出每次消元结果 puts("");//输出空格 } } if(flag){//判断无解和多解的情况(上面已经提到了) 已经懵bi求救qwq for(int i=1;i<=m;i++) if(!check(i)){//利用check判断是无解还是多解 printf("No solution\n"); return 0; } //每个答案行如果有一个等于0的数就多解?? printf("So many solutions\n"); } }
修改后(from rainbow):
#include<iostream> #include<algorithm> #include<cstdio> #include<cstring> #include<cmath> #include<cstdlib> using namespace std; typedef long long ll; typedef long double ld; typedef pair<int,int> pr; const double pi=acos(-1); #define rep(i,a,n) for(int i=a;i<=n;i++) #define per(i,n,a) for(int i=n;i>=a;i--) #define Rep(i,u) for(int i=head[u];i;i=Next[i]) #define clr(a) memset(a,0,sizeof a) #define pb push_back #define mp make_pair #define fi first #define sc second ld eps=1e-9; ll pp=1000000007; ll mo(ll a,ll pp){if(a>=0 && a<pp)return a;a%=pp;if(a<0)a+=pp;return a;} ll powmod(ll a,ll b,ll pp){ll ans=1;for(;b;b>>=1,a=mo(a*a,pp))if(b&1)ans=mo(ans*a,pp);return ans;} ll read(){ ll ans=0; char last=' ',ch=getchar(); while(ch<'0' || ch>'9')last=ch,ch=getchar(); while(ch>='0' && ch<='9')ans=ans*10+ch-'0',ch=getchar(); if(last=='-')ans=-ans; return ans; } //head int n,m; double a[100][100]; bool check(int k){ if(fabs(a[k][n+1])<eps)return 1; rep(i,1,n) if(fabs(a[k][i])>eps)return 1; return 0; } int main(){ n=read();m=n; // a_i,1 a_i,2 ... a_i,n a_i,n+1 rep(i,1,m) rep(j,1,n+1)a[i][j]=read(); int flag=0; rep(i,1,n){ int t=i; while(a[t][i]==0 && t<=n)t+=1; if(t==n+1){ flag=1; continue; } rep(j,1,n+1)swap(a[i][j],a[t][j]);//交换两行 double kk=a[i][i];//每一行对角线上的值 rep(j,1,n+1)a[i][j]/=kk; rep(j,1,m)//循环m个式子 开始消元 if(i!=j){ double kk=a[j][i]; rep(k,1,n+1) a[j][k]-=kk*a[i][k];//这样就能保证正好把第i列的数除了a[i][i] 都消成0 } } if(flag){ return printf("No Solution\n"),0; } rep(j,1,m){ printf("%.2lf",a[j][n+1]/a[j][j]); puts(""); } }
我们把对角线都整成1之后,为了消元,不能改变对角线的值
那么 i ! = j 确定了对角线不会被消掉
消元过程:
1. 把对角线消为1
比如这样: 1 x y z
a 1 b c
d g 1 h
p q j 1
2.消掉最外面一层,削成0
1 0 0 0
0 1 b c
0 g 1 h
0 q j 1
3.再消下一层
1 0 0 0
0 1 0 0
0 0 1 h
0 q j 1
4.以此类推吧
5.最后得到一个单位矩阵
特别鸣谢:rainbow Lz