高斯消元

Posted 2021-02-15 caijiaming

学习博客：https://blog.csdn.net/sslz_fsy/article/details/82021887

题目背景
Gauss消元

题目描述
给定一个线性方程组，对其求解

输入输出格式
输入格式：
 

第一行，一个正整数 n

第二至 n+1 行，每行 n+1 个整数，为 a1,a2...an ，b代表一组方程。

 

输出格式：
 

共n行，每行一个数，第 i行为 xi? （保留2位小数）

如果不存在唯一解，在第一行输出"No Solution".

 

输入输出样例
输入样例#1： 复制

3
1 3 4 5
1 4 7 3
9 3 2 2
输出样例#1： 复制

-0.97
5.18
-2.39
以样例为例,我们的矩阵为

1 3 4 5
1 4 7 3
9 3 2 2
 我们的最终目标是

1 0 0 ?
0 1 0 ?
0 0 1 ?
所以说,我们枚举第i行,然后把第i列保留为1

为了实现,有以下步骤

1.将第i个系数最大的放在第i行
什么意思,假如i=1,第i个系数最大的是9,于是我们交换第一行和第三行

9 3 2 2
1 4 7 3
1 3 4 5
代码

int k=i;
for(int j=i+1;j<=n;j++) if(fabs(a[j][i])>fabs(a[k][i])) k=j//找到这一行
for(int j=i;j<=n+1;j++) swap(a[i][j],a[k][j]);//将每一个元素交换

 

2.把第i行的第i个系数变为1 
我们要将9x+3y+2z=2,x的系数变为1,那每个数都除以9呗

1 1/3 2/9 2/9
1 4 7 3
1 3 4 5
代码

int res=a[i][i]
for(int j=i;j<=n+1;j++) a[i][j]/=res;

 

因为前面都是0了,所以从第i列开始除就行了 

3.把这一列回代到原方程
 就是将第i个系数消掉

怎么消呢

看看我们现在的矩阵

x+1/3y+2/9z=2/9

x+4y+7z=3

那么(x-x)+(4-1/3)y+(7-2/9)z=3-2/9,得

11/3y+61/9z=25/9

那万一系数不相同呢

假如说有个方程3x+2y+z=1

那我们就将要代入的方程成3(3x+y+2/3z=2/3),再依次相减,就能保证减后的系数为0了

1 1/3 2/9 2/9
0 4-1/3 7-2/9 3-2/9
0 3-1/3 4-2/9 5-2/9
代码

for(k=1;k<=n;k++)//回代的过程
if(k!=i)
{
res=a[k][i];//原方程每一项都要乘的,为了消掉一个系数
for(int j=i;j<=n+1;j++) a[k][j]-=a[i][j]*res;
}

 

附每次消元后矩阵的情况
 

 

 P.S 判断解的情况
如果我们到第i行,发现第i个的系数为0,就有无数多种解

因为这种情况,就是0x=m的情况,x为任意实数

代码

if(fabs(a[i][i])<1e-8) return 0;

 

最后看总代码：

#include<iostream>
#include<cmath>
using namespace std;
const int maxn=100+50;
const double eps=1e-8;
double a[maxn][maxn];
int Gauss(int N)
{
    for(int i=1;i<=N;i++)//遍历N*N矩阵中的 a[i][i]  也就是这个矩阵中的主对角线上的元素
    {
        int max_r=i;
        for(int j=i+1;j<=N;j++)//下面的行中 找到你绝对值最大的系数在哪一行
        {
            if(fabs(a[j][i])>fabs(a[i][i])) max_r=j;
        }
        if(fabs(a[max_r][i])<eps) return 0;//最大的一行系数都为0的话  那么代表这个方程有无数多个解
        if(max_r!=i)//最大的行不在当前行  交换使得在当前行  其实不交换也是可以的
        {
            for(int j=i;j<=N+1;j++)
            {
                swap(a[i][j],a[max_r][j]);
            }
        }
        double res=a[i][i];
        for(int j=i;j<=N+1;j++) a[i][j]/=res;//
        for(int j=1;j<=N;j++)//遍历每一行
        {
            if(j!=i)//不是当前行
            {
                res=a[j][i];
                for(int k=i;k<=N+1;k++)//遍历这一行从第i列开始的所有数 因为前面的列都是0
                {
                    a[j][k]-=a[i][k]*res;
                }
            }
        }
        for(int j=1;j<=N;j++) {for(int k=1;k<=N+1;k++) {cout<<a[j][k]<<" ";}cout<<endl;}
    }
    return 1;
}
int main()
{
    int N;
    cin>>N;//N 行 N+1列
    for(int i=1;i<=N;i++)
    {
        for(int j=1;j<=N+1;j++) cin>>a[i][j];
    }
    if(Gauss(N))
    {
        for(int i=1;i<=N;i++) cout<<a[i][N+1]<<endl;
    }
    else cout<<"No solution"<<endl;
    return 0;
}