acwing算法基础之数学知识--高斯消元法求解线性方程组

1 基础知识

高斯消元法，用来求解线性方程组的解，
{ a 11 x 1 + a 12 x 2 + ⋯ + a 1 n x n = b 1 a 21 x 1 + a 22 x 2 + ⋯ + a 2 n x n = b 2 ⋯ a n 1 x 1 + a n 2 x 2 + ⋯ + a n n x n = b n \left \{

\right. ⎩ ⎨ ⎧​a11​x1​+a12​x2​+⋯+a1n​xn​=b1​a21​x1​+a22​x2​+⋯+a2n​xn​=b2​⋯an1​x1​+an2​x2​+⋯+ann​xn​=bn​​

初等行列变换：

1. 某行乘以k倍，即row[i] = k * row[i]
2. 交换某两行，即row[i], row[j] = row[j], row[i]
3. 某行加上另一行的k倍，即row[i] = row[i] + k * row[j]

高斯消元法代码实现的关键步骤：

• 枚举每一列c，进行以下操作：
  – 找到绝对值最大的一行（注意，要从不固定的行中寻找）。
  – 将该行换到最上面。
  – 将该行的第1个数变成1。
  – 将下面所有行的第c列消成0。
• 经过上述步骤，行列式已转变为下三角行列式，需要从第n行（或者第n-1行）开始把未知数$x_i$的值代入方程进行求解，依次求解出$x_n$、$x_{n-1}$、……、$x_1$的值。

将上述操作转换为代码，如下，

#include 
#include 

using namespace std;

const int N = 110;
const double eps = 1e-6;

double a[N][N];
int n;

int gauss() {
    //从0到n-1遍历每一列，把行列式消成下三角形式
    int r, c;
    for (r = 0, c = 0; c < n; ++c) {
        //正在处理第r行第c列，a[r][c]
        //step1：找到r~n-1行中第c列绝对值最大的行，记为t
        int t = r;
        for (int i = r + 1; i < n; ++i) {
            if (fabs(a[i][c]) > fabs(a[t][c])) {
                t = i;
            }
        }
        //step2：交换第r行和第t行的第c~n列的元素值
        for (int j = c; j <= n; ++j) {
            swap(a[r][j], a[t][j]);
        }
        //step3：将第r行的变量x_c前面的系数变为1，即将a[r][c]变为1
        if (fabs(a[r][c]) < eps) {
            continue;
        }
        for (int j = n; j >= c; --j) {
            a[r][j] = a[r][j] / a[r][c];
        }
        //step4：将第r+1~n-1行中的第c列给去掉，即将第r行的某倍数加到下面的行中
        for (int i = r + 1; i < n; ++i) {
            //将第i行中的第c列给去掉（即x_c的系数为0）
            for (int j = n; j >= c; --j) {
                a[i][j] -= a[i][c] * a[r][j]; //将第r行的-a[i][c]倍加到第i行
            }
        }
        
        r += 1; //有效方程数加1
    }
    
    if (r < n) {
        for (int i = r; i < n; ++i) {
            if (fabs(a[i][n]) > eps) {
                return 2; //0等于非0，无解。
            }
        }
        return 1; //有效方程数小于未知量的个数，无穷多组解。
    }
    
    //step5：逆向求解x_{n-1},x_{n-2}...x_0
    for (int i = n - 1; i >= 0; --i) {
        //求解未知数x_i，第i行第i列
        for (int j = i + 1; j < n; ++j) {
            a[i][n] -= a[i][j] * a[j][n]; //a[j][n]表示x_j的值
        }
    }
    return 0;
}

int main() {
    cin >> n;
    //总共n行、n+1列
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n + 1; ++j) {
            cin >> a[i][j];
        }
    }
    
    int ans = gauss();
    
    if (ans == 0) {
        for (int i = 0; i < n; ++i) {
            printf("%.2lf\n", a[i][n]); //xi的解为a[i][n]
        }
    } else if (ans == 1) {
        cout << "Infinite group solutions" << endl;
    } else {
        cout << "No solution" << endl;
    }
    
    return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

2 模板

// a[N][N]是增广矩阵
int gauss()
{
    int c, r;
    for (c = 0, r = 0; c < n; c ++ )
    {
        int t = r;
        for (int i = r; i < n; i ++ )   // 找到绝对值最大的行
            if (fabs(a[i][c]) > fabs(a[t][c]))
                t = i;

        if (fabs(a[t][c]) < eps) continue;

        for (int i = c; i <= n; i ++ ) swap(a[t][i], a[r][i]);      // 将绝对值最大的行换到最顶端
        for (int i = n; i >= c; i -- ) a[r][i] /= a[r][c];      // 将当前行的首位变成1
        for (int i = r + 1; i < n; i ++ )       // 用当前行将下面所有的列消成0
            if (fabs(a[i][c]) > eps)
                for (int j = n; j >= c; j -- )
                    a[i][j] -= a[r][j] * a[i][c];

        r ++ ;
    }

    if (r < n)
    {
        for (int i = r; i < n; i ++ )
            if (fabs(a[i][n]) > eps)
                return 2; // 无解
        return 1; // 有无穷多组解
    }

    for (int i = n - 1; i >= 0; i -- )
        for (int j = i + 1; j < n; j ++ )
            a[i][n] -= a[i][j] * a[j][n];

    return 0; // 有唯一解
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

3 工程化

暂无。。。