浅谈高斯消元的实现和简单应用

标签：port   高斯消元   desc   模板   一个   rate   imp   script   AC   

一、高斯消元的原理

对于n元的m个线性方程组成的方程组，我们将其以矩阵的形式记录下来：

a11 a12 a13 ...... a1n b1
a21 a22 a23 ...... a2n b2
...
...
...
an1 an2 an3 ...... ann bn

然后进行初等行列变换，尝试构造出一个上三角矩阵，逐步使系数不为零的项减少；

等最后只剩下一个系数不为零时，进行回代，逐步求出已知解。（详解过程咨询小学老师）

二、高斯消元的实现

老老实实的回代代码参见其他人的博客，这里介绍一种比较毒瘤的不回代暴力消元法：

1.Process

对于每个方程，按照一定的规则（后话）挑选一个主元，记录该主元对应第几个方程，然后用初等行列变换消去其他所有该元的系数；

最后我们得到的是一个每行只有一个数不为零，每列只有一个数不为零的鬼畜矩阵（自己脑补）

此时令ans向量对应的数字出去该行的非零系数，即为对应该元的解。

2.Code

设a数组为已经记录系数的数组（格式见上方），那么a应该是n行n+1列的，最后一列代表方程的常数项；

设w数组记录每个方程的主元是第几项，v数组记录答案；

当多解时输出“Multiple solutions”，无解时输出”No solution”；

double a[max_n][max_n+1],v[max_n];

int w[max_n]; 

void gauss(){
    double eps=1e-6;
    for(int i=1;i<=n;++i){    //Enumerate the equation;
        int p=0;                //Record the position of the largest number; 
        double mx=0;        //Recording the largest number;
        for(int j=1;j<=n;++j)
            if(fabs(a[i][j])-eps>mx){
                mx=fabs(a[i][j]);p=j;    //fabs() returns the absolute value of float; 
            }
        if(!p){
            if(fabs(a[i][n+1]<eps))printf("Multiple solutions");
            else printf("No solution");
            return; 
        }
        w[i]=p;
        for(int j=1;j<=n;++j)
            if(i!=j){       //other equations
                double t=a[j][p]/a[i][p];
                for(int k=1;k<=n+1;++k)    //n+1 is important
                    a[j][k]-=a[i][k]*t;
            }
    }
    for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]];
} 

3.notice

（1）精度的设置

众所周知浮点数是有精度丢失的，在高斯消元中，精度的选择要依题目而定，精度过低会导致系数较小的数被误认为系数为零，而精度过高也有可能会导致误差大于精度，导致本应该系数消为0的项误认为系数不为零，所以精度的选择是很哲学的问题，要注意。

推荐范围：1e-4到1e-10

（2）主元的选取原则

接着（1）说，我们知道，用浮点数是有精度丢失的，那么让一个较大的数除以一个极小的数误差自然大的可想而知，所以我们想得到在精度允许的条件下系数最大的主元，所以对于每个方程，我们都应该选择最大系数的元作为主元。

（3）在模2意义下的高斯消元

使用bitset优化运行时间，详见相关应用中第三个例题的讲解；

三、相关应用

这里给出高斯消元的几道基础题目，难度适合初学者。

1.[Luogu P3389]【模板】高斯消元

Description

给定一个线性方程组，对其求解

输入格式：
第一行，一个正整数 n
第二至 n+1行，每行 n+1个整数，为 a1,a2?an和 b，代表一组方程。

输出格式：
共n行，每行一个数，第 i行为 xi(保留2位小数）
如果不存在唯一解，在第一行输出"No Solution".

Solution

如上所述。

Code

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;

const int max_n=110;
double a[max_n][max_n+1],v[max_n];
int n,w[max_n]; 

inline int rd(){
    int x=0;
    bool f=0;
    char c=getchar();
    while(!isdigit(c)){
        if(c=='-') f=1;
        c=getchar();
    }
    while(isdigit(c)){
        x=(x<<1)+(x<<3)+(c^48);
        c=getchar();
    }
    return f?-x:x;
}

void gauss(){
    double eps=1e-6;
    for(int i=1;i<=n;++i){//enumerate the equation;
        int p=0;          //Record the position of the largest number; 
        double mx=0;      //Recording the largest number;
        for(int j=1;j<=n;++j)
            if(fabs(a[i][j])-eps>mx){
                mx=fabs(a[i][j]);p=j;//fabs() returns the absolute value of float; 
            }
        if(!p){
            printf("No Solution");
            return; 
        }
        w[i]=p;
        for(int j=1;j<=n;++j)
            if(i!=j){       //other equations
                double t=a[j][p]/a[i][p];
                for(int k=1;k<=n+1;++k)//n+1 is important
                    a[j][k]-=a[i][k]*t;
            }
    }
    for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]];
    for(int i=1;i<=n;++i) printf("%.2lf\n",v[i]);
}

int main(){
    n=rd();
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n+1;++j)
            a[i][j]=rd();
    gauss();
    return 0; 
}

浅谈高斯消元的实现和简单应用

标签：port   高斯消元   desc   模板   一个   rate   imp   script   AC   

原文地址：https://www.cnblogs.com/COLIN-LIGHTNING/p/8981923.html