高斯消元
高斯消元
作用：把方程组化为上三角。类似下图
\(\left [ \begin{array}{} a & a & a \\ 0 & a & a \\ 0 & 0 & a \end{array}\right ]\)
退去华丽的外表，其实就是一个加减消元
步骤：
-
选择第 \(i\) 项系数尽可能大的，对其他行消元
-
把最大系数所在行换到当前行，系数变为 \(a_{i,i}\)
-
当前行是 \(i\) ，消去行是 \(j\) ，设当前行需要乘 \(tmp=a_{j,i}/a_{i,i}\)
则每一个 \(k\ge i,a_{j,k}\) 要减去 \(a_{i,k}*tmp\)
-
回代：从最后一行往前面代入，求得解
无解判定：消元时会用到除法，若除数为 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\) 的逆矩阵
思路：
- 把 \(I\) 放在 \(A\) 的右边
- 对 \(A\) 进行消元使得 \(A\) 称为单位矩阵，可直接使用约旦-高斯消元
- 原来的单位矩阵转换成逆矩阵
因为：\(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());
}
要用到 ↓
辗转相除高斯消元
因为如上题，由于精度问题和逆元不一定存在，
普通高斯消元会出现无法解决的问题。
因为可以任意相减，于是想到辗转相除，这样一定能消掉一个
可用于矩阵树定理
