zoukankan      html  css  js  c++  java
  • 数论小白都能看懂的线性方程组及其解法(高斯消元)

    此文章依 CC 4.0 BY-SA 版权协议转载自 ShineEternal 的博客

    -1. 序言

    说到线性方程组,大家第一反应大概就是高斯消元,本文将对其详细讲解并配合例题与相关方法为您呈现。

    本文因图文并茂有较多配图且讲解详细较多,再过多的放置代码会引起文章的冗长以及阅读的不适,故只将模板放在上面。

    0、由来

    遇到形如:

    [left{egin{array}{l}{ a_{11} x_{1}+a_{12} x_{2}+cdots+a_{1 n} x_{n}=b_{1}} \ {a_{21} x_{1}+a_{22} x_{2}+cdots+a_{2 n} x_{n}=b_{2}} \ {cdots} \ {cdots} \ {cdots} \ {a_{n 1} x_{1}+a_{n 2} x_{2}+cdots+a_{n n} x_{n}=b_{n}} end{array} ight. ]

    这样的方程组,你会怎么解呢?

    如果是在数学领域中,可能这样的方程组中方程的个数是屈指可数的,这样就特别简单,一通乱搞就出来了。

    但是如果在 OI 领域遇到这类问题呢?我们发现,计算机中必须由固定的算法,才能实现某一问题,而人脑能进行综合性的思考,但是却无法将大脑中的神经元活动表述到程序中,这是人脑与电脑的区别,也就由此诞生了各种解线性方程组的方法。

    说到这里,我想起了一本著名的书:计算机与人脑。冯诺依曼写的,有兴趣者可以自行购买阅读也能普及一些历史

    这种线性的方程组在工程中也运用广泛,这种模型中:

    • 未知量较多
    • 方程个数不同

    而在考虑解方程时,要考虑:

    • 是否有解?
    • 如果有解,是否唯一?
    • 如果不唯一,解有什么规律与结构?

    由此,线性方程组的解决方法就显得尤为重要。

    下面就来介绍一下两种常用的方法及相关例题

    1、高斯消元

    插播一些历史知识:

    该方法以数学家高斯命名,由拉布扎比。伊丁特改进,发表于法国但最早出现于中国古籍《九章算术》,成书于约公元前 (150) 年。

    高斯:

    德国著名数学家、物理学家、天文学家、几何学家,大地测量学家。享有「数学王子」的美誉,真的是多才多艺,著名的「自然数 (1sim100) 的求和问题」就是他 (9) 岁时计算出来的,虽然可能现在对我们来说轻而易举,但是的确体现了他那时的创新精神,为日后的研究打下了坚实的基础。

    这种是大家比较耳熟能详的解法了。

    那么,这种解法的核心思想是什么呢?

    我们先来构想一下,什么样子的线性方程式我们用电脑也能写出解法?

    显然,我们如果从某一个式子开始,每一个式子都能得到一个未知数的解,而得到了这个未知数的解,又能代入另一个方程,再得到一个未知数的解。

    以此类推。

    总结一下,就是一个「三角形分割」(自己瞎起的名字)

    还是拿上面那个普遍形式:

    [left{egin{array}{l}{ a_{11} x_{1}+a_{12} x_{2}+cdots+a_{1 n} x_{n}=b_{1}} \ {a_{21} x_{1}+a_{22} x_{2}+cdots+a_{2 n} x_{n}=b_{2}} \ {cdots} \ {cdots} \ {cdots} \ {a_{n 1} x_{1}+a_{n 2} x_{2}+cdots+a_{n n} x_{n}=b_{n}} end{array} ight. ]

    但是为了观察方便,我们假设一个三行四列的式子:(真香)

    [left{egin{array}{l}{x_{1}+3 x_{2}+4 x_{3}=5} \ {x_{1}+4 x_{2}+7 x_{3}=3} \ {9 x_{1}+3 x_{2}+2 x_{3}=2}end{array} ight. ]

    此时把系数写到矩阵里:(没学过矩阵也没关系,就当成一个表示就行啦)

    [left[egin{array}{llll}{1} & {3} & {4} & {5} \ {1} & {4} & {7} & {3} \ {9} & {3} & {2} & {2}end{array} ight] ]

    其实就是 高斯消元模板题的样例

    第一步:对左侧的系数进行消元

    [left[egin{array}{lll} {color{red}1} & {3} & {4} \ {color{green}1} & {color{blue}4} & {7} \ {color{green}9} & {color{orange}3} & {color{purple}2}end{array} ight] ]

    进行向「三角形分割」的转换

    即目标状态为

    [left[egin{array}{lll}{?} & {?} & {?} \ {0} & {?} & {?} \ {0} & {0} & {?}end{array} ight] ]

    此时,我们设标红的 (color{red}1) 为主元,对照目标状态可知,标绿的 (color{green}1)(color{green}9) 需要化 (0)

    所以我们只需要将第二行整体减去第一行,就实现了将绿 (color{green}1)(0)

    那么 (9) 如何转化呢?

    此时就对应的整体减去 (9) 倍的第一行就行(可能会有更简便的系数化 (0) 法,但是我们要给计算机设计一个统一的算法才能实现)

    于是就得到了:

    [left[egin{array}{ccc} {1} & {3} & {4} \ {0} & {color{blue}1} & {3} \ {0} & {color{orange}{-24}} & {-34}end{array} ight] ]

    然后再以蓝 (color{blue}1) 为主元,现在我们需要将 (color{orange}{-24}) 化为 0

    于是将第三行整体减去 (−24) 倍第二行(其实相当于加 (24) 倍)

    Q:这里不会把之前第一列化好的 (0) 给冲掉吗?
    A: 显然,在上一步已经保证第一列除主元外系数皆为 (0),所以不会出现上述情况

    Q: 万一上面那个绿 (color{green}1)(0) 怎么办?
    A: 这就是比较特殊的主元为 (0) 的情况,此时因为每一个方程地位都相等,所以只需将主元所在的行与下方主元位置非 (0) 的行互换即可。

    追加 Q:万一下面也没有呢?
    A: 那么方程就无解了,这种情况在后面也会提到。

    于是就和

    [left[egin{array}{lll}{?} & {?} & {?} \ {0} & {?} & {?} \ {0} & {0} & {?}end{array} ight] ]

    这个目标状态对应起来了。

    然后就想摩拳擦掌的解了。

    但是发现右边的 (b) 还没有代入矩阵。

    于是

    第二步:对增广矩阵消元

    所谓增广矩阵?

    增广矩阵(又称扩增矩阵)就是在系数矩阵的右边添上一列,这一列是线性方程组的等号右边的值。 -----百度百科

    其实就是多出 (b) 那一列啦

    [left[egin{array}{llll}{1} & {3} & {4} & {5} \ {1} & {4} & {7} & {3} \ {9} & {3} & {2} & {2}end{array} ight] ]

    如上

    增广矩阵的操作与左侧的系数矩阵完全一样,刚开始单列系数矩阵只是一个由简入难的过程,加深对高斯消元的理解。限于篇幅原因,就不再赘述。

    第三步:回带

    最后经过一系列的消元,就得到了如下矩阵:

    [left[egin{array}{cccc}{1} & {3} & {4} & {5} \ {0} & {1} & {3} & {-2} \ {0} & {0} & {38} & {-91}end{array} ight] ]

    然后把矩阵再写回方程组:

    [left{egin{array}{l}{x_{1}+3 x_{2}+4 x_{3}=5} \ {0 x_{1}+1 x_{2}+3 x_{3}=-2} \ {0 x_{1}+0 x_{2}+38 x_{3}=-91}end{array} ight. ]

    这样我们显而易见地就可以去掉系数为 (0) 的字母:

    [left{egin{array}{l}{x_{1}+3 x_{2}+4 x_{3}=5} \ {1 x_{2}+3 x_{3}=-2} \ {38 x_{3}=-91}end{array} ight. ]

    现在结果明了了,从下往上推:

    依次可以得到

    [x_3​=frac{-91}{38}​=−2.39 ]

    (保留了两位小数)

    然后已知 (x_3),将其代入中间式子中的 (x_3),同理解得

    [x_2=5.18 ]

    最后将 (x_2​,x_3)​ 都代回最上面的式子,得:

    [x_1=−0.97 ]

    于是,对高斯消元的手动模拟就结束了

    显然,这个在数学考试中会答不完题的,但是对于计算机来说,却是一个通用的方法。

    时间复杂度:(O(n^3))

    虽然时间复杂度较高,不管怎样,我们已经找到了解决线性方程组的计算机通用方法,下面就用代码实现:

    #include<cstdio>
    #include<cmath>
    #include<iostream>
    using namespace std;
    const double eps=1e-7 ;
    double a[105][105],ans[105];
    int main()
    {
        int n;
    	scanf("%d",&n);
        for(int i=1;i<=n;i++)
        { 
            for(int j=1;j<=n+1;j++)//这里千万别忘了右侧那一列了 
            {
                scanf("%lf",&a[i][j]);
            }
        }
        for(int i=1;i<=n;i++)
    	{
            int m=i;
            for(int j=i+1;j<=n;j++)
                if(fabs(a[m][i])<fabs(a[j][i]))//fabs 能取浮点数绝对值
                    m=j;//找到当前这一列最大的数字,作为主元消元,这样能最大限度的避免精度误差 
            if(fabs(a[m][i])<eps)//这里要判精度,其实就是如果该位置的系数为 $0$ 则无解(double 没法准确处理这种精度 
    		{
                printf("No Solution
    ");
                return 0;
            }
            if(i!=m)swap(a[i],a[m]);//如果巧了,当前行不是最大行,那得罪了,你换下去吧,让未来的主元上来,好直接枚举之后的方程就行了 
            double d=a[i][i];
            for(int j=i;j<=n+1;j++)
                a[i][j]/=d;//将该位置的系数变为 1 
            for(int j=i+1;j<=n;j++)
    		{
                d=a[j][i];
                for(int k=i;k<=n+1;k++)
                {
                    a[j][k]-=a[i][k]*d;//将其他的方程用两式相减的方法减去应当减掉的系数的值 
                }
            }
        }
        ans[n]=a[n][n+1];
        for(int i=n-1;i>=1;i--)//回带,最后一行直接写出答案,然后其他行还有等待前面处理出来,从后往前推即可 
    	{
            ans[i]=a[i][n+1];//这里不难理解 
            for(int j=i+1;j<=n;j++)
            {
                ans[i]-=a[i][j]*ans[j]; 
        	}
        }
        for(int i=1;i<=n;i++)
        { 
            printf("%.2lf
    ",ans[i]);
        }
        return 0;
    }
    

    这个模板可以在 P3389 【模板】高斯消元法 提交测试(话说貌似我之前已经放过一遍链接了

    其他的地方应该都不难理解,就是为什么要让主元是最大的来避免误差?

    这个是从期望上来证明的,有兴趣者阔以看看这一篇 日报,应该就能想出来了吧。

    大概就是:

    设当前要消元的式子系数为 (q_{i1},q_{i2},q_{i3},cdots,q_{in})

    并假定我们主要针对的系数是 (q_{i1})

    则有

    [q_{j n}=q_{j n}-frac{q_{j 1}}{q_{i 1}} imes q_{i w} ]

    此时我们求的最大主元就是 (q_{i1})
    ,它越大,就说明 (dfrac{q_{j1}}{q_{i1}}) ​​的期望越小,这时较小的数给整个结果带来的影响也小,那么就不容易造成精度误差了 QAQ

    到此,高斯消元就结束了,但是

    还有高斯消元的进阶版,能避免回带来求出答案。

    2、高斯——约旦消元

    说到避免回带来求出答案,大家应该能想到这种算法的目的,还是拿样例来说,就是把方程组变成:

    [left[egin{array}{llll} {?} & {0} & {0} & {?} \ {0} & {?} & {0} & {?} \ {0} & {0} & {?} & {?} end{array} ight] ]

    这样,每个方程就能独立的解了。

    那么,我们如何得到这种形式呢?

    首先,我们依照朴素的高斯消元不难得到:

    [left[egin{array}{llll} {?} & {?} & {?} & {?} \ {0} & {?} & {?} & {?} \ {0} & {0} & {?} & {?} end{array} ight] ]

    观察上下两个矩阵,不难得到,我们在消元时不仅仅消去主元所在式子下方的式子,而对于上方的式子也应当予以处理。

    于是就开始了,下方是原图:

    [left[egin{array}{llll}{1} & {3} & {4} & {5} \ {1} & {4} & {7} & {3} \ {9} & {3} & {2} & {2}end{array} ight] ]

    第一次与普通消元类似,即得:

    [left[egin{array}{rrrr} {1} & {3} & {4} & {5} \ {0} & {1} & {3} & {-2} \ {0} & {-24} & {-34} & {-43} end{array} ight] ]

    但是第二次也要将第一行进行消元,得到:

    [left[egin{array}{rrrr} {1} & {0} & {-5} & {11} \ {0} & {1} & {3} & {-2} \ {0} & {0} & {38} & {-91} end{array} ight] ]

    然后就是最后一步(这一步是比朴素高斯消元要多的):

    [left[egin{array}{rrrr} {1} & {0} & {0} & {-0.97} \ {0} & {1} & {0} & {5.18} \ {0} & {0} & {38} & {-91} end{array} ight] ]

    把矩阵写回方程组

    [left{egin{array}{l}{x_{1}+0 x_{2}+0 x_{3}=-0.97} \ {0 x_{1}+x_{2}+0 x_{3}=5.18} \ {0 x_{1}+0 x_{2}+38 x_{3}=-91}end{array} ight. ]

    系数为 (0) 去掉:

    [left{egin{array}{l}{x_{1}=-0.97} \ {x_{2}=5.18} \ {38 x_{3}=-91}end{array} ight. ]

    这下子结果显然:

    [left{egin{array}{l}{x_{1}=-0.97} \ {x_{2}=5.18} \ {x_{3}=-2.39}end{array} ight. ]

    下面是模板代码:

    #include<cstdio>
    #include<cmath>
    #include<iostream>
    using namespace std;
    double a[105][105];
    int main()
    {
    	int n;
        scanf("%d",&n);
        for(int i=1;i<=n;i++)
        {
            for(int j=1;j<=n+1;j++)
            {
                scanf("%lf",&a[i][j]);
            }
        }
        for(int i=1;i<=n;i++)//枚举列(项) 
        {
            int m=i;
            for(int j=i+1;j<=n;j++)//选出该列最大系数 
            {
                if(fabs(a[j][i])>fabs(a[m][i]))//fabs 是取浮点数的绝对值的函数
                {
                    m=j;
                }
            }
            for(int j=1;j<=n+1;j++)//交换
            {
                swap(a[i][j],a[m][j]);
            }
            if(a[i][i]==0)//最大值等于 $0$ 则说明该列都为 $0$,肯定无解 
            {
                printf("No Solution
    ");
                return 0;
            }
            for(int j=1;j<=n;j++)//每一项都减去一个数
            {
                if(j!=i)//不是主元那一项 
                {
                    double d=a[j][i]/a[i][i];
                    for(int k=i+1;k<=n+1;k++)
                    {
                        a[j][k]-=a[i][k]*d;
                    }
                }
            }
        }
        for(int i=1;i<=n;i++)//最后的结果系数可能不为 $1$,所以记得消去常数 
        {
            printf("%.2lf
    ",a[i][n+1]/a[i][i]);
        }
        return 0;
    }
    

    显然,这一份代码较上一份简练一些

    3、例题

    例题 (1)

    P2455 [SDOI2006] 线性方程组

    题意简叙:

    已知 (n) 元线性一次方程组,判断解的情况。

    • 无实数解输出-1
    • 无穷多实数解输出 0
    • 有唯一解,则输出解(小数点后保留两位小数)。具体格式见样例

    分析:

    这道题只是细化了前面模板题对无法得出唯一解情况的判断,于是顺便讲一下无解&&多解时的判断方法:

    (我们知道,只要一个未知数无解或多解就可以说明整个方程的情况)

    无解

    我们最直接的就是判断最后一行是否是

    [left[egin{array}{llll} {?} & {?} & {?} & {?} \ {0} & {?} & {?} & {?} \ {0} & {0} & {?} & {?} end{array} ight] ]

    左边系数为 (0) 而右边系数不为 (0)

    但在容易疏漏情况,我们考虑到每一行,都会在其前面的消元中将唯一一个未知数留下来,也就都等同于第一行。

    于是,每一行都可以按照最后一行的套路模板来判断:

    if(a[i][i]==0&&a[i][n+1]!=0) 无解

    多解

    与上面对应的,我们很容易想到:

    if(a[i][i]==0&&a[i][n+1]==0) 多解

    • 要注意先判断无解,如果是的话直接跳出,然后轮到多解。

    因为多解首先建立在有解的基础上,如果哪一个未知数无解,整个方程自然不存在解,更不必说多解了。

    然后其他的就和模板题完全一样了

    例题 2

    (这题我没有找到出处)

    题目:【维他命的配方】

    题意

    现有若干种维他命,问能否利用这些维他命配制出适合人需求的各种维生素

    输入格式

    第一行:人需补充的维生素种类数 (V(1le Vle 25))

    第二行:(V) 个数,第 (i) 个数为 ( V_i)​,表示人体对第 (i) 种维生素的需求量 ((1le Vile 1000))

    第三行:已知的维他命种类 (G(1le Gle15))

    接下来是一个 (G imes V) 的整数矩阵,(A_{ij}) ​表示第 (i) 种维他命中所含的第 (j) 中维生素的含量 ((1le A_{ij}le 000))

    输出格式

    第一行输出能否配制,能则「YES」, 否则「NO」

    第二行:若能配制则输出 (G) 个整数,Gi

    Gi​表示第 (i) 中维他命所取的数列,多种方案输出任意一组。不能配制酒无需输出此行。

    输入样例

    4
    100 200 300 400
    4
    50 50 50 50
    30 100 100 100
    20 50 150 250
    50 100 150 200
    

    输出样例:

    YES
    1 1 1 0
    

    分析:

    首先要弄清维他命与维生素的区别。

    然后就应该不难想到高斯消元。

    设需要配制第 (i) 种维他命数量为 (x_i(ile 15))

    根据题意可列出:

    [left{egin{array}{l}{50 x_{1}+30 x_{2}+20 x_{3}+50 x_{4}=100} \ {50 x_{1}+100 x_{2}+50 x_{3}+100 x_{4}=200} \ {50 x_{1}+100 x_{2}+150 x_{3}+150 x_{4}=300} \ {50 x_{1}+100 x_{2}+250 x_{3}+200 x_{4}=400}end{array} ight. ]

    • 这第一步很好想到,但是行列的顺序却容易弄错,仔细观察就会发现需要把 (V_i) 补到矩阵的第 (G+1) 行,然后每一列为一个方程计算。

    列成矩阵的形式:

    [left[egin{array}{lllll}{50} & {30} & {20} & {50} & {100} \ {50} & {100} & {50} & {100} & {200} \ {50} & {100} & {150} & {150} & {300} \ {50} & {100} & {250} & {200} & {400}end{array} ight]]

    一通消元变成:

    [egin{array}{ccccc}{1} & {2} & {1} & {2} & {4} \ {0} & {1} & {3 / 7} & {5 / 7} & {10 / 7} \ {0} & {0} & {1} & {1 / 2} & {1} \ {0} & {0} & {0} & {0} & {0}end{array} ]

    这时,我们惊奇的发现:最后一行都是 (0)

    (呀这怎么办!)

    (刚才的知识怎么学的!不就是个多解情况吗?!)

    好了,多不多解不是关键,关键在于我们要怎么处理这一情况

    但这是关键,也并不难。

    我们只需懒一些将可能多解的未知数设为 (0),来尽可能的减少运算就行。

    具体针对这个样例来说,就是设 (x_4=0),他的数量与答案无关,并不造成影响。

    4、高斯消元注意事项:

    • 模板不要背错
    • 为了避免精度问题,最好和 eps 比较而不是和 (0) 比较(虽然有时这个基本相同)
    • 无解多解的判断情况以及顺序
    • 遇到题目转换成高斯消元时不要兴奋的转换出错或弄反

    5、后记

    这一篇高斯消元的文章终于完工了,里面的矩阵配图和线性方程组的 LaTeX 书写不易,如果这篇文章对您有帮助,那请不吝赐赞。

  • 相关阅读:
    C# .Net WinForm控件GDI+重绘位置错乱
    查询视图对应的基表名以及视图字段和对应的基表字段名
    解决在高分屏下开发winform界面变形
    ping命令工具:同时ping多个IP
    SmartAssembly批处理用法
    C#二维数组的初始化和存取
    win7 X64 进程名称不一致,导致杀进程失效!
    在存储过程中声明局部游标以循环调用自身
    强制设置双缓冲DoubleBuffered 解决tableLayoutPanel 闪烁
    Using SmartAssembly with MSBuild
  • 原文地址:https://www.cnblogs.com/1024th/p/11622317.html
Copyright © 2011-2022 走看看