高斯消元法求解方程

Posted 2020-09-03 穿越银河系

#include<iostream>
#include<math.h>
#include<vector>
#include<bits/stdc++.h>
using namespace std;
typedef vector<double> vec;
typedef vector<vec> mat;
const double EPS=1E-8;

//求解Ax=b,其中A是方针
//当方程组无解或者有无穷多解时，返回一个长度为0的数组
vec gauss_jordan(const mat& A,const vec& b)
{
    int n=A.size();
    mat B(n,vec(n+1));
    for(int i=0; i<n; i++)
        for(int j=0; j<n; j++)
            B[i][j]=A[i][j];
    //把b存放在A的右边方便一起处理
    for(int i=0; i<n; i++) B[i][n]=b[i];

    for(int i=0; i<n; i++) //第i个未知数
    {
        //把正在处理的未知数绝对值最大的式子换到第i行
        int pivot=i;
        for(int j=i; j<n; j++)
        {
            if(abs(B[pivot][j]<abs(B[i][j]))) pivot=j;
        }
        swap(B[i],B[pivot]);

        //无解或者有无穷多解
        if(abs(B[i][i])<EPS) return vec();

        //把正在处理的未知数系数变成1
        for(int j=i+1; j<=n; j++) B[i][j]/=B[i][i];

        for(int j=0; j<n; j++)
        {
            if(i!=j)
            {
                //从第j个式子中消去第i个未知数
                for(int k=i+1; k<=n; k++) B[j][k]-=B[j][i]*B[i][k];
            }
        }
    }
    vec x(n);
    for(int i=0; i<n; i++) x[i]=B[i][n];
    return x;
}

int main()
{
    //等式右侧
    vec ve(3);
    ve[0]=6;
    ve[1]=12;
    ve[2]=21;
    //等式左侧系数矩阵
    mat ma(3,vec(3));
    ma[0][0]=1;
    ma[0][1]=-2;
    ma[0][2]=3;
    ma[1][0]=4;
    ma[1][1]=-5;
    ma[1][2]=6;
    ma[2][0]=7;
    ma[2][1]=-8;
    ma[2][2]=10;

    vec result=gauss_jordan(ma,ve);
    //存放在右边的b就是答案
    for(int i=0;i<3;i++) cout<<result[i]<<endl;
    return 0;
}

View Code

求解的过程中，在消去某个未知数时，打算保留该未知数的式子的对应未知数系数可能是0，在这种情况下，只需要调整方程的顺序，使得对应的系数不为0即可。

为了减小误差，应该选择要消去的未知数系数的绝对值尽可能大的方程。