高斯消元学习笔记

高斯消元学习笔记

——by sunzz3183

---

介绍

高斯消元是一种求解线性方程组的方法。线性方程组就是 $m$ 个 $n$ 元一次方程。如：

$$\left \{\begin{matrix} x_1+2x_2-x_3&=-6 \\2x_1+x_2-3x_3&=-9 \\-x_1-x_2+2x_3&=7 \end{matrix} \right.$$

求解

线性方程组可以写成一个由系数和常数组成的 $m$ 行 $n+1$ 列的增广矩阵。如：

$$\begin{bmatrix} 1& 2& -1& -6\\ 2& 1& -3& -9\\ -1& -1& 2& 7 \end{bmatrix}$$

增广矩阵有三个性质：

1. 任意交换矩阵的两行或两列，矩阵不变；

2. 矩阵一行乘上一个非零数，矩阵不变；

3. 把其中一行的若干倍加到另一行上，矩阵不变。

那我们运用矩阵的各种性质，来将矩阵消成对角线上的元素为 $1$，并且除了第 $n+1$ 列其余元素均为 $0$ 的矩阵，这样我们就很容易的得出每个未知数的值：分别是从上到下第 $n+1$ 列的值（因为这时候第 $i$ 行的第 $i$ 个未知数系数都为 $1$）。

那么怎么消呢？

高斯约旦消元法

以上面的矩阵为例子：

$$\begin{bmatrix} 1& 2& -1& -6\\ 2& 1& -3& -9\\ -1& -1& 2& 7 \end{bmatrix}$$

明确我们的目的：把矩阵消成对角线为 $1$，除了第 $n+1$ 列其余元素都为 $0$。

若有一列全为 $0$，那么肯定有第 $i$ 行第 $i$ 列消不成 $1$，此时无解（少了一个方程，未知数就无法确定了）。

知道了这个，我们就可以对这个矩阵进行初步判定：

for(int i=1,j,l;i<=n;i++){
    for(j=i+1,t=i;j<=n;j++)
        if(fabs(a[j][i])>fabs(a[t][i]))t=j;//选绝对值最大的作用一会讲
    for(j=1;j<=n+1;j++)swap(a[i][j],a[t][j]);
    if(fabs(a[i][i])<eps)return puts("No Solution"),0;    
}

然后利用性质 $2$ 让第 $i$ 行第 $i$ 列系数化为 $1$。

for(j=1,k=a[i][i];j<=n+1;j++)a[i][j]/=k;

最后，利用性质 $3$ 将其他行的系数化为 $0$

for(j=1;j<=n;j++)
    for(l=1,k=a[j][i];i!=j&&l<=n+1;l++)
        a[j][l]-=k*a[i][l];

然后就完成了！

最后第 $n+1$ 列就是答案了。

完整代码

#include<bits/stdc++.h>
// #define int long long
#define eps 0.000000001
using namespace std;
inline int read(){
    char ch=getchar();int x=0;bool f=1;
    while(ch<'0'||'9'<ch){if(ch=='-')f=0;ch=getchar();}
    while('0'<=ch&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
    return f?x:-x;
}
const int N=1e2+3;
int n,t=1;
double a[N][N],k;
signed main(){
    // freopen(".in","r",stdin);
    // freopen(".out","w",stdout);
    n=read();
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&a[i][j]);
    for(int i=1,j,l;i<=n;i++){
        for(j=i+1,t=i;j<=n;j++)
            if(fabs(a[j][i])>fabs(a[t][i]))t=j;
        if(fabs(a[t][i])<eps)return puts("No Solution"),0;
        for(j=1;j<=n+1;j++)swap(a[i][j],a[t][j]);
        for(j=1,k=a[i][i];j<=n+1;j++)a[i][j]/=k;
        for(j=1;j<=n;j++)
            for(l=1,k=a[j][i];i!=j&&l<=n+1;l++)
                a[j][l]-=k*a[i][l];
    }
    for(int i=1;i<=n;i++)
        printf("%.2f\n",a[i][n+1]);
    return 0;
}

遗留问题

$Q$：为什么要选绝对值最大的数？

$A$：在消元过程中，为了避免出现除数为 $0$ 的情况，需要在每一行中找到一个绝对值最大的非零数作为主元，然后将该行与其他行进行消元。这样可以避免出现较大的误差，提高计算精度。IEEE 754 浮点数标准中数越接近 $0$ 越精确。

例题

1.模板

【模板】高斯消元法

数据较水，正常就能过。

2. 进阶

[SDOI2006]线性方程组

区别：数据更强，并加入了判断无解和多解。

判断方法：如果存在一行全为 $0$ 显然就是多解，如果存在一行系数全为 $0$ 常数不为 $0$，那么是无解。

注意：无解的优先级比多解高，多解的优先级比唯一解高。

3. 异或消元

异或就是不进位的加法和减法，正常做就行。

优化：可以用一个数表示二进制数字。

代码：

#include<bits/stdc++.h>
// #define int long long
using namespace std;
inline int read(){
    char ch=getchar();int x=0;bool f=1;
    while(ch<'0'||'9'<ch){if(ch=='-')f=0;ch=getchar();}
    while('0'<=ch&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
    return f?x:-x;
}
const int N=31;
int T,n,t=1,ans;
int a[N];
int yudan(){
    for(int i=1,j,k;i<=n;i++){
        for(j=i+1,t=i;j<=n;j++)
            if(a[j]>a[t])t=j;//找主元最大的，实际上就是找a[j][i]==1
        swap(a[i],a[t]);
        if(!a[i])return 1<<(n-i+1);
        if(a[i]==1)return 0;
        k=floor(log2(a[i]));
        for(j=1;(a[i]>>k&1)&&j<=n;j++)
            if(j!=i&&(a[j]>>k&1))a[j]^=a[i];
    }
    return 1;
}
signed main(){
    // freopen(".in","r",stdin);
    // freopen(".out","w",stdout);
    T=read();
    while(T--){
        n=read();
        for(int i=1;i<=n;i++)a[i]=read(),a[i]|=1<<i;
        for(int i=1;i<=n;i++)a[i]^=read();
        int u=read(),v=read();
        while(u&&v)
            a[v]|=1<<u,u=read(),v=read();
        int ans=yudan();
        if(!ans)puts("Oh,it's impossible~!!");
        else printf("%lld\n",ans);
    }
    return 0;
}