《[理学]第二章 解线性方程组的直接法》由会员分享，可在线阅读，更多相关《[理学]第二章 解线性方程组的直接法（89页珍藏版）》请在金锄头文库上搜索。
1、第二章 解线性方程组的直接法,第二章 解线性方程组的直接法,引言 Gauss消元法 列主元素消元法 矩阵三角分解法 向量和矩阵的范数 误差分析,2.1 引言,小行星轨道问题：,天文学家要确定一小行星的轨道，在轨道平面建立以太阳为原点的直角坐标系。在坐标轴上取天文测量单位(一天文单位为地球到太阳的平均距离：9300万英里，约1.5亿千米)，对小行星作5次观察，测得轨道上5个点的坐标数据如下：,x 5.7640 6.2860 6.7590 7.1680 7.4800 y 0.6480 1.2020 1.8230 2.5260 3.3600,椭圆的一般方程: a1x2 + a2xy + a3y2 +
2、 a4x + a5y + 1 = 0,将数据逐个代入，可得五个方程的方程组，求解该线性方程组即可得行星轨道方程。,对一般线性方程组: A x = b, 其中,由以前所学内容知，当且仅当矩阵A行列式不为0时，即A非奇异时，方程组存在唯一解，可根据Cramer法则求解。,其算法设计如下：,（1） 输入系数矩阵A和右端向量b；,（2）计算系数矩阵A的行列式值D，如果D=0，则输出错误信息，结束，否则进行第(3)步；,（3） 对k=1,2,n，用b替换A的第k列数据，并计算替换后矩阵的行列式值Dk；,（4） 计算并输出x1 = D1 / D，x2 = D2 / D，xn=Dn/D， 结束。,但Cram
3、er法则只适用于低阶方程组，高阶方程组工作量太大，故一般用数值方法求解。数值方法分两类： 1. 直接法 2. 迭代法,2.2 Gauss消元法,第二步: 回代过程， 解上三角形方程组，得原方程组的解。,基本思想：逐步消去未知元，将方程组化为与其等价的上三角方程组求解。,分两步：,第一步: 消元过程，将方程组消元化为等价的上三角形方程组;,Gauss消元的目的：,原始方程组,约化方程组,消元过程(化一般方程组为上三角方程组),以四阶为例：,其系数增广矩阵为：,第一轮消元：,计算3个数: m21 m31 m41T = a21 a31 a41T / a11,用-m21乘矩阵第一行后加到矩阵第二行;,
4、用-m31乘矩阵第一行后加到矩阵第三行;,用-m41乘矩阵第一行后加到矩阵第四行;,其系数增广矩阵变为：,第二轮消元：,计算2个数：m32 m42T = a32(1) a42(1)T / a22(1),用-m32乘矩阵第二行后加到矩阵第三行;,用-m42乘矩阵第二行后加到矩阵第四行;,其系数增广矩阵变为：,第三轮消元：,计算: m43=a43(2)/a33(2),用-m43乘矩阵第三行后加到矩阵第四行;,其系数增广矩阵变为：,其对应的上三角方程组为,若对于一般的线性方程组Ax=b，其消元过程的计算公式为： （k=1,2,n-1）,回代过程(解上三角方程组),上三角方程组的一般形式为： 其中a1
5、1ann0,回代过程的计算公式：,工作量计算：,消去过程：,“”：第k步，n-k次，共 (n-1)+(n-2)+1=n(n-1)/2,“”：第k步，(n-k)(n-k+1)次，共 (n-1)n+(n-2)(n-1)+12= (n3-n)/3,总工作量： S1=n(n-1)/2+ (n3-n)/3,回代过程：,“”：n,“”：1+2+(n-1)=n(n-1)/2,总工作量：S2=n+ n(n-1)/2= n(n+1)/2,（k=1,2,n-1）,故Gauss消元法的总工作量为：,S=S1+S2 =n(n-1)/2+ (n3-n)/3+ n(n+1)/2 = n2+(n3-n)/3,若用克莱姆法则
6、求解，则工作量为：,“”：（n+1个n阶行列式的值）(n+1)(n-1)n!,“”：n,故总工作量为： (n+1)(n-1) n!+n,如当n=6时, Gauss消元法工作量为106 ；而克莱姆法则求解工作量为25206。,定理： 约化的主元素ak+1,k+1(k) 0 (k=0,1,n-1)的充分必要条件是矩阵A的各阶顺序主子式不为零。即,注：对角线上的元素ak+1,k+1(k)在Gauss消元法中作用突出，称约化的主元素。,推论: 如果A的顺序主子式Dk 0 (k=1,n-1)，则Gauss消元法中的约化主元可以表示为,x1= -13, x2 = 8, x3 = 2,m21=3/2 m31
7、=4/2,m32= -3/0.5,例 用高斯消元法求解方程组,矩阵的三角分解:,对线性方程组Ax=b的系数矩阵A施行初等行变换相当于用初等矩阵左乘A，故第一次消元后方程组化为A(1)x=b(1)，即L1Ax=A(1)x，L1b=b(1)，其中,同理 LkA(k-1)=A(k) Lkb(k-1)=b (k),其中,将A 分解为单位下三角矩阵 L 和上三角矩阵 U的乘积的算法称为矩阵A的三角分解算法。,重复该过程，最后得,记U=A(n-1)，则,其中,定理：设A为n阶矩阵，若A的顺序主子式Di 0（i=1,2,n-1），则A可分解为一个单位下三角矩阵L和一个上三角矩阵U的乘积，且这种分解是唯一的。
8、,由Gauss消元过程可推得,U即为Gauss消元后所得的上三角方程的系数矩阵。,例,解,由Gauss消元法可得，,m21=0， m31=2， m32= -1,故,如果已经有 A = L U 则 AX = b = L U X = b，,（1）求解方程组 LY = b 得向量 Y 的值； （L 是下三角矩阵，用顺代算法）,（2）求解方程组 UX = Y 得向量 X 的值。 （U 是上三角矩阵，用回代算法）,记UX = Y ， LY = b ，则求解方程组分两步进行：,基本思想：Gauss消元法中，若主元 akk(k) 太小会使误差增大，故应避免采用绝对值小的元素作主元。最好每一步选取系数矩阵中（
9、或消元后的低阶矩阵中）绝对值最大的元素作主元，以具较好的数值稳定性。,2.3 Gauss列主元素消元法,例：求解方程组,（用四位浮点数计算，精确解舍入到4位有效数字为：x1*=-0.4904,x2*=-0.05104,x3*=0.3675）,解：方法一Gauss消元法,（A/b）=,其中， m21=-1.000/0.001=-1000 m31=-2.000/0.001=-2000 m32=4001/2004=1.997 解为x1=-0.4000，x2=-0.09980，x3=0.4000 (x1*=-0.4904，x2*=-0.05104，x3*=0.3675) 显然，此解并不准确。,方法二交
10、换行，避免绝对值小的主元作除数。,（A/b）=,其中， m21=0.5000 m31=-0.0005 m32=0.6300,解为x1=-0.4900，x2=-0.05113，x3=0.3678,(x1*=-0.4904，x2*=-0.05104，x3*=0.3675),与方法一相比，此解显然要精确得多。,设Axb的增广矩阵为,在A的第一列中选绝对值最大的元素作主元，设该元素所在行为第i1行，交换第一行与第i1行，进行第一次消元；再在第2n行的第二列中选绝对值最大的元素作主元，设该元素所在行为第i2行，交换第二行与第i2行，进行第二次消元，直到消元过程完成为止。,Gauss列主元素消元法的基本思
11、想：,例：用列主元法解,解：第一列中绝对值最大是10，取10为主元,第二列的后两个数中选出主元 2.5,列主元矩阵的三角分解：,解：交换行变换,例：对矩阵A做列主元三角分解：,则列主元的Gauss变换可记为,A(2)=F2P2F1P1A,记 U=A(2)=F2(P2F1P2)P2P1A （因P2P2=I）,P = P2P1,则有,对于一般的n阶矩阵的列主元三角分解，通常令,定理：（列主元素的三角分解定理）若A非奇异，则存在排列阵P使PALU，其中L为下三角阵，U为上三角阵。,矩阵分解关系为,全主元素消元法：,定义：,则称,为全主元素。,经过行列互换，使得 位于经交换行和列 后的等价方程组中的
12、位置，然后再实施消元。,全主元素消元法的基本思想：,若,注：全主元素消元法有可能改变未知数的顺序。,直接三角分解法：若将A分解为LU的积，则求Axb等价于求解两个三角形方程组： （1）Lyb，求y； （2）Uxy，求x。,2.4 矩阵三角分解法,（一） Doolittle分解法,（二 ） Crout分解法,（三） 对称正定矩阵的Cholesky分解法,（四） 三对角方程组的数值解法,设A非奇异，且ALU，L为单位下三角阵，U为上三角阵，即,可得直接三角分解法解Axb的计算公式：,一、Doolittle分解法：,u11=a11, , u1n=a1n,m21u11=a21, , mn1u11=an
13、1,对A的元素aij ，当 jk 和 ik时,m21u12+ u22=a22, , m21u1n+ u2n=a2n,m21=a21/u11, , mn1=an1/u11,u22=a22 - m21u12, , u2n=a2n- m21u1n,m31u12+ m32u22=a32, , mn1u12+ mn2u22=an2,m32=(a32- m31u12)/u22, , mn2=(an2- mn1u12)/u22,矩阵L和矩阵U的元素计算公式：,计算秩序如图所示：,用Doolittle分解法求解线性方程组Axb的具体计算公式如下：,例：用Doolittle法解方程组,解:,由Doolittle
14、分解,解Lyb，得,解Uxy，得,直接三角分解的Doolittle方法可用以下过程表示：,此方法称紧凑格式的Doolittle法。,从上式最后一个矩阵中可知L，U，y， 然后解线性方程组Uxy。,例：用紧凑格式的Doolittle分解法解上例方程组,解：,所以,计算公式：,二、Crout分解法：,将矩阵A分解为如下形式：,例： 求矩阵A的Crout分解：,所以,若n阶实矩阵A为对称正定矩阵，则： （1）ATA； （2）对任意的X0，有XTAX0； （3） A的各阶顺序主子式大于零。,三、对称正定矩阵的Cholesky分解法（平方根法）,故A可进行LU分解：,记,则 A=LDU1,即 ALDLT
15、,称A的LDLT分解,若记,则,称A的LLT分解,= AT = U1TDLT,= U1=LT,LLT分解的计算公式：,对A的第i行元素, 有 ( j = 1,2,i ),( j =1,2,i-1 ),对于 i = 2,3, n,LDLT分解的计算公式：,d1 = a11, l21= a21/d1, , ln1 = an1/d1,( i = 2,3, n； j = 1,2,i-1 ),无需开方，故称改进的平方根法,应用Cholesky分解可将方程组Axb分解为两个三角形方程组：,而应用改进的Cholesky分解可将方程组Axb分解为下面的方程组：,例：用改进的平方根法解方程组Axb，其中,解：由,d1 = a11, l21= a21/d1, , ln1 = an1/d1,d1=1，l21=2，l31=1，l41=-3,得,又由,得,d2=1，l32=-2，l42=1， d3=9，l43=2/3，d4=1,因此，得到,解方程组Lyb，得,解方程组LTxD-1y，得,最终求得方程组Axb的解为 x=(1,1,1,1)T,四、 三对角方程组的数值解法,且 ai ci 0， (i = 2，3，n - 1),其中,三对角方程组形式如下：,计算公式：,Step
