数学板块学习之高斯消元

最新推荐文章于 2022-03-02 08:31:20 发布

Toniht

最新推荐文章于 2022-03-02 08:31:20 发布

阅读量152

点赞数

分类专栏：知识点

本文链接：https://blog.csdn.net/qq_41746268/article/details/90069999

版权

知识点专栏收录该内容

49 篇文章 1 订阅

订阅专栏

线性方程组

①方程组
$\begin{cases} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1 & \\ a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2\\ \cdots\cdots\\ a_{m1}x_1+a_{m2}x_2+\cdots+a_{mn}x_n=b_m\\ \end{cases}$
称为 $m$ 个方程， $n$ 个未知量的线性方程组
其中 $x_1,x_2,\cdots,x_n$ 是未知数， $m$ 为方程个数， $b_1,b_2,\cdots,b_m$ 为常数项
如果常数项 $b_1=b_2=\cdots=b_m=0$ ，则称该方程组为齐次线性方程组。相应来讲，如果 $b_1,b_2,\cdots,b_m$ 不全为0，则称该方程组为非齐次线性方程组
②线性方程组的矩阵表示。
由线性方程组
$\begin{cases} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1 & \\ a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2\\ \cdots\cdots\\ a_{m1}x_1+a_{m2}x_2+\cdots+a_{mn}x_n=b_m\\ \end{cases}$
该线性方程组的系数构成的 $m \times n$ 矩阵
$A=\left[ \begin{matrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots\\ a_{m1} & a_{m2} & \cdots & a_{mn} \end{matrix} \right]$
矩阵 $A$ 称为该线性方程组的系数矩阵。
由线性方程组的系数矩阵和常数项构成的 $m \times (n + 1)$ 矩阵
$\overline{A}=\left[ \begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1\\ a_{21} & a_{22} & \cdots & a_{2n} & b_2\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ a_{m1} & a_{m2} & \cdots & a_{mn} & b_m \end{array} \right]$
矩阵 $\overline{A}$ 称为该线性方程组的增广矩阵
如果令 $x=[x_1,x_2,\cdots,x_n]^T,b=[b_1,b_2,\cdots,b_n]^T$ ，利用矩阵的乘法，我们可以将原线性方程组简写为 $A x = b$ 。

高斯消元

高斯消元主要用来求线性方程组的解，也可用来求解矩阵的秩，矩阵的逆
时间复杂度为 $O(n^3)$ ，主要与方程组的个数与未知数个数有关
一直以来我们学习的方程组解法就是加减消元，消去未知数，然后求出一个未知数后回带。这也就是利用初等变换做到的，因为做完初等变换的方程组和原方程组是同解，所以，我们对其增广矩阵做初等变换也是同解，而因为齐次线性方程组的常数项全为0,所以我们仅需要对系数矩阵做初等变换即可

对线性方程组做初等变换是为了将其化为该种形式
$\begin{cases} a'_{11}x_1+a'_{12}x_2+a'_{13}x_3+\cdots+a'_{1n}x_n=b'_{1}\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ a'_{22}x_2+a'_{23}x_3+\cdots+a'_{2n}x_n=b'_2\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ a'_{33}x_3+\cdots+a'_{3n}x_n=b'_3\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ a'_{kk}x_k+\cdots+a'_{mn}x_n=b'_m\\ \end{cases}$
该方程组中每一个方程都比上一个方程少一个未知数，这种方程称为阶梯型方程组，每行第一个未知数称为主元，其他未知数称为自由变量。将线性方程组化为同解的阶梯型方程组的过程即为高斯消元。也容易得知，利用高斯消元求解线性方程组就等价于利用初等变换将线性方程组的增广矩阵化为阶梯型矩阵。
例如利用高斯消元解方程组
$\begin{cases} x_1+x_2+x_3=3 & \\ x_1+2x_2+4x_3=7\\ x_1+3x_2+9x_3=13\\ \end{cases}$
利用高斯消元对其增广矩阵进行变化：
$\left[ \begin{array}{ccc|c} 1 & 1 & 1 & 3\\ 1 & 2 & 4 & 7\\ 1 & 3 & 9 &13 \end{array} \right]\Rightarrow\left[ \begin{array}{ccc|c} 1 & 1 & 1 & 3\\ 0 & 1 & 3 & 4\\ 0 & 2 & 8 &10\\ \end{array} \right]\Rightarrow\left[ \begin{array}{ccc|c} 1 & 1 & 1 & 3\\ 0 & 1 & 3 & 4\\ 0 & 1 & 4 & 5\\ \end{array} \right]\Rightarrow\left[ \begin{array}{ccc|c} 1 & 1 & 1 & 3\\ 0 & 1 & 3 & 4\\ 0 & 0 & 1 & 1\\ \end{array} \right]$
再将增广矩阵化为线性方程组即可求出线性方程组的解

注意事项：
①消元后发现有一行系数为0，但是常数项不为0，则方程组无解例如
$\left[ \begin{array}{ccc|c} 1 & 1 & 1 & 3\\ 0 & 1 & 3 & 4\\ 0 & 0 & 0 & 1\\ \end{array} \right]$
②消元后发现有好几行的系数为0且常数项也为0，则方程组有多个解，其中有多少行全为0既有多少个自由元，例如
$\left[ \begin{array}{ccc|c} 1 & 2 & 1 & 3\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ \end{array} \right]$

大佬的高斯消元的模板:

    const int MAXN=50;
    int a[MAXN][MAXN];//增广矩阵
    int x[MAXN];//解集
    bool free_x[MAXN];//标记是否是不确定的变元
    int gcd(int a,int b){
        if(b == 0) return a; else return gcd(b,a%b);
    }
    inline int lcm(int a,int b){
        return a/gcd(a,b)*b;//先除后乘防溢出
    }
    // 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解，但无整数解，
    //-1表示无解，0表示唯一解，大于0表示无穷解，并返回自由变元的个数)
    //有equ个方程，var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var.
    int Gauss(int equ,int var){
        int i,j,k;
        int max_r;// 当前这列绝对值最大的行.
        int col;//当前处理的列
        int ta,tb;
        int LCM;
        int temp;
        int free_x_num;
        int free_index;

        for(int i=0;i<=var;i++){
            x[i]=0;
            free_x[i]=true;
        }

        //转换为阶梯阵.
        col=0; // 当前处理的列
        for(k = 0;k < equ && col < var;k++,col++){// 枚举当前处理的行.
        // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
            max_r=k;
            for(i=k+1;i<equ;i++){
                if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
            }
            if(max_r!=k){// 与第k行交换.
                for(j=k;j<var+1;j++) swap(a[k][j],a[max_r][j]);
            }
            if(a[k][col]==0){// 说明该col列第k行以下全是0了，则处理当前行的下一列.
                k--;
                continue;
            }
            for(i=k+1;i<equ;i++){// 枚举要删去的行.
                if(a[i][col]!=0){
                    LCM = lcm(abs(a[i][col]),abs(a[k][col]));
                    ta = LCM/abs(a[i][col]);
                    tb = LCM/abs(a[k][col]);
                    if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加
                    for(j=col;j<var+1;j++){
                        a[i][j] = a[i][j]*ta-a[k][j]*tb;
                    }
                }
            }
        }
        // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
        for (i = k; i < equ; i++){ // 对于无穷解来说，如果要判断哪些是自由变元，那么初等行变换中的交换就会影响，则要记录交换.
            if (a[i][col] != 0) return -1;
        }
        // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行，即说明没有形成严格的上三角阵.
        // 且出现的行数即为自由变元的个数.
        if (k < var){
            return var - k; // 自由变元有var - k个.
        }
        // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
        // 计算出Xn-1, Xn-2 ... X0.
        for (i = var - 1; i >= 0; i--){
            temp = a[i][var];
            for (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(void){
        int i, j;
        int equ,var;
        while (scanf("%d %d", &equ, &var) != EOF){
            memset(a, 0, sizeof(a));
            for (i = 0; i < equ; i++){
                for (j = 0; j < var + 1; j++){
                    scanf("%d", &a[i][j]);
                }
            }
            int free_num = Gauss(equ,var);
            if (free_num == -1) printf("无解!\n");
            else if (free_num == -2) printf("有浮点数解，无整数解!\n");
            else if (free_num > 0){
                printf("无穷多解! 自由变元个数为%d\n", free_num);
                for (i = 0; i < var; i++){
                    if (free_x[i]) printf("x%d 是不确定的\n", i + 1);
                    else printf("x%d: %d\n", i + 1, x[i]);
                }
            }else{
                for (i = 0; i < var; i++){
                    printf("x%d: %d\n", i + 1, x[i]);
                }
            }
            printf("\n");
        }
        return 0;
    }

浮点型

const int MAXN=220;
double a[MAXN][MAXN],x[MAXN];//方程的左边的矩阵和等式右边的值, 求解之后 x存的就是结果
// 注意a不是增广矩阵!!!,x开始存的是增广矩阵的b,消元后x是答案
//方程数和未知数个数
/*
* 返回 0 表示无解,1 表示有解
*/
int Gauss(int equ,int var){
    int i,j,k,col,max_r;
    for(k = 0,col = 0;k < equ && col < var; k++, col++){
        max_r = k;
        for(i = k+1;i < equ; i++) 
            if(fabs(a[i][col])>fabs(a[max_r][col]))
                max_r=i;
        if(fabs(a[max_r][col]) < eps) return 0;    //第一个最大的项为0,代表无解
        if(k != max_r){       //交换最大行和当前行
            for(j = col;j < var; j++)
                swap(a[k][j],a[max_r][j]);
            swap(x[k],x[max_r]);
        }
        
        x[k] /= a[k][col];    
        for(j = col+1;j < var; j++) a[k][j]/=a[k][col]; 
        a[k][col] = 1;
        for(i = 0;i < equ; i++)
            if(i != k){
                double rate = a[i][col];
                x[i] -= x[k]*rate;
                for(j = col;j < var; j++) a[i][j] -= a[k][j]*rate;
                    
        }
    }
    return 1;
}

异或方程组

int equ,var;
int a[maxn][maxn];
int x[maxn];
int free_x[maxn];
int free_num;

int Gauss(){
    int max_r,col,k;
    free_num = 0;
    for(k = 0,col = 0; k < equ && col < var; ++k,++col){
        max_r = k;
        for(int i = k+1; i < equ; ++i){
            if(abs(a[i][col]) > abs(a[max_r][col]))
                max_r = i;
        }
        if(a[max_r][col] == 0){
            k--;
            free_x[free_num++] = col;
            continue;
        }
        if(max_r != k){
            for(int j = col; j < var+1; ++j){
                swap(a[k][j],a[max_r][j]);
            }
        }
        for(int i = k + 1; i < equ; ++i){
            if(a[i][col] != 0){
                for(int j = col; j < var+1; ++j){
                    a[i][j] ^= a[k][j];
                }
            }
        }
        
    }
    for(int i = k; i < equ; ++i){
        if(a[i][col] != 0)return -1;
    }
    if(k < var)return var-k;
    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;
}