高斯消元法,是线性代数中的一个算法,可用来求解线性方程组,并可以求出矩阵的秩,以及求出可逆方阵的逆矩阵。
高斯消元法的原理是:
若用初等行变换将增广矩阵 化为 ,则AX = B与CX = D是同解方程组。
所以我们可以用初等行变换把增广矩阵转换为行阶梯阵,然后回代求出方程的解。
1、线性方程组
1)构造增广矩阵,即系数矩阵A增加上常数向量b(A|b)
2)通过以交换行、某行乘以非负常数和两行相加这三种初等变化将原系统转化为更简单的三角形式(triangular form)
注:这里的初等变化可以通过系数矩阵A乘上初等矩阵E来实现
3)从而得到简化的三角方阵组,注意它更容易解
z=-4/-4=1, y=4-2z=4-2=2, x= (1-y-z)/2=(1-2-1)/2=-1
总结上面过程,高斯消元法其实就是下面非常简单的过程
原线性方程组 ——> 高斯消元法 ——> 下三角或上三角形式的线性方程组 ——> 前向替换算法求解(对于上三角形式,采用后向替换算法)
补充1:
高斯-若尔当消元法(Gauss-Jordan Elimination)
相对于高斯消元法,高斯-若尔当消元法最后的得到线性方程组更容易求解,它得到的是简化行列式。其转化后的增高矩阵形式如下,因此它可以直接求出方程的解,而无需使用替换算法。但是,此算法的效率较低。
例子如下:
解为
个人感觉区别就是对每行进行了归一化处理
补充2:
介绍了最基本的高斯消元法,现在看看应用于实际问题的实用算法
因为实际应用中,我们总是利用计算机来分析线性系统,而计算机中以有限的数来近似无限的实数,因此产生舍入误差(roundoff error),进而对解线性系统产生很多影响。
一个t位(即精度为t)以为基的浮点数的表达形式为:,。对于一个实数x,其浮点近似值为最接近x的浮点数,必要时进行近似。
例1:对2位以10为基的浮点算法,。
例2:同样考虑,。
以下面系统为例,看看在高斯消元中采用浮点算法会产生什么效果。
当以精确解法时,通过将第一行乘以m=89/47,并从第二行中减去得到,进而利用后向替换算法得x=1,y=-1。
当以3位以10为基的浮点算法时,乘子变为,因为,因此第一步高斯消元后得
。此时,因为不能将第2行第1列位置变为0,所以不能将其三角化。从而,我们只能接受将这个位置值赋为0,而不管其实际浮点值。因此,3位浮点高斯消元的结果为,后向算法计算结果为。
尽管无法消除近似误差的影响,可以采用一些技术来尽量减小这类机器误差。部分主元消元法在高斯消元的每一步,都选择列上最大值为轴(通过行变换将其移动)。
下面给出列主元消去的代码(所谓列主元消去法是在系数矩阵中按列选取元素绝对值得最大值作为主元素,然后交换所在行与主元素所在行的位置,再按顺序消去法进行消元。)
1 列选主元消元法2 /*
3 *Gauss.cpp4 *功能:列选主元消元法5 */
6 #include
7 #include"Gass.h"
8
9 //高斯消元法(列选主元)
10 void Gauss (double a[][MAXNUM],intn)11 {12 inti,j;13 SelectColE(a,n); //列选主元并消元成上三角14 //回代求解
15 printf("上三角的结果\n");16 printM(a,3);17 for(i=n;i>=1;i--)18 {19 for(j=i+1;j<=n;j++)20 a[i][n+1]-=a[i][j]*a[j][n+1];21 a[i][n+1]/=a[i][i];22 }23 return;24 }25 //选择列主元并进行消元
26 void SelectColE(double a[][MAXNUM],intn)27 {28 inti,j,k,maxRowE;29 double temp; //用于记录消元时的因数
30 for(j=1;j<=n;j++)31 {32 maxRowE=j;33 for(i=j;i<=n;i++)34 if(fabs(a[i][j])>fabs(a[maxRowE][j]))35 maxRowE =i;36 if(maxRowE!=j)37 swapRow(a,j,maxRowE,n); //与最大主元所在行交换38 //消元
39 for(i=j+1;i<=n;i++)40 {41 temp =a[i][j]/a[j][j];42 for(k=j;k<=n+1;k++)43 a[i][k]-=a[j][k]*temp;44 }45 printf("第%d列消元后:\n",j);46 printM(a,3);47 }48 return;49 }50 void swapRow(double a[][MAXNUM],int m,int maxRowE,intn)51 {52 intk;53 doubletemp;54 for(k=m;k<=n+1;k++)55 {56 temp =a[m][k];57 a[m][k] =a[maxRowE][k];58 a[maxRowE][k] =temp;59 }60 return;61 }62
63 //测试函数
64 intmain()65 {66 double a[4][MAXNUM];67
68 inti,n,j;69
70 a[1][1] = 0.001; a[1][2]=2.000;a[1][3]=3.000;a[1][4]=1.000;71 a[2][1] = -1.000;a[2][2]=3.712;a[2][3]=4.623;a[2][4]=2.000;72 a[3][1] = -2.000;a[3][2]=1.070;a[3][3]=5.643;a[3][4]=3.000;73 Gauss (a,3);74 for(i=1;i<=3;i++)75 printf("X%d=%f\n",i,a[i][4]);76 return 0;77 }78 //输出矩阵
79 void printM(double a[][MAXNUM], intn)80 {81 inti,j;82 for(i=1;i<=n;i++)83 {84 for(j=1;j<=n+1;j++)85 {86 printf("%f ,",a[i][j]);87 }88 printf("\n");89 }90 }
测试结果:
2、逆矩阵
下面来说说高斯消元法在编程中的应用。
首先,先介绍程序中高斯消元法的步骤:
(我们设方程组中方程的个数为equ,变元的个数为var,注意:一般情况下是n个方程,n个变元,但是有些题目就故意让方程数与变元数不同)
1. 把方程组转换成增广矩阵。
2. 利用初等行变换来把增广矩阵转换成行阶梯阵。
枚举k从0到equ – 1,当前处理的列为col(初始为0) ,每次找第k行以下(包括第k行),col列中元素绝对值最大的列与第k行交换。如果col列中的元素全为0,那么则处理col + 1列,k不变。
3. 转换为行阶梯阵,判断解的情况。
① 无解
当方程中出现(0, 0, …, 0, a)的形式,且a != 0时,说明是无解的。
② 唯一解
条件是k = equ,即行阶梯阵形成了严格的上三角阵。利用回代逐一求出解集。
③ 无穷解。
条件是k < equ,即不能形成严格的上三角形,自由变元的个数即为equ – k,但有些题目要求判断哪些变元是不缺定的。
这里单独介绍下这种解法:
首先,自由变元有var - k个,即不确定的变元至少有var - k个。我们先把所有的变元视为不确定的。在每个方程中判断不确定变元的个数,如果大于1个,则该方程无法求解。如果只有1个变元