浅谈高斯消元-CSDN博客

现在有一个三元方程：

 2x+y+z=1;
 6x+2y+z=-1;
 -2x+2y+z=7;

很显然如果你有兴趣解这个方程组的解为

 x=-1.00
 y=2.00
 z=1.00

如果这道题出在数学试卷上相信你一定会做，三元一次方程组嘛！那么这里就写个程序模拟你的运算过程.

为了书写方便，这里化为矩阵的写法：

       x y z val
line1: 2 1 1  1
line2: 6 2 1 -1
line3:-2 2 1  7

用等式1的x消去等式2、3的第一个x

        x  y  z  val
line1:  2  1  1   7
line2:  0 -1 -2  -4
line3: -2  2  1   7(line1*(-3)+line2)
-->
       x  y   z   val
line1: 2  1   1   7
line2: 0 -1  -2  -4
line3: 0  3   2   8 （line1*1+line3）

接下来我们拿line2第二个数消掉下面的一个数

       x  y  z  val
line1: 2  1  1  7
line2: 0 -1 -2 -4
line3: 0  0 -4 -4 （line2*3+line3）

这样我们就得到一个左下角一块全是0的矩阵，

我们将之称之为上三角矩阵就是这个东东

       x  y  z  
line1: 2  1  1  
line2: 0 -1 -2 
line3: 0  0 -4 
(去掉val)

那么我们还是看这个矩阵

       x  y  z  val
line1: 2  1  1  7
line2: 0 -1 -2 -4
line3: 0  0 -4 -4
（加上val）

对于line3我可以轻易的知道z=(-4)/(-4)=1.00
那么知道z的值那么可以推出line2中的y,y=(-4-(-2)*1)/(-1)=2其实就是把z=1代入line2啦
知道y，z的值代入line我们就可以计算出x=-1啦（方法和上面一样）

An another question:

怎么判断是不是无解或是多组解呢？

对于处理处的上三角矩阵，如果全是0（含val）那么多组解

如果有一行全是0（不含val）但是val的值不为0那么无解

这里有个问题必须要说明，对于主元为0的情况因为除数为0，可能存在错误，我们采用交换行的方式避免错误，如

HACK!数据 Ghastlcon 2018-07-17 11:46
HACK.in

3
1 1 1 3
1 1 2 4
3 2 2 7

HCAK.out

1.00
1.00
1.00

这组数据,我们尝试分析：

      x y z val
line1 1 1 1 3 
line2 1 1 2 4 
line3 3 2 2 7 
--->      
      x  y  z val
line1 1  1  1  3 
line2 0  0 -1 -1 (line1-line2) 
line3 0 -1 -1 -2 (line1*(-3)+line3)

我们发现line2和line3位置颠倒，换一下刚好是上三角矩阵，那就换下吧

      x  y  z val
line1 1  1  1  3 
line2 0 -1 -1 -2 
line3 0  0 -1 -1

然后按照上述方案进行处理如果是多元的从上倒下依次交换像插入排序一样，就可以保证最后处理出的矩阵是上三角矩阵。

那么我们可以模拟上述过程，写出一个高斯消元的模板：

# include <bits/stdc++.h>
using namespace std;
const int MAXN=105;
      double eps=0.000001; //精度推荐1e-4到1e-10之间
double a[MAXN][MAXN],val[MAXN];
int n;
bool flag;
void swapline(int i)//交换i行和下面的行，保证i行及以上有序
{
    int j=i+1;
    while (fabs(a[j][i]) < eps && j<=n) ++j;
    if (j>n) {
        printf("No Solution\n");
        flag=true;
        return;
    }
    if (j<=n) {swap(a[i],a[j]);swap(val[i],val[j]);}
}
int main()
{
    scanf("%d",&n);
    for (int i=1;i<=n;i++)  {
        for (int j=1;j<=n;j++) scanf("%lf",&a[i][j]);
        scanf("%lf",&val[i]);
    }
    for (int i=1;i<=n-1;i++) {
        double num=a[i][i];
        if (fabs(num)<eps)  swapline(i);
        if (flag) return 0;
        for (int j=i+1;j<=n;j++) {
            if (fabs(a[j][i])<=eps) continue;//当此行主元a[i][i]=0时交换行
            double w=-a[j][i]/num; a[j][i]=0;
            for (int k=i+1;k<=n;k++) a[j][k]+=a[i][k]*w;
            val[j]+=val[i]*w;
        }
    }//求上三角矩阵
    for (int i=1;i<=n;i++) {
         flag=true;
        for (int j=1;j<=n;j++)
         if (fabs(a[i][j])>eps) {
            flag=false;break;
         }
        if (flag) {
            printf("No Solution\n");
            return 0; //判断是否无解或有多组解
        }
    }
    val[n]=val[n]/a[n][n];//求出第n个元
    for (int i=n-1;i>=1;i--) {
        double sum=0;
        for (int j=n;j>=i+1;j--) sum=sum+a[i][j]*val[j];
        val[i]=(val[i]-sum)/a[i][i];
    }//代入已知元求出未知的一个元
    for (int i=1;i<=n;i++) printf("%.2lf\n",val[i]);
    return 0;
}