高斯消元法

作用
可以快速求解n元线性方程组：

\[\begin{cases} a_{11}x_1+a_{12}x_2+a_{13}x_3+\dots+a_{1n}x_n=b_1\\ a_{21}x_1+a_{22}x_2+a_{23}x_3+\dots+a_{2n}x_n=b_2\\ \dots\\ a_{n1}x_1+a_{n2}x_2+a_{n3}x_3+\dots+a_{nn}x_{n}=b_n\\ \end{cases} \]

思路
利用线性代数知识，存下方程的增广矩阵

\[\left( \begin{array}{ccc} a_{11} & a_{12} & a_{13} & \dots & a_{1n} & b_1\\ a_{21} & a_{22} & a_{23} & \dots & a_{2n} & b_2\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \dots & a_{nn} & b_n \end{array} \right) \]

通过初等行变换将其化为行阶梯型矩阵，即可求解。
行变换步骤：
（1）对每一列进行遍历找到绝对值最大的元素所在的行。^[1]
（2）将这一行的所有元素与当前的第一行所有元素进行交换。
（3）将当前第一行的第一个元素化为1。
（4）从第二行开始把当前第一列的元素全部消为0。

代码实现

#include<iostream>
#include<algorithm>
#include<cmath>
const int N=110;
using namespace std;
const double eps=1e-6;//代表误差
int n;
double a[N][N];
int gauss()
{
    int c,r;//分别表示当前最小列，行
    for(c=0,r=0;c<n;c++)
    {
        int t=r;
        for(int i=r;i<n;i++ )
        {
            if(fabs(a[i][c])>fabs(a[t][c]))//注意浮点数是有误差的，不能直接写等于0。
            {
                t=i;
            }
        }
        if(fabs(a[t][c])<eps)
        {
            continue;
        }
        for(int j=c;j<=n;j++)
        {
            swap(a[t][j],a[r][j]);
        }
        for(int j=n;j>=c;j--)
        {
            a[r][j]/=a[r][c];
        }
        for(int j=r+1;j<n;j++)
        {
            if(fabs(a[j][c])>eps)
            {
                for(int k=n;k>=c;k--)
             {
                a[j][k]-=a[j][c]*a[r][k];
             }
            }
           
        }
        r++;
    }//其实每一次没有被continue的循环就是删去系数矩阵的一条长和一条宽，变为一个较小的正方形。
    if(r<n)//n-r等于几那么最后这几行就全为0
    {
        for(int i=r;i<n;i++)
        {
            if(fabs(a[i][n])>eps)
            {
                return 2;
            }
        }
        return 1;
    }
    for(int i=n-2;i>=0;i--)//从下往上消去每个元素上面的元素
    {
        for(int j=i+1;j<n;j++)
        {
            a[i][n]-=a[j][n]*a[i][j];
        }
    }
     return 0;
}
int main()
{
    cin>>n;
    for(int i=0;i<n;i++)
        {
            for(int j=0;j<=n;j++)
            {
                scanf("%lf",&a[i][j]);
            }
        }
    int res=gauss();
    if(res==0)
    {
        for(int i=0;i<n;i++)
        {
            printf("%.2lf\n",a[i][n]);
        }
    }
    else if(res==2)
    {
        printf("No solution");
    }
    else if(res==1)
    {
        printf("Infinite group solutions");
    }
    return 0;
}

-代码细节
c++里的浮点数储存是有误差的，判断一个数x等于a等价于fabs(x-a)<eps。

要找当前第一列绝对值最大的元素原因是：如果最大的元素的绝对值为1则不用进行下面的步骤直接跳过即可。 ↩︎

526互联

高斯消元法求线性方程组

高斯消元法