线性方程组-高斯消元

苜小夕

已于 2024-08-27 18:56:04 修改

阅读量1k

点赞数 31

分类专栏： ACM 文章标签：线性代数算法

于 2024-08-27 18:53:29 首次发布

本文链接：https://blog.csdn.net/2201_75771132/article/details/141610112

版权

ACM 专栏收录该内容

10 篇文章 0 订阅

订阅专栏

线性方程组-高斯消元

P2455 [SDOI2006] 线性方程组

已知 $n$ 元线性一次方程组。

$\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases}$

请根据输入的数据，编程输出方程组的解的情况。

高斯消元是什么呢

枚举每一个未知数也就是n列

利用第一个系数将其他的系数全部消为0 (只剩最后一个肯定是消不掉的)

$\begin{cases} ax_1+\cdots=B_1 \\ bx_2+\cdots=B_2 \end{cases}$

利用a消掉b

即第二行减去 $b / a$ 倍的第一行 --> $b - b / a * a = 0$

继续对第三行…第n行做同样的操作这样就可以消掉其他行的这个未知数的系数

容易发现需要用double 可能存在丢精度的问题

只需要使 $b / a$ 变小就可以让丢失的精度尽可能低即找出最大的 $a$

只需要做一遍for循环然后交换这两行即可

这样一定可以是最初矩阵转化为行阶梯型矩阵

最后只需要倒着遍历求解然后回代即可

需要注意的是判断为0不能 $== 0$ 因为小数会有精度误差的问题

判断解的个数

既然最后一定是行阶梯型矩阵

只需要在for循环的时候求矩阵的秩

我们用line表示当前访问到了哪个非0行即目前矩阵的秩+1

即 $a [l in e] [i]$ 如果 $< e p s$ 此时矩阵的秩应该减去1 即line此时不应该+1

无穷多解即零行是这样的：

$\ 0 \ 0 \ 0 \cdots 0$

无解即零行是这样的

$\ 0 \ 0 \ 0 \cdots b \ (b\neq0)$

只需要从line遍历到n 判断一下即可

详见代码

#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<iomanip>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<stack>
//#include<random>
//#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define pii pair<int,int>
#define double long double
const double eps=1e-9;
void solve(){
    int n;
    cin>>n;
    vector<vector<double>>a(n+10,vector<double>(n+10));
    for(int i=1;i<=n;++i){
        for(int j=1;j<=n+1;++j){
            int x;cin>>x;
            a[i][j]=x;
        }
    }
    int line=1;
    for(int i=1;i<=n;i++){
        int maxn=line;
        for(int j=line+1;j<=n;j++){
            if(fabs(a[j][i])>fabs(a[maxn][i]))maxn=j;
        }
        swap(a[maxn],a[line]);
        if(fabs(a[line][i])<eps)continue;
        for(int j=line+1;j<=n;j++){
            for(int k=n+1;k>=i;k--){
                a[j][k]-=a[j][i]/a[line][i]*a[line][k];
            }
        }
        line++;
    }
    for(int i=line;i<=n;i++){
        if(fabs(a[i][n+1])>eps){
            cout<<-1;//无解
            return ;
        }
    }
    if(line<=n){cout<<0;return ;}//无穷多解
    for(int i=n;i>=1;i--){
        for(int j=i+1;j<=n;j++){
            a[i][n+1]-=a[i][j]*a[j][n+1];//回代
        }
        a[i][n+1]/=a[i][i];
    }
    for(int i=1;i<=n;i++)cout<<"x"<<i<<"="<<fixed<<setprecision(5)<<a[i][n+1]<<'\n';
}
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0),cout.tie(0);
    int t=1;
    // cin>>t;
    while(t--)solve();
    return 0;
}