.高斯消元

151 阅读 0 评论 100 点赞

我是靠谱客的博主成就冬日，这篇文章主要介绍.高斯消元，现在分享给大家，希望可以做个参考。

假如要解n元1次方程组，将n个方程放入一个n*n的矩阵，矩阵的每个位置代表未知数的系数，再写一个答案矩阵，记录每个方程的解，然后将矩阵用方程组的性质（加，乘）化为一个阶梯矩阵,如同下图

$begin{bmatrix} 1 & * & *\ 0 & 1& *\ 0& 0 &1 end{bmatrix}=begin{bmatrix} x\ y\z end{bmatrix}$

这时可以发现，最后一个方程只有最后一个未知数的系数为1，其余都为0，所以可以知道最后的未知数的值，然后用已知的未知数的值一步一步向上递推消元，就可以求出方程组的解

当然，方程组可能无解或有无限多个解

无解：有一行出现了形如0=d（d为非0数）的式子

无数解：有一行为0=0，这是说明有另外一行是ax+by=d的形式，故此时方程有无数解

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include<bits/stdc++.h>
using namespace std;
int n;
double a[110][110],ans[110];
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++) scanf("%lf",&a[i][j]);
scanf("%lf",&ans[i]);
}
for(int i=1;i<=n;i++){
double mx=-1000000.00;
int id;
for(int j=i;j<=n;j++){
if(a[j][i]>mx){
mx=a[j][i];
id=j;
}
}//a矩阵记录了方程每个未知数的系数，ans记录了方程组的答案
for(int j=1;j<=n;j++) swap(a[i][j],a[id][j]);
swap(ans[i],ans[id]);
if(!a[i][i]){
if(ans[i]) printf("No Solution");//出现了0=d这种形式
else printf("wushujie");//出现了0=0这种形式
return 0;
}
ans[i]/=a[i][i];
for(int j=n;j>=i;j--) a[i][j]/=a[i][i];
for(int j=i+1;j<=n;j++){
double ty=a[j][i];
for(int k=1;k<=n;k++) a[j][k]-=ty*a[i][k];
ans[j]-=ty*ans[i];
}
}
for(int i=1;i<=n;i++){
bool ok=false;
for(int j=1;j<=n;j++){
if(a[i][j]) ok=true;
}
if(!ok){
printf("No Solution");
return 0;
}
}
for(int i=n-1;i>=1;i--){
for(int j=i+1;j<=n;j++) ans[i]-=ans[j]*a[i][j];//递推答案
}
for(int i=1;i<=n;i++){
printf("%.2lfn",ans[i]);
}
return 0;
}

例题

一个球面上的点到圆心的距离都相等，所以就有

$forall_{i} sum_{j=1}^{jleq n}(a_{ij}-x_{j})^{2}=d$

a[i][j]表示球面上第i个点第j维的坐标，x[j]表示圆心第j维的坐标，d是一个固定的数

因为d是固定的，并且有有n+1个这种式子，所以可以让这些式子相减得到n个为0的式子

$because sum_{j=1}^{jleq n}(a_{i,j}-x_{j})^{2}-(a_{i,j+1}-x_{j+1})=0$

$therefore sum_{j=1}^{n}a_{i,j}^{2}+2*x_{j}*a_{i,j}+x_{j}^{2}-a_{i,j+1}^{2}+2*x_{j+1}*a_{i,j+1}-x_{j+1}^{2}=0$

$therefore sum_{j=1}^{jleq n} 2*x_{j}*(a_{i,j}-a_{i,j+1})=sum_{j=1}^{jleq n} a_{i,j}^{2}-a_{i,j+1}^{2}$

这就是高斯消元的形式了，第i行第j位的值为2*x[j]*(a[i][j]-a[i][j+1])，答案的值为a[i][j]的平方减a[i][j+1]的平方的差的和

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include<bits/stdc++.h>
using namespace std;
int n;
double k[110][110],a[110][110],ans[110];
int main(){
scanf("%d",&n);
for(int i=1;i<=n+1;i++){
for(int j=1;j<=n;j++){
scanf("%lf",&a[i][j]);
}
}
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
k[i][j]=2.00*(a[i][j]-a[i+1][j]);//每一位的系数
ans[i]=ans[i]+a[i][j]*a[i][j]-a[i+1][j]*a[i+1][j];
}
}
for(int i=1;i<=n;i++){
double mx=-10000000.00;
int id;
for(int j=i;j<=n;j++){
if(mx<k[j][i]){
mx=k[j][i];
id=j;
}
}
for(int j=1;j<=n;j++) swap(k[i][j],k[id][j]);
swap(ans[i],ans[id]);
ans[i]/=k[i][i];
for(int j=n;j>=i;j--) k[i][j]/=k[i][i];
for(int j=i+1;j<=n;j++){
double bs=k[j][i];
for(int h=1;h<=n;h++) k[j][h]-=bs*k[i][h];
ans[j]-=bs*ans[i];
}
}
for(int i=n-1;i>=1;i--){
for(int j=i+1;j<=n;j++){
ans[i]-=ans[j]*k[i][j];
}
}
for(int i=1;i<=n;i++){
printf("%.3lf ",ans[i]);
}
return 0;
}

异或运算的高斯消元

比如说有3个未知数和几个式子，那么矩阵每一位就表示在这个式子中xi是取还是不取

可以将矩阵压成一个数，这样可以缩小时间复杂度

每次找到最大的数来更新

和普通高斯消元一样，使第i位为1，1~i-1位为0

题面上说对2取模，相当于就是这一堆数异或后最后一位的值

n达到了1000，只能用bitset来表示每一位的状态

对于第i位，找到第一个第i位为1的式子来更新答案

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
/*bitset对高斯消元的优化
对于只有异或运算的高斯消元，可以用bitset将状态压为一个数，使得时间复杂度降倒1/32
找最大值时候就枚举要找的位数就行
*/
#include<bits/stdc++.h>
using namespace std;
bitset<1010>a[2010];
int m,n,c,ans,pd[2010];
int main(){
scanf("%d%d",&m,&n);
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++){
char ch=getchar();
while(ch!='0'&&ch!='1') ch=getchar();
a[i][j]=ch-'0';
}
int w;
scanf("%d",&w);
a[i][m+1]=w;
}
for(int i=1;i<=m;i++){
int j=i;
while(!a[j][i]&&j<=n) j++;
if(j>n){
printf("Cannot Determine");
return 0;
}
ans=max(ans,j);
swap(a[i],a[j]);
for(int j=1;j<=n;j++){
if(i!=j&&a[j][i]) a[j]^=a[i];
}
}
for(int i=m;i>=1;i--){
for(int j=i+1;j<=m;j++){
if(a[i][j]) a[i]^=a[j];
}
pd[i]=a[i][m+1];
}
printf("%dn",ans);
for(int i=1;i<=m;i++){
pd[i]%2?printf("?y7M#n"):printf("Earthn");
}
return 0;
}