题解 P3389 【【模板】高斯消元法】

· · 题解

题目链接： 洛谷P3389高斯消元【模板】

提醒：在cnblog博客或是luogu博客看此文更优哦，戳这里

现在有一个三元方程：

 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

step1: 用等式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）

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;
}

截止2018.8.11此程序通过讨论区内的所有HACK数据