高斯消元

高斯消元¶

高斯消元法（Gaussian Elimination）是一种线性代数中的方法，用于求解变量数量与方程数量相等的线性方程组。它通过一系列的代数运算，将方程组转换为一个简单的方程，然后求解这个简单的方程来得到原始方程组的解。

高斯消元法的步骤如下：

将方程组写成增广矩阵（又称扩增矩阵），即在原方程组的右边添上一列，这一列是每个方程右边的和。
从左到右，从上到下，依次消去增广矩阵中的每一列，使得最下面的行成为仅含有一个未知数的方程。
重复第二步的操作，直到所有的未知数都出现在最下面的行中。
最后，将最下面的行中的系数除以-1，以求得最后一个未知数的值。

533 线性方程组高斯消元法_哔哩哔哩_bilibili

www.luogu.com.cn

例题 #1 求解线性方程组¶

给定一个线性方程组，对其求解

输入格式

第一行，一个正整数 n

第二至 n+1 行，每行 n+1 个整数，为a1,a2⋯an 和 b，代表一组方程。

输出格式

共 n 行，每行一个数，第 i 行为 xi （保留 2 位小数）

如果不存在唯一解，在第一行输出 No Solution.

增广矩阵¶

我们对方程组进行以上操作。由方程组特性可以得到以下变换结论：

矩阵的初等行变换，具体内容为： 1.交换两行 2.把某一行乘一个非0的数 3.把某行的若干倍加到另一行上去

以上三个操作都不会对方程组的解产生影响

高斯消元法¶

消元，懂吧？

e.g.

$\begin{cases}a_1x+b_1y=c_1\\a_2x+b_2y=c_2\end{cases}$

针对以上，根据初等行变换，我们可以把下式中的$a_2x$消去（用$下式-n \times 上式$ 即可）

递推到$n$个方程的方程组，也是一样的。虽然我们可能不能将两个方程消到单项式个数相同，但每个之间相差一个是一定可以做到的。

进一步地，我们还可以将消去后的方程组的每一个方程的第一个单项式（我们把这个单项式里的未知数$x$作为这个方程的主元）的系数化为$1$.然后从下往上，依次求解带入上面一个方程，就完工啦！

概括一下，就是

找
换
主元系数变为1
其他方程的该系数变为0

下面就先举例解释：

图

初始状态
从第i=1行问下，寻找第一个“主元$x_1$（第i=1个单项式）系数不为0的方程”并且将它交换到第i=1行
将主元系化为1
将第i=1行下方的方程（或者是说未**处理过的方程**）如果有$x_1$的系数不为0的，就用**变换3**将其系数变为0
i=2，重复以上2~4步操作
...
...所有方程组都处理完毕后，从下往上，逐步计算出$x_3,x_2,x_1$的值（从下往上代入）

代码¶

代码如下

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

const double eps=1e-6;
const int N=1e3+5;
double a[N][N];
int n;

bool guass(){
    for(int i=1;i<=n;i++){
        int r=i;

        for(int j=i+1;j<=n;j++)
            if(fabs(a[r][i])<fabs(a[j][i]))r=j;// fabs不能省

      /*以下被注释代码是错误的，无法处理
        1 1 2
        3 3 6
      */

        // for(int k=i;k<=n;k++){
        //  if(a[k][i]!=0){
        //      r=k;
        //      break;
        //  }
        // }

        if(r!=i)swap(a[r],a[i]);
        if(a[i][i]-0.0==0)return 0;

        for(int j=n+1;j>=i;j--){
            a[i][j]/=a[i][i];
        }

        for(int j=i+1;j<=n;j++){
            if(a[j][i]==0)continue;
            for(int k=n+1;k>=i;k--)a[j][k]-=a[j][i]*a[i][k];
        }
      //可以直观感受下流程：
        // for(int j=1;j<=n+1;j++)cerr<<a[i][j]<<' ';
        // cerr<<endl;

    }


    for(int i=n-1;i;i--){
        for(int j=i+1;j<=n;j++){
            a[i][n+1]-=a[i][j]*a[j][n+1];
        }
    }

    return 1;
}

signed main(){
    cin>>n;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n+1;j++){
            cin>>a[i][j];
        }
    }


    if(!guass())puts("No Solution");
    else {
        for(int i=1;i<=n;i++){
            printf("%.2lf\n",a[i][n+1]);
        }
    }
}

难点解析

for(int k=n=1;j>=i;k--)a[j][k]-=a[j][i]*a[i][k];//变其它数字的系数为0

a[j][k]-=a[j][i]*a[i][k]什么意思??

考察以下方程组

$\begin{cases}x+b_1y=c_1\\ax+b_2y=c_2\end{cases}$

我们怎么样将方程2的第一个单项式系数化为0呢?

当然是将它减去 $n$$ \times 方程1$ 了

那么这个比$n$是多少?

恰恰是a[j][i],即方程组中的$a$（注意，这是因为此时第i=1个方程组的主元系数已经化为了1）

那么对于**方程2**，我们将$a$减去$1\times a$ ，将$b_2$减去$ b_1 \times a $，将$c_2$减去$c_1 \times a$，就相当于将**方程2**减去了$a \times 方程1$，满足变换3

高斯-约旦消元法¶

模拟一下过程

图

和前面的一样
和前面的一样，但此处我们**不需要**将主元系数化为1了
i=2，用**方程2**将其他方程的第i=2个单项式的系数化为0（此处凑巧将第3个单项式也化为了0，纯属巧合）.因为在第2步已经把**方程2**的第一项系数化为0，因此不会对已经处理好的**方程1**的**主元项**造成影响。同样，我们**不需要**将主元系数化为1
同3
将$c_i$除以主元的系数，得到的新$c_i$就是方程组的解$x_i$

代码¶

注释：fabs()是对于小数的abs()

/*                                                                                
                      Keyblinds Guide
                    ###################
      @Ntsc 2024

      - Ctrl+Alt+G then P : Enter luogu problem details
      - Ctrl+Alt+B : Run all cases in CPH
      - ctrl+D : choose this and dump to the next
      - ctrl+Shift+L : choose all like this
      - ctrl+K then ctrl+W: close all
      - Alt+la/ra : move mouse to preHash/nxt pos'

*/

//头文件
#include <bits/stdc++.h>

//数据类型
#define ll long long
#define ull unsigned long long
#define db double
#define endl '\n'
//命名空间
using namespace std;
//常量
const int N = 2e5 + 5;
const int M = 1e3;
const int MOD = 903250223;
const int INF = 1e9;
//变量//变量
long double a[M][M];
int n,b,c,x[N],y[N],ans,res,tmp,cnt,m;

int gauss_jordan(){

    int idx=1;
    for(int i=1;i<=n;i++){

        int id=idx;
        for(int j=idx;j<=m;j++)if(fabs(a[j][i])>fabs(a[id][i]))id=j;
        if(fabs(a[id][i])>1e-10)for(int k=i;k<=n+1;k++)swap(a[id][k],a[idx][k]);//将当前系数最大的方程交换到第i行
        else continue;//当前主元系数为0,此时不应该加idx !!
        for(int j=n+1;j>=i;j--)a[idx][j]/=a[idx][i];
        for(int j=1;j<=m;j++)if(j!=idx){//将其他方程的这一元消去,最后只留下对角线的元
            for(int k=n+1;k>0;k--)a[j][k]-=a[j][i]*a[idx][k];
        }
        idx++;

    }

    for(int i=idx;i<=m;i++){
        if(fabs(a[i][n+1])>1e-10)return -1;
    }
    for(int i=1;i<=m;i++){
        int flg=0;
        for(int j=1;j<=n;j++)if(fabs(a[i][j])<1e-10&&fabs(a[i][n+1])>1e-10)flg++;//先判无解
        if(flg==n)return -1;
    }
    for(int i=1;i<=n;i++)if(fabs(a[i][i])<1e-10)return 0;//再判多解
    return 1;
}

signed main(){
    scanf("%d%d",&n,&m);
    for(int i=1;i<=m;i++){
        for(int j=1;j<=n+1;j++)cin>>a[i][j];
    }

    int flg=gauss_jordan();
    if(flg==-1)cout<<"No solutions";
    if(flg==0)cout<<"Many solutions";
    if(flg==1){
        for(int i=1;i<=n;i++)cout<<(int)(a[i][n+1]+0.5)<<endl;//printf("%.2Lf\n",a[i][n+1]);
    }


    return 0;
}

关于一组hack

因此，以下写法是错误的！

for(int i=1;i<=m;i++)if(fabs(a[i][i])<1e-10&&fabs(a[i][n+1])>1e-10)return -1;//先判无解

Hack¶

Many solutions

No solutions

无解当且仅当左边方程的所有项系数都为0，且右边常数项不为0

需要判断当前当前消元消去了几个主元，而非第几个元，若只有x个主元，则只固定前x个方程。

例题 #2 [JSOI2008] 球形空间产生器¶

题目描述

有一个球形空间产生器能够在 $n$ 维空间中产生一个坚硬的球体。现在，你被困在了这个 $n$ 维球体中，你只知道球面上 $n+1$ 个点的坐标，你需要以最快的速度确定这个 $n$ 维球体的球心坐标，以便于摧毁这个球形空间产生器。

输入格式

第一行是一个整数 $n$ $(1\le N\le 10)$。接下来的 $n+1$ 行，每行有 $n$ 个实数，表示球面上一点的 $n$ 维坐标。每一个实数精确到小数点后 $6$ 位，且其绝对值都不超过 $20000$。

输出格式

有且只有一行，依次给出球心的 $n$ 维坐标（ $n$ 个实数），两个实数之间用一个空格隔开。每个实数精确到小数点后 $3$ 位。数据保证有解。你的答案必须和标准输出一模一样才能够得分。

提示：给出两个定义：

球心：到球面上任意一点距离都相等的点。
距离：设两个 $n$ 维空间上的点 $A,B$ 的坐标为 $(a_1, a_2, \cdots , a_n), (b_1, b_2, \cdots , b_n)$，则 $A,B$ 的距离定义为：$dist = \sqrt{ (a_1-b_1)^2 + (a_2-b_2)^2 + \cdots + (a_n-b_n)^2 }$。

设球心$O(x_1,\dots,x_n)$，那么我们对于每一个$i\in[1,n+1]$有

$\sum_{j\in[1,n+1]}(a_{i,j}-x_i)^2=C$，其中C是一个定值。

因为我们发现这里是一个关于x的二次式，所以我们考虑相邻两项相减，这样我们就可以得到关于$x_i$的一次方程组了。

/*                                                                                
                      Keyblinds Guide
                    ###################
      @Ntsc 2024

      - Ctrl+Alt+G then P : Enter luogu problem details
      - Ctrl+Alt+B : Run all cases in CPH
      - ctrl+D : choose this and dump to the next
      - ctrl+Shift+L : choose all like this
      - ctrl+K then ctrl+W: close all
      - Alt+la/ra : move mouse to preHash/nxt pos'

*/

//头文件
#include <bits/stdc++.h>

//数据类型
#define ll long long
#define ull unsigned long long
#define db double
#define endl '\n'
//命名空间
using namespace std;
//常量
const int N = 2e5 + 5;
const int M = 1e3;
const int MOD = 903250223;
const int INF = 1e9;



#define rd read()
int read(){
    int xx = 0, ff = 1;
    char ch = getchar();
    while (ch < '0' || ch > '9') {
        if (ch == '-')
            ff = -1;
        ch = getchar();
    }
    while (ch >= '0' && ch <= '9')
      xx = xx * 10 + (ch - '0'), ch = getchar();
    return xx * ff;
}
void write(int out) {
    if (out < 0)
        putchar('-'), out = -out;
    if (out > 9)
        write(out / 10);
    putchar(out % 10 + '0');
}




//变量//变量
long double a[M][M];
int n,b,c,x[N],y[N],ans,res,tmp,cnt,m;

int gauss_jordan(){

    int idx=1;
    for(int i=1;i<=n;i++){

        int id=idx;
        for(int j=idx;j<=m;j++)if(fabs(a[j][i])>fabs(a[id][i]))id=j;
        if(fabs(a[id][i])>1e-10)for(int k=i;k<=n+1;k++)swap(a[id][k],a[idx][k]);//将当前系数最大的方程交换到第i行
        else continue;//当前主元系数为0,此时不应该加idx !!
        for(int j=n+1;j>=i;j--)a[idx][j]/=a[idx][i];
        for(int j=1;j<=m;j++)if(j!=idx){//将其他方程的这一元消去,最后只留下对角线的元
            for(int k=n+1;k>0;k--)a[j][k]-=a[j][i]*a[idx][k];
        }
        idx++;

    }

    for(int i=idx;i<=m;i++){
        if(fabs(a[i][n+1])>1e-10)return -1;
    }
    for(int i=1;i<=m;i++){
        int flg=0;
        for(int j=1;j<=n;j++)if(fabs(a[i][j])<1e-10&&fabs(a[i][n+1])>1e-10)flg++;//先判无解
        if(flg==n)return -1;
    }
    for(int i=1;i<=n;i++)if(fabs(a[i][i])<1e-10)return 0;//再判多解
    return 1;
}

long double ip[M][M];


signed main(){
    n=rd,m=n;
    for(int i=1;i<=n+1;i++){
        for(int j=1;j<=n;j++){
            cin>>ip[i][j];
        }
    }



    for(int i=1;i<=n;i++){

        long double sum=0;
        for(int j=1;j<=n;j++){
            a[i][j]=(ip[i][j]-ip[i+1][j])*2;
            sum+=ip[i][j]*ip[i][j]-ip[i+1][j]*ip[i+1][j];
        }

        a[i][n+1]=sum;
    }

    int flg=gauss_jordan();
    for(int i=1;i<=n;i++)printf("%.3Lf ",a[i][n+1]);



    return 0;
}

例题 #3 [USACO09NOV] Lights G¶

给出一张 $n$ 个点 $m$ 条边的无向图，每个点的初始状态都为 $0$。

你可以操作任意一个点，操作结束后该点以及所有与该点相邻的点的状态都会改变，由 $0$ 变成 $1$ 或由 $1$ 变成 $0$。

你需要求出最少的操作次数，使得在所有操作完成之后所有 $n$ 个点的状态都是 $1$。

输入格式

第一行两个整数 $n, m$。

之后 $m$ 行，每行两个整数 $a, b$，表示在点 $a, b$ 之间有一条边。

输出格式

一行一个整数，表示最少需要的操作次数。

本题保证有解。

对于 $100\%$ 的数据，$1\le n\le35,1\le m\le595, 1\le a,b\le n$。保证没有重边和自环。

本题也可以使用搜索优化算法折半搜索。

这里我们介绍一个01异或方程组。那么本题我们怎么样列出方程组呢？

我们需要用一个式子表示出每个灯的最终状态。我们发现每个灯的最终状态就是这个灯和与之相连的灯的状态的异或值。因此我们可以列出方程组

$\begin{cases}a_1\oplus a_{\dots}=1 \\ a_2\oplus a_{\dots}=1\\ \dots\end{cases}$

欸，这个矩阵好像有一点不一样啊。

对于一个灯我只有按或者不按两种状态，如果按也只会按一次，最终我们希望所有的灯的状态都为1，那么对于每一个灯i的状态就是与它相连的每一个灯包括它自己按或者不按的状态d[i]的异或和(a[1][i]∗d[1])xor(a[2][i]∗d[2])⋯(a[n][i]∗d[n])=1待求的元素就是d[]，就可以得到一个n*(n+1)方程组。

然后用高斯消元之后会得到一个上三角矩阵，但你可能会产生一些自由元，此时我们采用DFS的方法，消完元后的矩阵相当于重新定义了两个灯之间的关系以及每个灯要到达的状态，由于该矩阵是一个上三角矩阵，意味着此时你要做的就是让所有的灯都符合等号右边的状态，由于新的关系可以保证每一个灯都只会被他后边的灯所影响(注意在消元之后两个灯之间的新关系不一定是双向的)，因此我们从后往前搜，如果这个不是自由元那么我们就可以根据后边灯的状态来判断它按或者不按，如是自由元的话，我们就要枚举一下按或者不按，再分别搜索。