蒟蒻林荫小复习——高斯消元

在无数次逃避之后，林荫终于鼓起了勇气，迈出了向数论进击的第一步！

高斯消元：众所周知就是高斯解方程用的消元方法。

至于啥是消元法：将方程组中的一方程的未知数用含有另一未知数的代数式表示，并将其代入到另一方程中，这就消去了一未知数，得到一解；或将方程组中的一方程倍乘某个常数加到另外一方程中去，也可达到消去一未知数的目的。消元法主要用于二元一次方程组的求解。

听着是不是很强？

然而这就是我们在初一学的解二元一次方程的扩展：尝试将一个未知量用其他的未知量表示，直到最后得到一元一次方程，解开后再反向带入回去即可。

前置芝士：无！

实际上这个是需要懂得一些矩阵的知识的，但是我们可以用普通的方式来解释它（NOIP范围内）

题目：给定N个N元一次多项式，求出解

先找个样例：

3
1 3 4 5
1 4 7 3
9 3 2 2

每一行前3个数分别是3个元（x,y,z）的系数，最后一个数是等式的值

下面我们可以列出一个方阵（不是矩阵）：

1    3   4    5

1    4   7    3

9    3    2   2

观察如下方阵，我们可知：

1. 上下交换任意两整个横行对方阵没有影响（相当于把方程式给出的顺序变一下）
2. 整个横行中所有数同时乘k（实数）对方阵没有影响（相当于把方程式左右同时乘k）

那么由于我们的目的是将前面n*n的方阵消成只有从左上到右下的对角线为1，其他位置为0，

那么这个时候第i行的第n+1列就是第i个未知数的值

现在我们来确定一下消元的步骤:

1. 现在是确定第i个元,找到一个之前没用过的方程记为第i个等式（你总不能通过一个等式同时将两个元用其他元表示吧，如果这样你一定会得到一个恒等式。这也就是为啥要解n元一次方程一定要有n个不同的n元一次方程），并将其用之前提到的性质1转换到第i横行的位置（方便下面继续寻找）
2. 将这个等式同时除以第i个元的系数，使得可以用其他元表示第i个元（假设2x+3y+4z=10，当前元是y，那么等式会变成2/3x+y+4/3z=10/3，实际上这种情况是不可能存在的，因为在消去y之前，x一定已经被消去，x的系数会变为0，这里举出只是为了说明怎么除）
3. 然后开始用这个等式对其他等式进行操作（实际上相当于将第i+1个等式的第i+1个元用其他元表示）
4. 当对每一个元进行过操作之后，矩阵的第n+1列就是每个元的值。

好了，现在小朋友们可能有一个问题，这样算来算去元消在哪里了呢？

我们现在开始考虑这个问题：

当我们对第一个元进行过操作之后，除了第一个式子以外，其他式子中的第一个元都已经由其他元的组合代替。

第二个元的操作也是一样，只有第二个式子中第二个元的系数为1，其他式子（包括第1个式子）的第二个元已经全部有3，4，5，6......等元替代

因为在操作之前第二个式子中的第一个元的系数已经为0(被第1个式子消去了）.这样的话，第i个式子在算法最后只有第i个元的系数为1，其余都是0，这样就达到了我们的目的。

1 3 4 5
0 1 3 -2
0 -24 -34 -43

1 0 -5 11
0 1 3 -2
0 0 38 -91

1 0 0 -0.973684
0 1 0 5.18421
0 0 1 -2.39474

把样例运行之后，输出对每个元操作结束后的方阵就得到了上面的东西

可以看到，算法每次将一个纵列上除了当前第i个元所在位置以外的常数变成0

这样的话，我们最后就可以求出解。

但是，既然是个多元方程，就有可能出现无解和有无穷多个解的情况！

那么怎么判断无解呢？

当算法进行到某个元i时，无法找到一个没有用过且该元系数不为0的方程

换句话说，如果给了n个n元一次方程，其中有一个系数为0，就相当于给了一个n-1元一次方程，可以直接判定无解。

先来份代码！

#include<iostream>
#include<cstdio>
using namespace std;
int n;
double wws[110][110];
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n+1;j++)
        {
            cin>>wws[i][j];
        }
    }
    for(int i=1;i<=n;i++)//枚举每个元
    {
        int pl=i;
        while(pl<=n&&wws[pl][i]==0)
        {
            pl++;
        }
        if(pl>n)
        {
            cout<<"No Solution"<<endl;
            return 0;
        }
        for(int j=1;j<=n+1;j++)//把第一个合法的行翻到这个元的位置 
        {
            swap(wws[pl][j],wws[i][j]);
        }
        double k=wws[i][i];
        for(int j=1;j<=n+1;j++)
        {
            wws[i][j]=wws[i][j]/k;//对自己行的处理，确定当前目标元的常数=1 
        }
        for(int m=1;m<=n;m++)
        {
            if(m==i)
                continue;
                double ki=wws[m][i];
            for(int k=1;k<=n+1;k++)
            {
                wws[m][k]-=ki*wws[i][k];//因为ws[i]中存有当前元常数为1时每个元对应的常数
                //因此消去时枚举每个方程，ki为该方程中当前元常数，对其他元依次消去 
            } 
        }
    }
    for(int i=1;i<=n;i++)
    {
        printf("%.2lf
",wws[i][n+1]);
    }
    return 0;
}

无解判完了那怎样判断无穷解？

这个先等我研究研究，马上回来填坑！

 再填个小坑

wws[m][k]-=ki*wws[i][k];

代码中这一行可能有点难度哈

假定我们现在有个方程：x+2y+3z=10

那么x用y和z表示就是x=10-2y-3z

那么我们现在尝试将上述x带入2x+3y+z=15中

可得20-4y-6z+3y+z=15

转化一下（3-4）y+（1-6）z=15-20

然后和上面代码结合一下：

 for(int m=1;m<=n;m++)
 {//枚举每个横行
            if(m==i)//是第i行就不管(前面已经处理过）
                continue;
                double ki=wws[m][i];//ki是在第m行中元i的系数（相当于上面2x+3y+z=15中的2，此时元i为x）
            for(int k=1;k<=n+1;k++)
            {//wws[i][k]可以认为是在用其他元表示元i的表达式中元k的个数的相反数，也是式子i中元k的系数（假设当前元k为上面的y，那么wws[i][k]==2）
                wws[m][k]-=ki*wws[i][k];//因为前面存的是相反数，所以要用减法。
                //因此消去时枚举每个方程，ki为该方程中当前元常数，对其他元依次消去
            }
}

判断误解实际上可能还有一种情况,就是式子足够n个,但是有些式子本质是一样的

那么我们消元消完后,可以发现,如果有某一行系数全部是0了,但是常数并不是0,这个时候就代表无解

至于多解,在确定有解后,找一下看有多少个式子系数和常数全是0

这样就代表这个元可以随便取值,怎么取都可以.

那自然就多解了.

UPD: 之前的消元方法不能很好的判断是否有多解或无解

下面给出新的代码

#include<iostream>
#include<cstdio>
using namespace std;
double A[55][55],EPS=1e-8;
int n,t;
double ABS(double x)
{
    return max(x,-x);
}
void Work()
{
    //freopen("gaess.in","r",stdin);
    //freopen("gaess.out","w",stdout);
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&A[i][j]);
    }
    int row,i;
    for(i=1,row=1;row<=n&&i<=n;i++,row++)
    {
        //row代表行,i代表决策元 
        int t=row;
        for(int j=row+1;j<=n;j++)
        {
            if(ABS(A[j][i])>ABS(A[t][i]))
                t=j;
        } 
        if(t!=row)
        for(int j=1;j<=n+1;j++)
        {
            swap(A[t][j],A[row][j]);
        }
        if(ABS(A[row][i])<EPS)
        {
            row--;continue;
        }
        //代表这个元没能找到合适求解的式子,先放一放,式子下面还能用 
        for(int j=row+1;j<=n;j++)
        {
            if(ABS(A[j][i])>EPS)
            {
                double Kel=A[j][i]/A[row][i];
                for(int k=1;k<=n+1;k++)
                {
                    A[j][k]-=A[row][k]*Kel;
                }
            }
        }
    }
    row--;
    for(int j=row;j<=n;j++)
    {
        if(ABS(A[j][n])<EPS&&ABS(A[j][n+1]>EPS))
        {
            printf("%d
",-1);
            return;
        }
    }
    if(row<n)
    {
        printf("%d
",0);
        return;
    }
    for(int j=row;j>=1;j--)
    {
        for(int k=j+1;k<=n;k++)
        {
            A[j][n+1]-=A[j][k]*A[k][n+1];
        }
        A[j][n+1]/=A[j][j];
    }
    for(int i=1;i<=n;i++)
    {
        if(A[i][n+1])
        {    
            cout<<"x"<<i<<"=";
            printf("%.2lf
",A[i][n+1]);
        }
        else
        {
            cout<<"x"<<i<<"="<<0<<endl;
        }
    }
    return;
} 
int main()
{
    Work();
    return 0;
}

多解输出0,无解输出-1

完结撒花！