zoukankan      html  css  js  c++  java
  • 浅谈高斯消元

    题目链接: 洛谷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

    用等式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
    ");
            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
    ");
                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
    ",val[i]);
        return 0;
    }
  • 相关阅读:
    客户端命令行查看dubbo服务的生产者和消费者
    灰度发布
    灰度发布方法了解
    spring事务隔离级别
    小东西
    从excel表中生成批量SQL,将数据录入到数据库中
    CSV导出
    详解Java的自动装箱与拆箱(Autoboxing and unboxing)
    java多线程理解
    spring事务传播行为讲解转载
  • 原文地址:https://www.cnblogs.com/ljc20020730/p/9454997.html
Copyright © 2011-2022 走看看