高斯消元

高斯消元

作用：把方程组化为上三角。类似下图

\(\left [ \begin{array}{} a & a & a \\ 0 & a & a \\ 0 & 0 & a \end{array}\right ]\)

退去华丽的外表，其实就是一个加减消元

步骤：

1. 选择第 \(i\) 项系数尽可能大的，对其他行消元

2. 把最大系数所在行换到当前行，系数变为 \(a_{i,i}\)

3. 当前行是 \(i\) ，消去行是 \(j\) ，设当前行需要乘 \(tmp=a_{j,i}/a_{i,i}\)

   则每一个 \(k\ge i,a_{j,k}\) 要减去 \(a_{i,k}*tmp\)

4. 回代：从最后一行往前面代入，求得解

无解判定：消元时会用到除法，若除数为 0，则无解

例：一个模板

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
int n,m;
double x[105][105],ans[105];
int main() {
	scanf("%d",&n),m=n+1;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=m;j++)
			scanf("%lf",&x[i][j]);
	for(int i=1,r;i<=n;i++) {
		r=i;
		for(int j=i+1;j<=n;j++)
			if(fabs(x[r][i]-x[j][i])<=eps)
				r=j;
		if(fabs(x[i][i])<=eps)
			return printf("No Solution"),0;
		if(r^i)swap(x[i],x[r]);
		double tmp;
		for(int j=i+1;j<=n;j++) {
			tmp=x[j][i]/x[i][i];
			for(int k=i;k<=m;k++)
				x[j][k]-=x[i][k]*tmp;
		}
	}
	for(int i=n;i;i--) {
		for(int j=i+1;j<=n;j++)
			x[i][m]-=x[j][m]*x[i][j];
		x[i][m]/=x[i][i];
	}
	for(int i=1;i<=n;i++)
		printf("%.2lf\n",x[i][m]);
} 

约旦-高斯消元

精度好像更优秀，还不用回代

作用：把一个方程组化成对三角。类似下图

\(\left [ \begin{array}{} a & 0 & 0 \\ 0 & a & 0 \\ 0 & 0 & a \end{array}\right ]\)

其步骤和普通高斯消元类似

就是第 3 步消的元不再是 \(k\ge i,a_{j,k}\) 而是 \(k\ne i,a_{j,k}\)

求解：直接用唯一剩下的系数求就可以了

例：还是一个模板

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
int n,m;
double x[105][105];
int main() {
	scanf("%d",&n),m=n+1;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=m;j++)
			scanf("%lf",&x[i][j]);
	for(int i=1,r;i<=n;i++) {
        // 找到最大的系数
		r=i;
		for(int j=i+1;j<=n;j++)
			if(fabs(x[r][i]-x[j][i])<=eps)
				r=j;
        // 无解判定
		if(fabs(x[i][i])<=eps)
			return printf("No Solution"),0;
		if(r^i)swap(x[i],x[r]);
        // 加减消元
		double tmp;
		for(int j=1;j<=n;j++)
			if(i^j) {
				tmp=x[j][i]/x[i][i];
				for(int k=i+1;k<=m;k++)
					x[j][k]-=x[i][k]*tmp;
			}
	}
	for(int i=1;i<=n;i++)
	x[i][m]/=x[i][i];
	for(int i=1;i<=n;i++)
		printf("%.2lf\n",x[i][m]);
} 

逆矩阵

定义：假设 \(A\) 是一个方阵，如果存在 \(A^{-1}\) 使得 \(A^{-1}A=I\) 那么矩阵 \(A\) 可逆，\(A^{-1}\) 称为 \(A\) 的逆矩阵

思路：

1. 把 \(I\) 放在 \(A\) 的右边
2. 对 \(A\) 进行消元使得 \(A\) 称为单位矩阵，可直接使用约旦-高斯消元
3. 原来的单位矩阵转换成逆矩阵

因为：\(A^{-1}*[AI]=[IA^{-1}]\)

例：矩阵求逆。注意求逆元即可。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod=1000000007;
inline LL Pow(LL x,LL y) {
	register LL res=1;
	for(;y;y>>=1,x=x*x%mod)
	if(y&1)res=res*x%mod;
	return res;
}
int n,m;
LL x[405][805];
int main() {
	scanf("%d",&n),m=n<<1;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			scanf("%lld",&x[i][j]);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			x[i][j+n]=(i==j);
	for(int i=1,r;i<=n;i++) {
		r=i;
		for(int j=i+1;j<=n;j++)
			if(x[j][i]>x[r][i])r=j;
		if(!x[i][i])return puts("No Solution"),0;
		if(r^i)swap(x[i],x[r]);
		register LL inv=Pow(x[i][i],mod-2),tmp;
		for(int j=1;j<=n;j++)
			if(i^j) {
				tmp=x[j][i]*inv%mod;
				for(int k=i+1;k<=m;k++)
					x[j][k]=(x[j][k]-(x[i][k]*tmp%mod)+mod)%mod;
			}
		for(int j=1;j<=m;j++)
			x[i][j]=x[i][j]*inv%mod;
	}
	for(int i=1;i<=n;i++,putchar(10))
		for(int j=1;j<=n;j++)printf("%lld ",x[i][j+n]);
}

求解行列式

\(\det(A)=\begin{vmatrix}A\end{vmatrix}=\sum_p (-1)^{\tau(p)}\prod_{i=1}^n a_{i,p_i}\)

其中 \(p\) 是排列， \(\tau(p)\) 是排列 \(p\) 的逆序对数

行列式性质：

交换任意两行（列）结果变为相反数，

一行加上另一行乘常数结果不变

一行同乘 \(k\) 结果也乘 \(k\)

这些性质和高斯消元十分类似，所以行列式求值用高斯消元

例：一个以下代码无法通过的题

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
int n,m;
double x[105][105],ans[105];
inline double gauss() {
	register double res=1;
	for(int i=1,r;i<=n;i++) {
		r=i;
		for(int j=i+1;j<=n;j++)
			if(fabs(x[r][i]-x[j][i])<=eps)
				r=j;
		if(fabs(x[i][i])<=eps)
			return 0;
		if(r^i)swap(x[i],x[r]),res=-res;
		res*=x[i][i];
		double tmp;
		for(int j=i+1;j<=n;j++) {
			tmp=x[j][i]/x[i][i];
			for(int k=i;k<=m;k++)
				x[j][k]-=x[i][k]*tmp;
		}
	}
	return res;	
} 
int main() {
	scanf("%d",&n),m=n;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=m;j++)
			scanf("%lf",&x[i][j]);
	printf("%.2lf\n",gauss());
} 

要用到 ↓

辗转相除高斯消元

因为如上题，由于精度问题和逆元不一定存在，

普通高斯消元会出现无法解决的问题。

因为可以任意相减，于是想到辗转相除，这样一定能消掉一个

可用于矩阵树定理