1、用高斯消元法解线性方程组北京景山学校 何江舟.GPA排名系统(CTSC2001)高等院校往往采用GPA来评价学生的学术表现。传统的排名方式是求每一个学生的平均成绩,以平均成绩作为依据进行排名。对于不同的课程,选课学生的平均成绩会受到课程的难易程度等因素的影响,因此这种排名方式不够合理。为此,我们需要对排名系统进行这样的改进:对第i门课的每一个学生的成绩加上一个特定的修正值di(调整后的成绩不按照百分制),使得经过调整后,该课的平均分等于选该课的所有学生的所有课的平均分。对每一门课都这样调整,使得上述条件对所有课程都满足。你的任务是根据一个年级学生某学年的成绩,通过上述调整,得出他们的排名。.简
2、要分析Ai:选修第i门课的学生的集合Bj:第j个学生选修课程的集合Gi,j:第j个学生第I门课的成绩di:第i门课的修正值对于第p门课,可列出如下关系式:pjppAjjpAjBiii,jAjBpp,jA)d(GdG11 这是关于di(i=1,2,n)的线性方程,我们可以整理出n个这样的方程。.线性方程组的一般形式a1,1x1+a1,2x2+a1,nxn=b1a2,1x1+a2,2x2+a2,nxn=b2an,1x1+an,2x2+an,nxn=bn下面是n元线性方程组的一般形式:我们可以把它表示为增广矩阵的形式:a1,1a1,2 a1,nb1a2,1a2,2 a2,nb2an,1an,2 an
3、,nbn.先看一个例子2-131425412072-1314-122.5-1.5 6.52-1314-12 -0.875 5.252 0.52.5得出:x3=5.25/(-0.875)=-6x2=(2-(-1)x3)/4=-1x1=(1-(-1)x2-3x3)/2=9.消元过程a1,1(1)a1,2(1)a1,n(1)b1(1)a2,1(1)a2,2(1)a2,n(1)b2(1)an,1(1)an,2(1)an,n(1)bn(1)注:用上标(k)表示第k次消元前的状态第1次消元,第1行的乘数:(i=2,3,n))1(1,1)1(1,1,aaiima1,1(1)a1,2(1)a1,n(1)b1(
4、1)a2,2(2)a2,n(2)b2(2)an,2(2)an,n(2)bn(2)得到新的增广矩阵:ai,j(2)=ai,j(1)-mi,1a1,j(1)bi(2)=bi(1)-mi,1b1(1)(i,j=2,3,n).第k次消元,第k行的乘数:(i=k+1,k+2,n)消元过程a1,1(1)a1,2(1)a1,n(1)b1(1)a2,2(2)a2,n(2)b2(2)ak,k(k)ak,n(k)bk(k)an,k(k)an,n(k)bn(k)第k次消元前的增广矩阵:)(,)(,kkkkkiaakimai,j(k+1)=ai,j(k)-mi,kak,j(k)bi(k+1)=bi(k)-mi,kbk
5、(k)增广矩阵的变化:(i,j=k+1,k+2,n)第k步消元的主行第k步消元的主元素.回代过程a1,1(1)a1,2(1)a1,n(1)b1(1)a2,2(2)a2,n(2)b2(2)an,n(n)bn(n)最后得到的增广矩阵:)(,1)(,)(iiinijjijiiiiaxabx最终结果的计算:.为什么要选主元素 前面介绍的消元法都是按照自然顺序,即x1、x2、xn的顺序消元的。有:)(,)(,kkkkkiaakim 所以每一次消元的主元素都不能为0。如果按照自然顺序消元的过程中出现的ak,k(k)=0,那么消元无法继续进行下去。或者|ak,k(k)|很小,也会严重影响计算精度。.为什么要
6、选主元素例如(假设运算过程中使用单精度实数):10-101111210-1011-1010-1010解得:x1=0,x2=1 这个解与第二个方程差异很大。究其原因,因为消元过程中第一个方程所乘的系数过大,使得上式“吃掉”了下式,所以在结果中根本无法体现下式。但如果调整一下顺序:11210-101111211解得:x1=1,x2=1,这个解基本符合原方程 所以每次消元的主元素的绝对值应该尽可能大,使得与主行相乘的乘数尽可能小。.选主元素a1,1(1)a1,2(1)a1,n(1)b1(1)a2,2(2)a2,n(2)b2(2)ak,k(k)ak,n(k)bk(k)al,k(k)al,n(k)bl(
7、k)an,k(k)an,n(k)bn(k)进行第k次消元时,将ak,k一下各元素(包括ak,k)进行比较,将其中的最大者所在行与第k行交换。.无解的情况 如果在消元的过程中,增广矩阵出现这样一行:左侧各未知数的系数都为0,而右侧的常数项不为0,则意味着方程组无解。.无数组解的情况 在消元过程中,出现这样一行:各未知数的系数和常数项都为0。这相当于少了一个方程,也就是接下来的消元过程中,方程的个数少于未知数的个数,方程要么无解,要么有无数组解。下面讨论对于这样的方程,如何得到一组解。先看这样一个方程:423921142125423900-0.5-0.5000.50.5 如果继续消元(消第2列),
8、必须保证a2,20,可是第2列中不存在非0的项。.无数组解的情况423900-0.5-0.5000.50.5 只能够把第3列的元素作为第2次消元的主元素,进行消元:423900-0.5-0.50000 第2次消元得到的元素全部为0,所以第三行元素已失去意义。x2没有做过主元素,可随意取值,再进行回代,得到一组可行解。如令x2=0,x3=1,x1=1.5。对于一般的线性方程组,先进行消元,每次消元前找到系数不完全为0的列,相应的元素作为此次消元的主元素,直至第k次消元时,得到的新元素全部为0,这时把各未知数分为两种:第k+1列至第n列对应的未知数,可以将这些未知数随意取值;第1列至第k列对应的未
9、知数,这些未知数的值在回代过程中确定。.性能分析时间复杂度:O(n3)消元O(n3)选主元素:O(n2)回代O(n2)空间复杂度:O(n2)增广矩阵O(n2)如使用全选主元素,还需一个存储列与元素对应信息的表,为O(n)精度:由于采用实数运算,另外每一次(第一次除外)消元都要使用以前消元产生的结果,每一次回代都要使用消元结果和其它回代结果,所以累积误差比较严重,该方法只能够求得近似解。但是可以根据具体需要进行相应改进。.整数线性方程组的精确解法 前面讨论了对于一般线性方程组通过实数运算得到近似解的算法。而在一些问题中,常常要求精确解,这里讨论一下系数、常数项和解均为整数的线性方程组的精确解法。
10、前面是用这种方法消元的:)(,)(,kkkkkiaakimai,j(k+1)=ai,j(k)-mi,kak,j(k)bi(k+1)=bi(k)-mi,kbk(k)显然这里进行的是实数运算。.整数线性方程组的精确解法 由于不能够保证ai,k(k)是ak,k(k)的倍数,要想消元,必须使两行分别乘以一个乘数。)(,)(,)(,)(,)(,)(,kkikkkkkikkkkkkkkiaaakiaaakimmai,j(k+1)=mi,kai,j(k)-mi,kak,j(k)bi(k+1)=mi,kbi(k)-mi,kbk(k)方程较多时,系数有可能越来越大,到一定程度有可能导致系数越界,因此要随时对各行
11、化简,即把这一行中所有元素除以这些元素的最大公约数。但是,无论如何,这也保证不了不会发生越界,因此这种算法一般适用于系数、未知数范围较小,未知数个数较少的方程。.齿轮 你有一套玩具,包括许多不同尺寸的齿轮(至多20种,假定每一种齿轮有无限多个),每个齿轮最多100齿。你希望用它们构造不同比例的传动装置。一个传动装置包括偶数个齿轮,这些齿轮两两一组互相咬合,每一组齿轮都与下一组用轴承相连。用c1、c2、cm表示每组第一个齿轮的齿数,用d1、d2、dm表示每组第二个齿轮的齿数。例如你有3种齿轮:6齿、12齿、30齿,你需要实现4:5的传动比例,一种可行的方案是:使用4个齿轮,分2组,第1组的两个分
12、别为12齿、6齿,第2组的两个分别为12齿、30齿。miiiDCdc1:.简要分析 把这些齿轮的齿数设为a1、a2、an,设它们作为C类齿轮的数量分别为e1、e2、en,作为D类齿轮的数量分别为f1、f2、fn。有如下关系:nifeinifieiiiiiaaadc11:这时候我们不难发现,一种齿轮同时当作C类、D类使用是一种浪费。设xi=ei-fi,xi0表示这种齿轮只作为C类,xi0表示这种齿轮只作为D类。这就转化为解xi问题。我们可以将c、d、ai这些值分解质因数。由于ai不超过100,所以a1an能够分解为的质因数不超过25种。另外,如果c或d中包括这以外的质因数,显然问题无解。.简要分
13、析 设gr,i为质数r在ai的质因数分解中的指数,cr、dr分别为质数r在c、d中的质因数分解中的指数。有如下关系:2(x1g2,1+x2g2,2+xng2,n)=2(c2-d2)3(x1g3,1+x2g3,2+xng3,n)=3(c3-d3)这完全可以表示为关于指数的等式,即:g2,1x1+g2,2x2+g2,nxn=c2-d2 g3,1x1+g3,2x2+g3,nxn=c3-d3 g97,1x1+g97,2x2+g97,nxn=c97-d97 当然还有一个约束条件:x1+x2+xn=0 这就完全转化为了解线性方程组的问题,而且这需要精确地解出这个整数线性方程组,并且还要面临不定方程的问题,如果能够得出整数解,则问题有解。.小结 高斯消元法是一种比较简单、适用范围较广的有效算法,但在实际应用中,我们往往需要具体问题具体分析,对这样的标准算法进行改进,才能满足我们的需要。.谢谢请多提宝贵意见.