线性代数 —— 高斯消元法

最新推荐文章于 2024-02-01 03:00:00 发布

Alex_McAvoy

最新推荐文章于 2024-02-01 03:00:00 发布

阅读量6w

点赞数 103

分类专栏： ——————线性代数—————— # 线性代数——高斯消元法

本文链接：https://blog.csdn.net/u011815404/article/details/88890702

版权

——————线性代数—————— 同时被 2 个专栏收录

6 篇文章 7 订阅

订阅专栏

线性代数——高斯消元法

6 篇文章 0 订阅

订阅专栏

【概述】

高斯消元法主要用于求解线性方程组，也可以求矩阵的秩、矩阵的逆等，是一个重要的数学方法。

其时间复杂度主要与方程组个数、方程组未知数个数有关，一般来说，时间复杂度为 O(n^3)

【线性方程组】

线性方程组：有多个未知数，且每个未知数的次数均为一次，这样多个未知数所组成的方程组。

其形式为： $\left\{\begin{matrix} a_{11}x_1+a_{12}x_2+...+a_{1n}x_n=b_1 \\ a_{21}x_1+a_{22}x_2+...+a_{2n}x_n=b_2 \\ ... \\ a_{n1}x_1+a_{n2}x_2+...+a_{nn}x_n=b_n \end{matrix}\right.$

记为矩阵形式，有：

$Ax=B,A=\begin{bmatrix}a_{11} &a_{12} &... &a_{1n} \\ a_{21}& a_{22} & ... & a_{2n}\\ ...& ... & ... & ...\\ a_{n1} & a_{n2} &... &a_{nn} \end{bmatrix},B= \begin{bmatrix} b_1\\b_2\\...\\b_n\end{bmatrix}$

【高斯消元法】

高斯消元法的基本思想是：通过一系列的加减消元运算，直到得到类似 kx=b 的式子，然后逐一回代求解 x 向量

1.解方程组过程

以以下线性方程组为例

$\left\{\begin{matrix} 3x+2y+z=6 \:\:\:\: (1)\\ 2x+2y+2z=4 \:\:\:\:(2)\\4x-2y-2z=2\:\:\:\:(3) \end{matrix}\right.$

首先进行消元操作：

将 (1) 除以 3，把 x 的系数化为 1，得： $x+\frac{2}{3}y+\frac{1}{3}z=2 \:\:\:\: (1)'$

再令 (2)-2*(1)'、(3)-4*(1)'，将 (2)、(3) 的 x 消去，得： $\left\{\begin{matrix} x+\frac{2}{3}y+\frac{1}{3}z=2 \:\:\:\: (1)'\\ 0x+\frac{2}{3}y+\frac{4}{3}z=0 \:\:\:\:(2)'\\0x-\frac{14}{3}y-\frac{10}{3}z=-6\:\:\:\:(3)' \end{matrix}\right.$

将令 (2)' 除以 2/3，将 y 系数化为 1，得： $y+2z=0 \:\:\:\:(2)''$

再令 (3)'-(-14/3)*(2)''，将 (3)' 的 y 消去，得： $\left\{\begin{matrix} x+\frac{2}{3}y+\frac{1}{3}z=2 \:\:\:\: (1)'\\ 0x+y+2z=0 \:\:\:\:(2)''\\0x+0y+\frac{18}{3}z=-6\:\:\:\:(3)'' \end{matrix}\right.$

由 (3)'' 可得： z=-1

最后进行回代操作：

将 z=-1 带回 (2)'' 可得： y=2

将 y=2 、 z=-1 带回 (1)' 可得： x=1

即结果为： $\left\{\begin{matrix} x=1\\ y=2\\z=-1 \end{matrix}\right.$

2.矩阵运算消元过程

同样以以下线性方程组为例

$\left\{\begin{matrix} 3x_1+2x_2+x_3=6 \:\:\:\: (1)\\ 2x_1+2x_2+2x_3=4 \:\:\:\:(2)\\4x_1-2x_2-2x_3=2\:\:\:\:(3) \end{matrix}\right.$

将其记为矩阵形式，有： $\begin{bmatrix} x & y & z & val \\ 3&2 &1 &6 \\2 &2 &2 &4 \\4 & -2&-2 &2 \end{bmatrix}$

消去 x，有： $\begin{bmatrix} x & y & z & val \\ 1&\frac{2}{3} &\frac{1}{3} &2 \\0 &\frac{2}{3} &\frac{4}{3} &0 \\0 & -\frac{14}{3}&-\frac{10}{3} &-6 \end{bmatrix}$

消去 y，有： $\begin{bmatrix} x & y & z & val \\ 1&\frac{2}{3} &\frac{1}{3} &2 \\0 &1&2 &0 \\0 & 0&\frac{18}{3} &-6 \end{bmatrix}$

回代得解： $\begin{bmatrix}1 \\ 2 \\ -1 \end{bmatrix}$

3.解的判断

无解：当消元完毕后，发现有一行系数都为 0，但是常数项不为 0，此时无解

例如： $\begin{bmatrix} x & y & z & val \\ 1&\frac{2}{3} &\frac{1}{3} &2 \\0 &1&2 &0 \\0 & 0&0 &-6 \end{bmatrix}$

多解：当消元完毕后，发现有多行系数、常数项均为 0，此时多解，有几行为全为 0，就有几个自由元，即变量的值可以任取，有无数种情况可以满足给出的方程组

例如： $\begin{bmatrix} x & y & z & val \\ 1&\frac{2}{3} &\frac{1}{3} &2 \\0 &0&0 &0 \\0 & 0&0 &0 \end{bmatrix}$ ，此时自由元个数为 2

【实现】

1.线性方程组整数类型解

int a[N][N];//增广矩阵
int x[N];//解集
bool freeX[N];//标记是否为自由变元
int GCD(int a,int b){
    return !b?a:GCD(b,a%b);
}
int LCM(int a,int b){
    return a/GCD(a,b)*b;
}
int Gauss(int equ,int var){//返回自由变元个数
    /*初始化*/
    for(int i=0;i<=var;i++){
        x[i]=0;
        freeX[i]=true;
    }

    /*转换为阶梯阵*/
    int col=0;//当前处理的列
    int row;//当前处理的行
    for(row=0;row<equ&&col<var;row++,col++){//枚举当前处理的行
        int maxRow=row;//当前列绝对值最大的行
        for(int i=row+1;i<equ;i++){//寻找当前列绝对值最大的行
            if(abs(a[i][col])>abs(a[maxRow][col]))
                maxRow=i;
        }
        if(maxRow!=row){//与第row行交换
            for(int j=row;j<var+1;j++)
                swap(a[row][j],a[maxRow][j]);
        }
        if(a[row][col]==0){//col列第row行以下全是0，处理当前行的下一列
            row--;
            continue;
        }

        for(int i=row+1;i<equ;i++){//枚举要删去的行
            if(a[i][col]!=0){
                int lcm=LCM(abs(a[i][col]),abs(a[row][col]));
                int ta=lcm/abs(a[i][col]);
                int tb=lcm/abs(a[row][col]);
                if(a[i][col]*a[row][col]<0)//异号情况相加
                    tb=-tb;
                for(int j=col;j<var+1;j++) {
                    a[i][j]=a[i][j]*ta-a[row][j]*tb;
                }
            }
        }
    }

    /*求解*/
    //无解：化简的增广阵中存在(0,0,...,a)这样的行，且a!=0
    for(int i=row;i<equ;i++)
        if (a[i][col]!=0)
            return -1;

    //无穷解: 在var*(var+1)的增广阵中出现(0,0,...,0)这样的行
    int temp=var-row;//自由变元有var-row个
    if(row<var)//返回自由变元数
        return temp;

    //唯一解: 在var*(var+1)的增广阵中形成严格的上三角阵
    for(int i=var-1;i>=0;i--){//计算解集
        int temp=a[i][var];
        for(int j=i+1;j<var;j++){
            if(a[i][j]!=0)
                temp-=a[i][j]*x[j];
        }
        if(temp%a[i][i]!=0)//有浮点数解，无整数解
            return -2;
        x[i]=temp/a[i][i];
    }
    return 0;
}

int main(){
    int equ,var;//equ个方程，var个变元
    while(scanf("%d%d",&equ,&var)!=EOF) {

        memset(a,0,sizeof(a));
        for(int i=0;i<equ;i++)//输入增广矩阵
            for(int j=0;j<var+1;j++)
                scanf("%d",&a[i][j]);


        int freeNum=Gauss(equ,var);//自由元个数
        if(freeNum==-1)//无解
            printf("无解\n");
        else if(freeNum==-2)//有浮点数解无整数解
            printf("无整数解\n");
        else if(freeNum>0){//有无穷多解
            printf("有无穷多解，自由变元个数为%d\n",freeNum);
            for(int i=0;i<var;i++){
                if(freeX[i])
                    printf("x%d是自由变元\n",i+1);
                else
                    printf("x%d=%d\n",i+1,x[i]);
            }
        }
        else{//有唯一解
            for(int i=0;i<var;i++)
                printf("x%d=%d\n",i+1,x[i]);
        }
        printf("\n");
    }
    return 0;
}

2.线性方程组浮点类型解

double a[N][N];//增广矩阵
double x[N];//解集
bool freeX[N];//标记是否为自由变元

int Gauss(int equ,int var){//返回自由变元个数
    /*初始化*/
    for(int i=0;i<=var;i++){
        x[i]=0;
        freeX[i]=true;
    }

    /*转换为阶梯阵*/
    int col=0;//当前处理的列
    int row;//当前处理的行
    for(row=0;row<equ&&col<var;row++,col++){//枚举当前处理的行
        int maxRow=row;//当前列绝对值最大的行
        for(int i=row+1;i<equ;i++){//寻找当前列绝对值最大的行
            if(abs(a[i][col])>abs(a[maxRow][col]))
                maxRow=i;
        }
        if(maxRow!=row){//与第row行交换
            for(int j=row;j<var+1;j++)
                swap(a[row][j],a[maxRow][j]);
        }
        if(fabs(a[row][col])<1e6){//col列第row行以下全是0，处理当前行的下一列
            row--;
            continue;
        }

        for(int i=row+1;i<equ;i++){//枚举要删去的行
            if(fabs(a[i][col])>1e6){
                double temp=a[i][col]/a[row][col];
                for(int j=col;j<var+1;j++) 
                    a[i][j]-=a[row][j]*temp;
                a[i][col]=0;
            }
        }
    }

    /*求解*/
    //无解
    for(int i=row;i<equ;i++)
        if(fabs(a[i][col])>1e6)
            return -1;

    //无穷解: 在var*(var+1)的增广阵中出现(0,0,...,0)这样的行
    int temp=var-row;//自由变元有var-row个
    if(row<var)//返回自由变元数
        return temp;

    //唯一解: 在var*(var+1)的增广阵中形成严格的上三角阵
    for(int i=var-1;i>=0;i--){//计算解集
        double temp=a[i][var];
        for(int j=i+1;j<var;j++)
            temp-=a[i][j]*x[j];
        x[i]=temp/a[i][i];
    }
    return 0;
}

3.模线性方程组

在解模线性方程组组时，当有唯一解时，需要对每个方程的唯一解不断的循环取模判断，直到找出能整除的为止

int a[N][N];//增广矩阵
int x[N];//解集
bool freeX[N];//标记是否为自由变元
int GCD(int a,int b){
    return !b?a:GCD(b,a%b);
}
int LCM(int a,int b){
    return a/GCD(a,b)*b;
}
int Gauss(int equ,int var){//返回自由变元个数
    /*初始化*/
    for(int i=0;i<=var;i++){
        x[i]=0;
        freeX[i]=true;
    }

    /*转换为阶梯阵*/
    int col=0;//当前处理的列
    int row;//当前处理的行
    for(row=0;row<equ&&col<var;row++,col++){//枚举当前处理的行
        int maxRow=row;//当前列绝对值最大的行
        for(int i=row+1;i<equ;i++){//寻找当前列绝对值最大的行
            if(abs(a[i][col])>abs(a[maxRow][col]))
                maxRow=i;
        }
        if(maxRow!=row){//与第row行交换
            for(int j=row;j<var+1;j++)
                swap(a[row][j],a[maxRow][j]);
        }
        if(a[row][col]==0){//col列第row行以下全是0，处理当前行的下一列
            row--;
            continue;
        }

        for(int i=row+1;i<equ;i++){//枚举要删去的行
            if(a[i][col]!=0){
                int lcm=LCM(abs(a[i][col]),abs(a[row][col]));
                int ta=lcm/abs(a[i][col]);
                int tb=lcm/abs(a[row][col]);
                if(a[i][col]*a[row][col]<0)//异号情况相加
                    tb=-tb;
                for(int j=col;j<var+1;j++) {
                    a[i][j]=((a[i][j]*ta-a[row][j]*tb)%MOD+MOD)%MOD;
                }
            }
        }
    }

    /*求解*/
    //无解：化简的增广阵中存在(0,0,...,a)这样的行，且a!=0
    for(int i=row;i<equ;i++)
        if (a[i][col]!=0)
            return -1;

    //无穷解: 在var*(var+1)的增广阵中出现(0,0,...,0)这样的行
    int temp=var-row;//自由变元有var-row个
    if(row<var)//返回自由变元数
        return temp;

    //唯一解: 在var*(var+1)的增广阵中形成严格的上三角阵
    for(int i=var-1;i>=0;i--){//计算解集
        int temp=a[i][var];
        for(int j=i+1;j<var;j++){
            if(a[i][j]!=0)
                temp-=a[i][j]*x[j];
            temp=(temp%MOD+MOD)%MOD;//取模
        }
        while(temp%a[i][i]!=0)//外层每次循环都是求a[i][i]，它是每个方程中唯一一个未知的变量
            temp+=MOD;//a[i][i]必须为整数，加上周期MOD
        x[i]=(temp/a[i][i])%MOD;//取模
    }
    return 0;
}

4.异或方程组

异或方程组是指形如 $\left\{\begin{matrix}a_{11}x_1\bigoplus a_{12}x_2 \bigoplus ...\bigoplus a_{1n}x_n=b_1 \\ a_{21}x_1\bigoplus a_{22}x_2 \bigoplus ...\bigoplus a_{2n}x_n=b_2 \\ ... \\a_{n1}x_1\bigoplus a_{n2}x_2 \bigoplus ...\bigoplus a_{nn}x_n=b_n \end{matrix}\right.$ 的方程组

对于 k=1...n，找到一个 a[i][k] 不为 0 的行 i，把它与第 k 行交换，用第 k 行去异或下面所有 a[i][j] 不为 0 的行 i，消去它们的第 k 个系数，这样就将原矩阵化成了上三角矩阵

由于最后一行只有一个未知数，这个未知数就已经求出来了，然后用它跟上面所有含有这个未知数的方程异或，以此类推即可以自下而上求出所有未知数。

要注意的是，对于开关问题，a[i][j]=1 表示第 i 个开关受第 j 个开关影响

int a[N][N];//增广矩阵
int x[N];//解集
int freeX[N];//自由变元
int Gauss(int equ,int var){//返回自由变元个数
    /*初始化*/
    for(int i=0;i<=var;i++){
        x[i]=0;
        freeX[i]=0;
    }

    /*转换为阶梯阵*/
    int col=0;//当前处理的列
    int num=0;//自由变元的序号
    int row;//当前处理的行
    for(row=0;row<equ&&col<var;row++,col++){//枚举当前处理的行
        int maxRow=row;//当前列绝对值最大的行
        for(int i=row+1;i<equ;i++){//寻找当前列绝对值最大的行
            if(abs(a[i][col])>abs(a[maxRow][col]))
                maxRow=i;
        }
        if(maxRow!=row){//与第row行交换
            for(int j=row;j<var+1;j++)
                swap(a[row][j],a[maxRow][j]);
        }
        if(a[row][col]==0){//col列第row行以下全是0，处理当前行的下一列
            freeX[num++]=col;//记录自由变元
            row--;
            continue;
        }

        for(int i=row+1;i<equ;i++){
            if(a[i][col]!=0){
                for(int j=col;j<var+1;j++){//对于下面出现该列中有1的行，需要把1消掉
                    a[i][j]^=a[row][j];
                }
            }
        }
    }

    /*求解*/
    //无解：化简的增广阵中存在(0,0,...,a)这样的行，且a!=0
    for(int i=row;i<equ;i++)
        if(a[i][col]!=0)
            return -1;

    //无穷解: 在var*(var+1)的增广阵中出现(0,0,...,0)这样的行
    int temp=var-row;//自由变元有var-row个
    if(row<var)//返回自由变元数
        return temp;

    //唯一解: 在var*(var+1)的增广阵中形成严格的上三角阵
    for(int i=var-1;i>=0;i--){//计算解集
        x[i]=a[i][var];
        for(int j=i+1;j<var;j++)
            x[i]^=(a[i][j]&&x[j]);
    }
    return 0;
}
int enumFreeX(int freeNum,int var){//枚举自由元，统计有解情况下1最少的个数
    int sta=(1<<(freeNum));//自由元的状态总数
    int res=INF;
    for(int i=0;i<sta;++i){//枚举状态
        int cnt=0;
        for(int j=0;j<freeNum;j++){//枚举自由元
            if(i&(1<<j)){
                cnt++;
                x[freeX[j]]=1;
            }else
                x[freeX[j]]=0;
        }
        for(int k=var-freeNum-1;k>=0;k--){//没有自由元的最下面一行
            int index=0;
            for(index=k;k<var;index++){//在当前行找到第一个非0自由元
                if(a[k][index])
                    break;
            }
            x[index]=a[k][var];
            for(int j=index+1;j<var;++j){//向后依次计算出结果
                if(a[k][j])
                    x[index]^=x[j];
            }
            cnt+=x[index];//若结果为1，则进行统计
        }
        res=min(res,cnt);
    }
    return res;
}
int main(){
    memset(a,0,sizeof(a));

    int equ,var;
    scanf("%d%d",&equ,&var);
    for(int i=0;i<equ;i++){//构造初始状态

    }
    for(int i=0;i<equ;i++)//最终状态
        scanf("%d",&a[i][var]);

    int freeNum=Gauss(equ,var);//获取自由元
    if(freeNum==-1)//无解
        printf("inf\n");
    else if(freeNum==0){//唯一解
        int res=0;
        for(int i=0;i<var;i++)
            res+=x[i];
        printf("%d\n",res);
    }
    else{//多个解
        int res=enumFreeX(freeNum,var);
        printf("%d\n",res);
    }
    return 0;
}