高斯-约旦消元法、列主元消去法 及其差异

高斯消元法求解线性方程组：

对增广矩阵做相应的初等行变换，将右边系数矩阵化为单位矩阵后，最右边一列就是原线性方程组的解。

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

1、高斯-约旦(G-J)消元法

直接消成对角矩阵
（1）从第一行开始，先顺序寻找到主元（不能为0），并将主元通过行变换移动到主对角线上
（2）将主元所在的行内的所有元素除以主元，使得主元化为1（也可以放在第四部之后）
（3）其他行减去 主元所在行 乘以一定的倍数，使得主元所在的列内除主元外的其他元素为0
（4）在剩下的矩阵范围内寻找下一个主元，重复以上步骤

#include<algorithm>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
using namespace std;
//#define int long long
#define ll long long
#define uql unsigned long long
#define pii pair<int,int>
#define mid ((l + r)>>1)
#define chl (root<<1)
#define chr (root<<1|1)
#define lowbit(x) ( x&(-x) )
const int manx = 100+10;
const int manx2 = 4e7 + 10;
const ll INF = 1e18;
const ll mod = 1000000007;
const double eps=1e-10;
//运算过程中精度丢失导致本来应当为0的值的机器表示不为0,所以用abs(x)<eps 来判0
int n;
double f[manx][manx];
int Gauss()
{
    for(int i=1;i<=n;i++){
        for(int j=i;j<=n;j++){/********对应步骤（1）********/
            if(fabs(f[j][i])>eps){
                swap(f[i],f[j]);
                break;
            }
        }
        if(fabs(f[i][i])<eps){//若该列都为0，则无解
            puts("No Solution");
            return 0;
        }
        for(int j=1;j<=n;j++){/********对应步骤（3）********/
            if(i==j)continue;
            double temp=f[j][i]/f[i][i];
            for(int k=i;k<=n+1;k++)
                f[j][k]-=temp*f[i][k];
        }
    }
    for(int i=1;i<=n;i++)//步骤（2）也可以最后再将主对角线上的数化为1
        f[i][n+1]/=f[i][i];
    return 1;
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&f[i][j]);
    if(!Gauss())return 0;
    for(int i=1;i<=n;i++)
        printf("%.2f\n",f[i][n+1]);
    return 0;
}

2、列主元消去法

先消成上三角矩阵 ，然后回带求解
（1）从第一列行开始，选出当前列中绝对值最大的元素作为这一列的主元，并将主元通过行变换移动到主对角线上
（2）将主元所在的行内的所有元素除以主元，使得主元化为1（也可以放在第四部之后）
（3）将剩下矩阵范围的其他行减去 主元所在行 乘以一定的倍数，使得剩下矩阵 主元所在的列内除主元外的其他元素为0
（4）在剩下的矩阵范围内寻找下一个主元，重复以上步骤
（5）回带求解

#include<algorithm>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
using namespace std;
//#define int long long
#define ll long long
#define uql unsigned long long
#define pii pair<int,int>
#define mid ((l + r)>>1)
#define chl (root<<1)
#define chr (root<<1|1)
#define lowbit(x) ( x&(-x) )
const int manx = 100+10;
const int manx2 = 4e7 + 10;
const ll INF = 1e18;
const ll mod = 1000000007;
const double eps=1e-18;
int n;
double f[manx][manx];
int Gauss()
{
    int Max;
    for(int i=1;i<=n;i++){
        Max=i;
        for(int j=i+1;j<=n;j++){/********对应步骤（1）********/
            if(f[Max][i]>fabs(f[j][i]))
                Max=j;
        }
        if(fabs(f[Max][i])<eps){
            puts("No Solution");
            return 0;
        }
        swap(f[Max],f[i]);
        for(int j=i+1;j<=n;j++){/********对应步骤（3）********/
            double temp=f[j][i]/f[i][i];
            for(int k=i;k<=n+1;k++)
                f[j][k]-=temp*f[i][k];
        }
    }
    for(int i=n;i>=1;i--){
        for(int j=i+1;j<=n;j++)/********对应步骤（5）回带求解********/
            f[i][n+1]-=f[i][j]*f[j][n+1];
        f[i][n+1]/=f[i][i];/********对应步骤（2）********/
    }
    return 1;
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&f[i][j]);
    if(!Gauss())return 0;
    for(int i=1;i<=n;i++)
        printf("%.2f\n",f[i][n+1]);
    return 0;
}

3、几种消元法对比

• 高斯-约旦消元法：没有回带过程 ，代码简单；容易理解。但因为每一次都要将其他行消元，效率比列主元低；且顺序寻找主元会出现主对角线元素的绝对值比同列对角线以下元素的绝对值还小的情形，导致消元过程用绝对值极小的数做分母的极端现象，以至求解过程中的误差不断积累并扩大（所以也可以在这里选择较大的数作为主元减小误差）。
• 列主元消去法：多一个回带过程。但效率和精度比高斯-约旦消元法高
• 全选主元消元法：精度比列主元消去法更高，但全选主元对矩阵还进行了列变换，列变换使得回代时未知数的位置变换，需要在进行列变换的同时将未知量的顺序同时进行变换，增加了较多的计算量，实际用时接近列选主元的两倍。

4、double的精度问题

（1）double的存储方式：

（2）有效位：15~16位
因为最后一位有舍入误差，所以能确保至少正确记录前15位：
（根据存储格式，64位浮点数底数有52位，加一个隐藏位共53位，64位浮点数最大值有效是:2⁵³= 9007199254740992）

十进制数转二进制数：整数部分的转换：除2取余；小数部分的转换：乘2取余

有时候，在转换中，由于整数部分的二进制数位已经大于53位，或二进制小数的某些位会周而复始地重复，以致无穷。但计算机的表示是有限的，所以在计算机内只能截取到某个精度
（误差大小：浮点数的有效数字）

（3）高斯消元中主元绝对值大小对精度的影响

若选用较小数作为主元，可能会导致其它行所乘的系数temp过大，而此时计算机中所存储的f[i][k]是有误差的，temp个f[i][k]相加，temp越大，因为f[i][k]带来的误差就会越大

参考：
1、浅析高斯消元法如何求解线性方程组
2、高斯列主元消去法解线性方程组
3、列主元消去法与全主元消去法
4、选不同主元的高斯消元
5、什么是高斯消元法的精度？
6、浮点数在计算机中的存储方式