高斯消元学习体会

· · 个人记录

高斯消元就像初中学的多元一次方程组（且当n元时有n个方程）的一种解法。

一、解法

其中经过多次(n^3)消元后可呈倒三角形式

\begin{bmatrix} 1 & -1 & 3\\ 0 & 8 & 5\\ 0 & 0 & 4 \end{bmatrix}$ $\begin{bmatrix} 2\\ 3\\ 8 \end{bmatrix}

这样将4x_3=8\Rightarrow x_3=2带入②和①得

\begin{bmatrix} 1 & -1 & 0\\ 0 & 8 & 0\\ 0 & 0 & 0 \end{bmatrix}$ $\begin{bmatrix} -4\\ -7\\ 0\end{bmatrix}

可知8x_2=-7\Rightarrow x_2=-\frac{8}{7}

再次迭代可得

\begin{bmatrix} 1 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{bmatrix}$ $\begin{bmatrix} -\frac{36}{7}\\ 0\\ 0\end{bmatrix} 为了解决问题方便，我们可以将矩阵坐标为$i,j$的位置化为1，即 $\begin{bmatrix} 1 & -1 & 3\\ 0 & 8 & 5\\ 0 & 0 & 4 \end{bmatrix}$ $\begin{bmatrix} 2\\ 3\\ 8 \end{bmatrix}\rightarrow\begin{bmatrix} 1 & -1 & 3\\ 0 & 1 & \frac{5}{8}\\ 0 & 0 & 1 \end{bmatrix}$ $\begin{bmatrix} 2\\ \frac{3}{8}\\ 2 \end{bmatrix}

就可以看得更清楚，让程序直接从i+1做计算

二、精度

为了保证精度，我们可以让小的去扩大，而不是去对大的除，同时要保证被除数不能是0，防止RE（这里是因为当小数点后8位均为0时已无意义，甚至会爆double）

三、判断

判断是否有唯一解 如果做到某一列发现这一列的系数全部为0，说明所在这一列的x为“自由元”，则不为唯一解，这个x可以取任意值。

如果做完后发现系数有矛盾即消不掉的系数，则该矩阵无解。

四、其他

如上这种做法，因为计算机不能存储分数（可以手动实现）精度丢失过大，但在一般的题中已经够用。（见代码高斯消元v1.0）

为防止变态题精度误差太大，做除法之前可以求出这一列的最小公倍数，互相消元时就不会出现小数了（还是要避免0的情况出现），但是这样就不能保证i,j位置上的数是1了，迭代回溯时要注意这一点。（代码待补充）

高斯消元v1.0

#include<cstdio>
#include<algorithm>
double abs(double x)
{
    if(x<0)
        return -x;
    return x;
}
double a[111][111];
int n;
bool gs()
{
    for(int i=1;i<=n;i++)
    {
        double maxx=a[i][i];
        int tmp=i;
        for(int j=i;j<=n;j++)
            if(abs(a[j][i])>abs(maxx))
            {
                maxx=a[j][i];
                tmp=j;
            }
        if(abs(maxx)<=0.0000001)
            return false;
        if(i<n)
            std::swap(a[i],a[tmp]);
        double upd=a[i][i];
        for(int j=i;j<=n+1;j++)
            a[i][j]/=upd;
        for(int k=i+1;k<=n;k++)
        {
            upd=a[k][i];
            for(int j=i;j<=n+1;j++)
                a[k][j]-=upd*a[i][j];
        }
    }
    for(int i=n;i>=1;i--)
    {
        double upd=a[i][n+1];
        for(int j=i-1;j>=1;j--)
        {
            a[j][n+1]-=a[j][i]*upd;
            a[j][i]=0;
        }
    }
    return true;
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&a[i][j]);
    if(gs())
    {
        for(int i=1;i<=n;i++)
            printf("%.2lf\n",a[i][n+1]);
    }
    else
        puts("No Solution");
    return 0;
}