zoukankan      html  css  js  c++  java
  • 【线性代数】高斯消元


    普通代数式的消元法

    考虑一个普通的二元一次方程组

    [left { egin{aligned} x_1+x_2&=6 otag \ x_1-x_2&=2 otag end{aligned} ight. ]

    根据正常的数学知识,我们直到应该先消去一个未知数,解出另一个未知数,再去用另一个未知数表示消去的未知数

    故对于上述((1))式加上((2))式,得(2x_1=8),即(x_1=4),代回原式得(x_2=6-x_1=x_1-2=2)

    消元法遵循下列三个规则(懂的都懂):

    • 两方程互换,解不变
    • 某一方程乘上非零数(k),解不变
    • 用某一方程加上(或减去)另外一方程,解不变

    所以普通的消元法主要用于二元一次方程组的求解上,但对于多元方程组而言,普通消元法则表现出明显的力不从心(以前的数学题怎么都那么难)




    线性代数的高斯消元法

    德国数学家高斯分析了普通消元法并得出以下四个结论:

    • 在消元法中,参与计算和发生改变的是方程中各变量的系数
    • 各变量并未参与计算,且没有发生改变
    • 可以利用系数的位置表示变量,从而省略变量
    • 在计算中将变量简化省略,方程的解不变

    根据这些结论,便能通过对增广矩阵的运用来求出多元方程组的通解,这便是高斯消元法

    高斯消元法大致可分为以下五个步骤

    1. 增广矩阵行初等行变换为行最简形
    2. 还原线性方程组
    3. 求解第一个变量
    4. 补充自由未知量
    5. 列表示方程组通解

    举例说明,考虑下列多元一次方程组

    [left { egin{aligned} x_1+3x_2+x_4&=17 otag \ x_2+2x_4&=10 otag \ 2x_2+4x_4&=20 otag end{aligned} ight. ]

    1、增广矩阵行初等行变换为行最简形

    先将方程组以矩阵的方式表示

    [left ( egin{matrix} 1 & 3 & 0 & 1\ 0 & 1 & 0 & 2\ 0 & 2 & 0 & 4\ end{matrix} middle | egin{matrix} 17 \ 10 \ 20 end{matrix} ight ) ]

    根据增广矩阵的变换规则,将其变成行最简形

    [xrightarrow{r_3-2r_2} left ( egin{matrix} 1 & 3 & 0 & 1\ 0 & 1 & 0 & 2\ 0 & 0 & 0 & 0\ end{matrix} middle | egin{matrix} 17 \ 10 \ 0 end{matrix} ight ) \ xrightarrow{r_1-3r_2} left ( egin{matrix} 1 & 0 & 0 & -5\ 0 & 1 & 0 & 2\ 0 & 0 & 0 & 0\ end{matrix} middle | egin{matrix} -13 \ 10 \ 0 end{matrix} ight ) ]

    2、还原线性方程组

    即将增广矩阵形式重新转回方程组形式

    [left { egin{aligned} x_1-5x_4&=-13 otag \ x_2+2x_4&=10 otag end{aligned} ight . ]

    3、求解第一个变量

    对于每个方程,将其第一个变量用其余变量表示出来

    [left { egin{aligned} x_1&=5x_4-13 otag \ x_2&=-2x_4+10 otag end{aligned} ight. ]

    4、补充自由未知量

    由于自由未知量不受约束,所以只能是(x_i=x_i)

    所以将未出现的变量以上述形式表示出来

    [left { egin{aligned} x_1&=5x_4-13 otag \ x_2&=-2x_4+10 otag \ x_3&=x_3 otag \ x_4&=x_4 otag end{aligned} ight. ]

    5、列表示方程组通解

    [egin{aligned} left ( egin{matrix} x_1\x_2\x_3\x_4 end{matrix} ight ) =& left ( egin{matrix} 0\0\1\0 end{matrix} ight ) x_3 + left ( egin{matrix} 5\-2\0\1 end{matrix} ight ) x_4 + left ( egin{matrix} -13\10\0\0 end{matrix} ight ) \ =& left ( egin{matrix} 0\0\1\0 end{matrix} ight ) C_1 + left ( egin{matrix} 5\-2\0\1 end{matrix} ight ) C_2 + left ( egin{matrix} -13\10\0\0 end{matrix} ight ) end{aligned} ]

    其中(C_1,C_2in R)


    经过以上步骤,便得出了多元一次方程组的通解




    高斯消元法的运用


    行列式求值

    对于行列式,存在公式

    [D=|A|=sum_{p_1p_2dots p_n} (-1)^r a_{(1,p_1)} a_{(2,p_2)} a_{(3,p_3)} dots a_{(n,p_n)} ]

    其中(r)表示({p_1,p_2,p_3dots p_n})内逆序对对数

    由于高斯消元法可以让我们得到行最简形矩阵,即一个对角线矩阵

    此矩阵的行列式由对角线元素之积所决定

    其正负符号可由交换行的数量来确定

    故可以在(O(n^3))时间内求出给定行列式的值


    例题:NOJ 1035 - 行列式求值

    #include<bits/stdc++.h>
    using namespace std;
    
    struct matrix
    {
        typedef double Type;
        Type mat[1050][1050];
        int n,m;
        void init(int _n,int _m)
        {
            for(int i=1;i<=_n;i++)
                for(int j=1;j<=_m;j++)
                    mat[i][j]=0;
            n=_n,m=_m;
        }
        Type getVal()
        {
            Type det=1;
            for(int i=1;i<=n;i++)
            {
                int k=i;
                for(int j=i+1;j<=n;j++)
                    if(abs(mat[j][i])>abs(mat[k][i]))
                        k=j;
                if(mat[k][i]==0)
                {
                    det=0;
                    break;
                }
                swap(mat[i],mat[k]);
                if(i!=k)
                    det=-det;
                det*=mat[i][i];
                for(int j=i+1;j<=n;j++)
                    mat[i][j]/=mat[i][i];
                for(int j=1;j<=n;j++)
                    if(j!=i&&mat[j][i]!=0)
                        for(int k=i+1;k<=n;k++)
                            mat[j][k]-=mat[i][k]*mat[j][i];
            }
            return det;
        }
    }mat;
    
    int main()
    {
        int n;
        while(scanf("%d",&n)!=EOF)
        {
            mat.n=mat.m=n;
            for(int i=1;i<=n;i++)
                for(int j=1;j<=n;j++)
                    scanf("%lf",&mat.mat[i][j]);
            printf("%.2f
    ",mat.getVal());
        }
        return 0;
    }
    


    生成树计数

    生成树的计数一般采用Matrix-Tree Theorem解决——

    对于一张无向图,其生成树个数就等于这副图的拉普拉斯矩阵中任意一个代数余子式的值

    建立拉普拉斯矩阵,通过高斯消元化成上三角型再进行求值

    (行吧我也只会套模板了......)


    例题:SPOJ HIGH-Highways

    #include<bits/stdc++.h>
    using namespace std;
    
    struct matrix
    {
        typedef long long Type;
        Type mat[15][15];
        int n,m;
        void init(int _n,int _m)
        {
            for(int i=1;i<=_n;i++)
                for(int j=1;j<=_m;j++)
                    mat[i][j]=0;
            n=_n,m=_m;
        }
        void addedge(int a,int b)
        {
            mat[a][a]++;
            mat[b][b]++;
            mat[a][b]--;
            mat[b][a]--;
        }
        Type gauss()
        {
            Type res=1;
            for(int i=1;i<n;i++)
            {
                for(int j=i+1;j<n;j++)
                {
                    while(mat[j][i])
                    {
                        Type t=mat[i][i]/mat[j][i];
                        for(int k=i;k<n;k++)
                            mat[i][k]=mat[i][k]-t*mat[j][k];
                        swap(mat[i],mat[j]);
                        res=-res;
                    }
                }
                res=res*mat[i][i];
            }
            return res;
        }
    }mat;
    
    void solve()
    {
        int n,m,u,v;
        scanf("%d%d",&n,&m);
        mat.init(n,n);
        while(m--)
        {
            scanf("%d%d",&u,&v);
            mat.addedge(u,v);
        }
        printf("%lld
    ",mat.gauss());
    }
    int main()
    {
        int T;
        scanf("%d",&T);
        while(T--)
            solve();
        return 0;
    }
    



    上面用到的板子

    struct matrix
    {
        typedef long long Type;
        Type mat[1050][1050];
        int n,m;
        void init(int _n,int _m)
        {
            for(int i=1;i<=_n;i++)
                for(int j=1;j<=_m;j++)
                    mat[i][j]=0;
            n=_n,m=_m;
        }
        Type getVal()
        {
            Type det=1;
            for(int i=1;i<=n;i++)
            {
                int k=i;
                for(int j=i+1;j<=n;j++)
                    if(abs(mat[j][i])>abs(mat[k][i]))
                        k=j;
                if(mat[k][i]==0)
                {
                    det=0;
                    break;
                }
                swap(mat[i],mat[k]);
                if(i!=k)
                    det=-det;
                det*=mat[i][i];
                for(int j=i+1;j<=n;j++)
                    mat[i][j]/=mat[i][i];
                for(int j=1;j<=n;j++)
                    if(j!=i&&mat[j][i]!=0)
                        for(int k=i+1;k<=n;k++)
                            mat[j][k]-=mat[i][k]*mat[j][i];
            }
            return det;
        }
        void addedge(int a,int b) //生成树计数用
        {
            mat[a][a]++;
            mat[b][b]++;
            mat[a][b]--;
            mat[b][a]--;
        }
        Type gauss() //求矩阵K的n-1阶顺序主子式
        {
            Type res=1;
            for(int i=1;i<n;i++)
            {
                for(int j=i+1;j<n;j++)
                {
                    while(mat[j][i])
                    {
                        Type t=mat[i][i]/mat[j][i];
                        for(int k=i;k<n;k++)
                            mat[i][k]=mat[i][k]-t*mat[j][k];
                            //mat[i][k]=(mat[i][k]-t*mat[j][k]+mod)%mod;
                        swap(mat[i],mat[j]);
                        res=-res;
                    }
                }
                res=res*mat[i][i];
                //res=res*mat[i][i]%mod;
            }
            return res;
            //return (res+mod)%mod;
        }
    }mat;
    



    参考及引用:OI Wiki / Math / Gauss


  • 相关阅读:
    python一次性解压多层嵌套zip压缩包
    zabbix自定义监控项、添加图形、设置触发器、远程执行命令
    制作自己的docker镜像
    vue-router的Import() 异步加载模块问题的解决方案
    Idea Maven配置镜像地址总出问题,麻烦检查一下这里!(最新)
    杂谈:开发人员如何进行复杂业务的学习?让boss刮目相看
    Spring Boot 2.x 基础案例:整合Dubbo 2.7.3+Nacos1.1.3(配置中心)
    整理一些大厂的开源平台及github,向他们看齐...
    Spring Boot 2.x 基础案例:整合Dubbo 2.7.3+Nacos1.1.3(最新版)
    Java图片处理:ico格式转 PNG/JPG等格式
  • 原文地址:https://www.cnblogs.com/stelayuri/p/13503259.html
Copyright © 2011-2022 走看看