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

    这里暂时并不会涉及到复杂的线性代数中的情况
    化为上三角行列式之,有三种情况

    1. n个变元,n个方程:唯一解
    2. n个变元,小于n个方程,且其余方程不存在矛盾:无穷解
    3. 存在矛盾:无解

    算法流程

    1.从第一列开始,找到每一列绝对值最大的那一行
    2.把上一步找到的那一行放到最上面
    3.把最上面一行该列的系数变为1
    4.把该列其他行系数变为0
    5.如果有唯一解,从最后一个式子开始向上逐渐消元直至得到最简上三角行列式

    代码实现

    #include <iostream>
    #include <cmath>
    #include <algorithm>
    #include <stdio.h>
    
    using namespace std;
    
    const int N = 110;
    const double eps = 1e-6; // 这个数字的定法是个比较玄学的问题,因为不好说多小才算是0,0.01就不应该算作0,所以1e-6应该说是经验值
    
    int n;
    double a[N][N];
    
    int guass()
    {
        // 完成高斯消元的模拟过程并返回解的情况
        int c, r;
        for (c = 0, r = 0; c < n; ++ c) // 列,每次循环确定一个未知量的解
        {
            // 第一步:找到c列绝对值最大的行数
            int t = r; // c列系数最大的一行
            for (int i = r; i < n; ++ i) // 行
                if (fabs(a[i][c]) > fabs(a[t][c]))
                    t = i;
           
            /**
             * 实际目的为判断是否为0,但是由于精度问题,可能出现0.00000001这种类似的情况,所以采用这种方式,在浮点数二分时也采用过这种方式
             * 这里的contineu不太好理解
             * 我们需要从r的意义入手,我们希望每次r都能确定一个解,如果当前列的系数均为0,我们是没办法确定一个解的,所以去看下一列
             * 所以最终如果是唯一解,那么r应该是等于方程个数的,因为每个方程确定一个解,如果不等于,就要看方程右边是不是存在某个常数不为0,如果存在,说明存在0 = 非0这种式子,所以是无解,否则是无穷多解
             */
            if (fabs(a[t][c]) < eps) continue;
           
            // 第二步:将t行换到最上方还未确定的一行
            for (int i = 0; i < n + 1; ++ i) swap(a[r][i], a[t][i]);
           
            // 第三步:将首位系数变为1
            for (int i = n; i >= 0; -- i) a[r][i] /= a[r][c];
           
            // 第四步:将下方c列系数变为0
            for (int i = r + 1; i < n; ++ i)
                for (int j = n; j >= 0; -- j) // 需要从后向前更新,以保留该行第一个系数
                    a[i][j] -= a[i][c] * a[r][j];
            ++ r;
        }
        
        /**
        * 如果是唯一解,那么每一行确定一个一个解
        * 假设有3个方程,3个未知量,那么r=0时处理第一个未知量,=1时第二个,=2时第三个,=3时结束了,也就是r最终停止的位置是不确定解的位置
        * 所以如果是无解或者无穷解,r-1是最后一个确定方程解得位置,r是第一个全零行,所以判断是无解还是无穷解是从r开始的,而非r+1
        */
        // 第五步:判断是否有解并处理可消去可消系数
        if (r < n)
        {
            for (int i = r; i < n; ++ i)
                if (fabs(a[r][n]) > eps)
                    return 2;
            return 1;
        }
        
        // 用i行下面每一行来把i行对应的系数消为0,常数列对应改变
        // 这里i必须从下向上走,因为消i行的某个系数时,必须保证j行其它系数为0了才能保证得到正确的结果,所以要先消下面的系数
        for (int i = n - 1; 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 + 1; ++ j)
                cin >> a[i][j];
       
        int t = guass();
       
        if (!t)
        {
            // 此时有唯一解,最终方程的目标格式为每个式子只有一个x且该x前系数为1,所以最后的常数就是解
            for (int i = 0; i < n; ++ i)
                printf("%.2lf
    ", a[i][n]);
        }
        else if (t == 1) cout << "Infinite group solutions" << endl;
        else cout << "No solution" << endl;
       
        return 0;
    }
    
  • 相关阅读:
    贪心法之活动安排问题
    动态规划算法之最优二叉搜索树
    动态规划之最大字段和问题
    动态规划算法之图像压缩问题
    动态规划算法之0-1背包问题
    动态规划算法之投资问题
    平面点集的凸包问题
    动态规划(DP)之多边形游戏问题
    凸多边形最优三角划分
    最长公共子序列问题
  • 原文地址:https://www.cnblogs.com/G-H-Y/p/14380286.html
Copyright © 2011-2022 走看看