python高斯数组去噪

• -知识点整理-高斯消元-

• 典型题目
• 知识点
• 代码实现

->知识点整理-高斯消元<-

典型题目：

XJOI 1822：Civilization

知识点：

高斯消元其实在小学初中解多元一次方程的时候已经接触过了。其实，高斯消元就是建立在方程中加减消元和乘除消元之上的。只不过，高斯消元法把这两种方法应用于矩阵之中，使得高斯消元的复杂度达到O（n³）（相比于真正的去解方程可是要快的多了，想一想你手解100000元一次方程需要多久就行了）。
然后就是高斯消元的实现了。首先，我们要把高斯消元的每一个系数都加入到N*（N+1）的矩阵中（第N行第N+1列填的是这个方程的答案，但在最后要被我们转化为第N元的解）。

举个例子：假设我们要解如下一个三元一次方程组。

首先，把他们的系数都加入一个矩阵当中：

然后，我们要确定我们所期望的解决后的矩阵是怎样的。这样的矩阵无非两种，一种是上三角矩阵，第i行的元素个数为n-i+1个，并且第i行的前i-1个元素都必须为0。然后从最后一行开始解方程，不断地代入上一行的方程，知道解除答案。而另一种方法则是第i行只有第i个元素不为0，这样可以直接解出每个数的解。然而前者实现较为简单，所以我们今天就用前者。

然后开始我们的消元。首先先从第一列开始，先把该列中绝对值最大的元素所在的哪一行提至第一列（目的是为了提高数据的稳定性，数学运算都是含舍入误差的计算，选主元进行消去可以极大降低舍入误差）。然后就想简单的乘除加减消元一样，把下面的每一列的值都消为0。

下面是第一步举例：

最后是结果的矩阵：

现在只需要从下往上一个一个的解出每一元的解即可。最后得出的解为x=2，y=5，z=8。
然后该怎么判断方程是否有解呢？
只需要按下面的方法判断即可：
① 无解
当方程中出现(0, 0, …, 0, a)的形式，且a != 0时，说明是无解的。
② 唯一解
条件是k = equ，即行阶梯阵形成了严格的上三角阵。利用回代逐一求出解集。
③ 无穷解。
条件是k < equ，即不能形成严格的上三角形，自由变元的个数即为equ – k。

代码实现：

#include<bits/stdc++.h>
using namespace std;
typedef long double ld;//这里用了long double，是为了精度，但是XJOI原题需要用另一种方法来记录消元时的数，有兴趣可自己去试试
const int maxn=205;

ld a[maxn][maxn];
char in[maxn];//由于XJOI好像不能读long double，所以用了读字符串的方式

void guass(ld x[maxn][maxn],int n)
{
    int r;
    for(int i=0;i<n;i++)
    {
        r=i;
        for(int j=i+1;j<n;j++)
        {
            if(fabs(x[j][i])>fabs(x[r][i]))r=j;
        }
        if(r!=i)for(int j=0;j<=n;j++)swap(x[r][j],x[i][j]);//找出绝对值最大的那一列并进行交换

        for(int j=n;j>=i;j--)
        {
            for(int k=i+1;k<n;k++)
            {
                x[k][j]-=x[k][i]/x[i][i]*x[i][j];//如果中间用一个中间值进行计算，有可能会造成精度误差，所以我直接进行了计算
            }
        }
    }

    for(int i=n-1;i>=0;i--)
    {
        for(int j=i+1;j<n;j++)
        {
            x[i][n]-=x[j][n]*x[i][j];
        }
        x[i][n]/=x[i][i];
    }
}

int main()
{
    int n;
    scanf("%d",&n);
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<=n;j++)
        {
            scanf("%s",in);
            int r=0;
            for(int k=strlen(in)-1;k>=0;k--)
            {
                a[i][j]+=(ld)(in[k]-'0')*pow(10,r);
                r++;
            }
        }
    }

    guass(a,n);
    double res;
    for(int i=0;i<n;i++)
    {
        res=(double)a[i][n];//硬生生转回了double输出
        printf("%.0lf\n",res);
    }
}