高斯消元

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

  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\n", a[i][n]);
    }
    else if (t == 1) cout << "Infinite group solutions" << endl;
    else cout << "No solution" << endl;
   
    return 0;
}
posted @ 2021-02-05 23:49  0x7F  阅读(63)  评论(0编辑  收藏  举报