高斯消元(求解n元一次方程组)

最新推荐文章于 2023-11-06 11:13:39 发布

是哆啦D梦

最新推荐文章于 2023-11-06 11:13:39 发布

阅读量3.6k

点赞数 4

分类专栏：数学

本文链接：https://blog.csdn.net/weixin_43602607/article/details/111935577

版权

数学专栏收录该内容

8 篇文章

订阅专栏

一.高斯消元求解线性方程组

1. 问题： $m\times(n+1)$ 描述一个 n元线性方程组方程式，Rij为第 i 个方程未知数 xj 的系数，求解未知数 x1->xn

~~- 如果是我们人，我们自然会消元法来消元，通过将某个多元一次方程化为一元一次方程，求出结果。~~

2.方法：高斯消元的转化形式

我们定义第i行(第i个方程)为Ri，那么消元的过程就相当于以下三种操作：
1.交换Ri与Rj
2.Ri乘以常数k
3.Ri加上(Rj乘以常数k) ~~其实就是线性代数中的初等行变化~~
高斯消元的思路就是化为以下形式的矩阵(首非0元为1的阶梯矩阵)

$\begin{bmatrix} 1 & a & b &c\\ 0 &1 &d &e \\ 0 & 0 &1 &g \\ \end{bmatrix}$

这样就可以倒着回代求出所有解了

3.高斯-约旦消元

$\begin{bmatrix} 1 & 0 & 0 &a\\ 0 &1 &0 &b \\ 0 & 0 &1 &c \\ \end{bmatrix}$

要是能转化成这种形式，那就能避免回带了

4.如何转化成该形式？？？

步骤如下：

0.用 m 个方程解 n 个未知数，方程下标为 $[1, m]$ ；解下标为 $[1, n]$ ，n+1为等号右边数字的下标
1.一列一列处理，在用第 h 个方程处理第 i 个未知数
2.在未处理行中找到一个第 i 列为非 0 元素的行，并将其置换到第 h 行
3.将第 h 行第 i 列元素更改为 1 ，并修改行中的其他元素
4.用第 h 行第 i 列的 1 把其他行的第 i 列非 0 元素均变为 0
5.234成功后，表示第 h 个方程成功处理了第 i 个未知数，则 h++
6.结束后，表示用了 h-1 个方程解 n 个未知数。

注意事项

【注意】由于过程中会损失精度，所以我们通常会给定一个精度eps，当元素值小于eps时，将其视为 0

解的情况（无解，多解，自由元）

问：求 n 元解需要多少个线性方程呢？
答：最多需要 n 个，但不是前 n 个方程，而是置换后的前 n 个。倘若处理完 n 个未知数了，而只使用了 h-1（h-1<n）个方程，那么该方程组可能要出现无解或者多解了。
问：什么时候只使用了 h-1（h-1<n）个方程？
答：步骤 2 在未处理行中找不到一个第 i 列元素为非 0 的行。
那么它就会变成类似下面这个式子：（可以看到在处理 x3 时，没有找到一个非0行与其置换）
$\begin{bmatrix} 1 & 0 & e &0&a\\ 0 &1 &d &0&b \\ 0 & 0 &0 &1 &c \\ 0 &0 &0 &0 &p \end{bmatrix}$
可以看出，从第 h 个方程开始的所有方程的系数都会变为 0 。
问：什么时候无解？什么时候有解？
答： p!=0 时无解，p=0时有多解。我们可以去判断 $[h, n]$ 的所有方程的右边数字是否恒等于 0 来判断无解还是多解。
问：什么是自由元？
答：当出现多解时候，就会出现自由元。即：有些未知数的取值可以任意取，其取值可能影响着其他非自由元的取值也可能不影响。
问：有多少个自由元？是哪几个？
答： n-(h-1)个，是那些无法在步骤 2 中找到非 0 元素的未知数。

时间复杂度： O(n³)

4.模板代码

m个n元线性方程组（无解返回0，多解返回-1，唯一解返回1）

#include<bits/stdc++.h>
using namespace std;
double a[105][105];
double eps=1e-6;
int m,n;
int gauss() {
	int h=1;
	for(int i=1; i<=n; i++) { //处理xi 
		int maxj=h;
		for(int j=h+1; j<=m; j++) { //从未处理行中找到一个xi系数不为0的方程 
			if(fabs(a[maxj][i])<fabs(a[j][i]))maxj=j;
		}
		if(fabs(a[maxj][i])<eps)continue;//无法确定xi的值，去考虑下一个未知数 
		if(maxj!=h)swap(a[maxj],a[h]);//置换到第 h 个方程来 
		
		for(int j=i; j<=n+1; j++)a[h][j]/=a[h][i]; //将xi系数变为1，同时修改其所在方程其他系数的值
		for(int k=1; k<=m; k++) { //将其他方程的xi的系数全变为0，同时修改对于方程其他未知数的系数的值 
			if(k==h)continue;
			double div=a[k][i];
			for(int j=i; j<=n+1; j++) {
				a[k][j]-=div*a[i][j];
			}
		}
		h++; 
	}
	if(h<=n){
		for(int i=h;i<=n;i++){//有0就无解了 
			if(a[i][n+1]<eps)return 0;
		}
		return -1;//都不是0，有多解（自由元数量为n-h+1） 
	}else return 1;
}
int main() {
	cin>>n;
	for(int i=1; i<=m; i++) {
		for(int j=1; j<=n+1; j++)cin>>a[i][j];
	}

	int t=gauss();
	if(t==1) {
		for(int i=1; i<=n; i++)printf("%.2lf\n",a[i][n+1]);
	} else if(t==0)printf("无解");
	else if(t==-1)printf("有多解");
}

二.高斯消元求解异或方程式

1.问题： $m\times(n+1)$ 描述一个 n元异或方程组，Rij (0或1) 为第 i 个方程未知数 xj 的系数，求解未知数 x1->xn

2.方法：高斯-约旦消元

与一般的高斯消元相似，只不过加变成了异或

3.bitset优化

由于是行与行之间的所有元素相互异或，我们我们可以直接使用 bitset O(1) 代替 O(n)
由于方程的系数是0或1，并且是由方程之间的相互异或操作，我们可以直接使用 bitset 来表示一整个方程系数。这样异或操作就从O(1)变成了O(n)。
时间复杂度： $O (n m)$

#include<bits/stdc++.h>
using namespace std;
const int N=2005;
bitset<2005>a[N];
int m,n;
int gauss() {
	int h=1;
	for(int i=1; i<=n; i++) { //处理xi 
		int pos=-1;
		for(int j=h; j<=m; j++) { //从未处理行中找到一个xi系数不为0的方程 
		    if(a[j][i]){
		    	pos=j;
		    	break;
			}
		}
		if(pos==-1)continue;//无法确定xi的值，去考虑下一个未知数 
		if(pos!=h)swap(a[pos],a[h]);//置换到第 h 个方程来 
		
		for(int j=1; j<=m; j++) { //将其他方程的xi的系数全变为0，同时修改对于方程其他未知数的系数的值 
			if(a[j][i]==1&&j!=h)a[j]^=a[h];
		}
		h++; 
	}
	if(h<=n){
		for(int i=h;i<=n;i++){//有0就无解了 
			if(a[i][n+1]==0)return 0;
		}
		return -1;//都不是0，有多解（自由元数量为n-h+1） 
	}else return 1;
}
int main() {
	cin>>n;
	int x;
	for(int i=1; i<=m; i++) {
		for(int j=1;j<=n+1;j++){
			scanf("%1d",&x);
		    a[i][j]=x;
		}
	}

	int t=gauss();
	if(t==1) {
		for(int i=1; i<=n; i++)printf("%d\n",a[i][n+1]);
	} else if(t==0)printf("无解");
	else if(t==-1)printf("有多解");
}

三.高斯消元求逆矩阵

1.问题：给定n阶方阵，判断矩阵是否可逆，并求其 mod p下的逆矩阵

2.方法：高斯消元

仍然用初等行变化进行求解
设求解的矩阵为A，单位矩阵E。结论： (A,E)->(E,A^-1）
那么问题就转化为，如何将一个A转化为E

步骤如下：

1.一行一行处理，(假设处理到第i行)
2.通过交换两行使得第i行首元非0
3.将首非0元素变为1,(j就是第i行的第i个元素),并处理该行其他元素，ps:因为这涉及除法，所以要求逆元
4.用第i行去将该“1”所在其他列的其他元素消为0
【注意】若在求解的过程中无法实现第二步，则矩阵不可逆

矩阵最后会变成这个样子，右边就是该矩阵的逆矩阵

$\begin{bmatrix} 1 & 0 & 0 &a &b &c\\ 0 &1 &0 &d &e &f \\ 0 & 0 &1 &g &h &i\\ \end{bmatrix}$

3.模板代码

n阶行列式

#include<bits/stdc++.h>
using namespace std;

const int N=405;
const long long p=1e9+7;
int n,m;
long long a[N][N<<1];
long long ksm(long long a,long long b) { //求逆元
	ret=1;
	while(b) {
		if(b&1)ret=ret*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ret;
}

int main() {
	cin>>n;
	m=n*2;
	for(int i=1; i<=n; ++i) {
		for(int j=1; j<=n; j++) {
			cin>>a[i][j];
		}
		a[i][n+i]=1;//矩阵右边构造一个单位矩阵
	}
	for(int i=1; i<=n; i++) { //高斯消元板子
		for(int k=i; k<=n; k++) {//保证首非0元素不为0 
			if(a[k][i]) {
				for(int j=1; j<=m; j++)swap(a[i][j],a[k][j]);
				break;
			}
		}
		if(!a[i][i]) {//若首非0元素为0，则矩阵无解 
			cout<<"No Solution";
			return 0;
		}
		long long x=ksm(a[i][i],p-2);  //求逆元
		for(int j=i; j<=m; j++)a[i][j]=a[i][j]*r%p;//更改当前行
		for(int k=1; k<=n; k++) { //更改其他行信息
			if(k!=i) {
				long long div=a[k][i];
				for(int j=i; j<=m; j++)a[k][j]=(a[k][j]-div*a[i][j]%p+p)%p;
			}
		}
	//最后的矩阵的样子大概如下
	//100abc
	//010def
	//001ghi 
	}
	for(int i=1; i<=n; i++) {
		for(int j=n+1; j<=m; j++){
			cout<<a[i][j]<<" ";
		}
		cout<<endl;
	}

	return 0;
}

例题汇总

异或方程组

例题1

题目描述： 有 n 只虫子，每只虫子要么一只腿要么两只腿。有 m 次操作，每次操作取若干只虫子放入瓶中，瓶子会返回1或0来表示一共有奇数条腿还是偶数条腿。用 a[i][j] 描述第 j 只虫子取与不取，取为1不取为0 。求最少前需要前多少次操作，才能判断出每只虫子是一条腿还是两条腿。如果判断不了，就输出无法确定。
问题分析： 相加对 2 取余刚好对应了异或操作，即变成了最少需要前多少次操作，才能解除这个异或方程组。我们在置换的时候维护一下用到的方程的最大下标即可。

#include<bits/stdc++.h>
using namespace std;
const int N = 2005;
bitset<2005>a[N];
string s;
int ans=0,n,m;

int gssy() {
	int h=1;
	for(int i=1; i<=n; i++) { //处理xi
		int pos=-1;
		for(int j=h; j<=m; j++) { //从未处理行中找到一个xi系数不为0的方程
			if(a[j][i]) {
				pos=j;
				break;
			}
		}
		if(pos==-1)return 0;//无法确定xi的值，去考虑下一个未知数
		ans=max(ans,pos);
		if(pos!=h)swap(a[pos],a[h]);//置换到第 h 个方程来

		for(int j=1; j<=m; j++) { //将其他方程的xi的系数全变为0，同时修改对于方程其他未知数的系数的值
			if(a[j][i]==1&&j!=h)a[j]^=a[h];
		}
		h++;
	}
	return 1;
}
int main() {
	cin>>n>>m;
	int x;
	for(int i=1; i<=m; i++) {
		for(int j=1; j<=n+1; j++) {
			scanf("%1d",&x);
			a[i][j]=x;
		}
	}
	if(gssy()==0)cout<<"Cannot Determine";
	else {
		cout<<ans<<endl;
		for(int i=1; i<=n; i++) {
			if(a[i][n+1]==1)cout<<"?y7M#"<<endl;
			else cout<<"Earth"<<endl;
		}
	}
	return 0;
}