嘘~ 正在从服务器偷取页面 . . .

浅谈高斯消元

算法数学

发布日期: 2022-10-05

更新日期: 2022-10-05

文章字数: 680

浅谈高斯消元

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

本文采用高斯-约旦消元法（英语：Gauss-Jordan Elimination）。

相较传统的高斯消元，它不需要回代的过程，所以比较方便。

高斯消元通常用于求解 $n$ 元线性方程组的解（或返回无解），下面是 $n=3$ 的例子

$\begin{cases} 3 x+2 y+z=10 \\[6pt]5 x+y+6z=25 \\[6pt]2 x+3 y+4 z=20 \end{cases}$

我们可以把它写成矩阵的形式

$\begin{bmatrix} &3&2&1&10& \\[6pt]&5&1&6&25& \\[6pt]&2&3&4&20& \end{bmatrix}$

一般地，考虑 $n$ 元线性方程组

$\begin{cases} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 \\[6pt]a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n = b_2 \\[6pt]\quad\vdots \\[6pt]a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n = b_2 \end{cases}$

其矩阵形式为

$\begin{bmatrix} &a_{11}&a_{12}&\cdots&a_{1n}&b_1& \\[6pt]&a_{21}&a_{22}&\cdots&a_{2n}&b_2& \\[6pt]&\vdots&\vdots&\ddots&\vdots&\vdots \\[6pt]&a_{n1}&a_{n2}&\cdots&a_{nn}&b_n& \end{bmatrix}$

依次考虑每个 $a_{ii}$ ，将除第 $i$ 行的每一行第 $i$ 列都变为 $0$ 。

具体操作为，对于每个 $a_{jk}$ ，使其变为 $a_{jk} - a_{ik} \times \frac{a_{ji}}{a_{ii}}$ 。

最后 $x_i$ 的解就是 $\frac{b_i}{a_{ii}}$ ，且最后的矩阵应该长这样（空白处为 $0$ ）

$\begin{bmatrix} &a_{11}&&&&b_1& \\[6pt]&&a_{22}&&&b_2& \\[6pt]&&&\ddots&&\vdots \\[6pt]&&&&a_{nn}&b_n& \end{bmatrix}$

实现时有一些细节，每次把第 $i \sim n$ 行的第 $i$ 列的系数最大值找出来

把它与第 $i$ 行交换，这样可以保证算法的稳定性。（如果最后 $a_{ii}=0$ 则无解）

时间复杂度 $\mathcal{O}(n^3)$

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)(150))

const double eps = 1e-10;
int dcmp(int x)
{
    if(abs(x) <= eps) return 0;
    return x > eps ? 1 : -1;
}
int n; double a[N][N];
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    cout << fixed << setprecision(2);
    cin >> n;
    for(int i=1; i<=n; i++)
        for(int j=1; j<=n+1; j++)
            cin >> a[i][j];
    for(int i=1; i<=n; i++) // col
    {
        int p = i;
        for(int j=i+1; j<=n; j++)
            if(abs(a[j][i]) > abs(a[p][i])) p=j; // row
        for(int j=1; j<=n+1; j++)
            swap(a[i][j],a[p][j]); // swap p to i for proc
        if(!dcmp(a[i][i])) return cout << "No Solution\n",0;
        for(int j=1; j<=n; j++)
            if(j != i)
            {
                double tmp = a[j][i] / a[i][i];
                for(int k=1; k<=n+1; k++)
                    a[j][k] -= a[i][k] * tmp;
            }
    }
    for(int i=1; i<=n; i++)
        cout << a[i][n+1] / a[i][i] << '\n';
    return 0;
}