线性代数方程组.ppt

第二章,线性方程组的数值解法,确定小行星轨道,以太阳为原点在轨道平面内建立直角坐标系,取天文测量单位,在五个不同时间观察小行星,测得坐标数据：,通过计算确定椭圆方程a1x2+2a2xy+a3 y2+2a4 x+2a5 y+1=0,a1xj2+2a2xjyj+a3 yj2+2a4 xj+2a5 yj+1=0,将五个点的坐标(xj,yj)(j=1,2,3,4,5)代入二次曲线方程,得关于a1,a2,a3,a4,a5 的方程组,在科学计算中，经常需要求解含有n个未知量 的n个方程构成的线性方程组,方程组还可以用矩阵形式表示为:Ax=b,(2.1),克莱姆法则需(n+1)(n-1)n!+n次乘除法运算(1)输入系数矩阵A和右端向量b;,(4)计算并输出x1=D1/D;,xn=Dn/D,结束。,(3)对k=1,2,n用b替换A的第k列数据,并计算替换后矩阵的行列式值Dk;,(2)计算系数矩阵A的行列式值D,如果D=0,则输出错误信息结束,否则进行第(3)步;,线性方程组数值解法的分类：,线性方程组数值解法的分类：,直接法 Gauss消去法及其变形 矩阵的三角分解法,迭代法 Jacobi迭代法 Gauss-Seidel迭代法 松弛迭代,2.1 高斯消元法,1三角形方程组的解法-回代法,2.顺序Gauss（高斯）消元法是一种规则化的加减消元法。,基 本 思 想,通过逐次消元计算把需要求解的线性方程组转化成上三角形方程组，也就是把线性方程组的系数矩阵转化为上三角矩阵，从而使一般线性方程组的求解转化为等价（同解）的上三角形方程组的求解。,Gauss消元法由消元和回代两个过程组成，先讨论一个具体的线性方程组的求解。,例1 用Gauss消元法解方程组,用增广矩阵进行进算,这样，对于方程组,(2.1),用增广矩阵表示，并给出Gauss消元法的具体步骤：,或者 Ax=b,顺序Gauss消元法的消元过程可表述如下：,第一步：设 a11(1)0，将第一列a11(1)以下各元素消成零,乘以矩阵A（1），b（1）的第一行再加到第i 行，得到矩阵,(i2，3，n),即依次用,其中,第二步：设 a22(2)0，将第二列a22(2)以下各元素消成零,（i3，4，n）,即依次用,乘以矩阵A（2），b（2）的第二行再加到第i行，得到矩阵,其中,如此继续消元下去第n1步结束后，得到矩阵,增广矩阵A（n），b（n）对应如下上三角形方程组,这是与原线性方程组（2.1）等价的方程组.,对于等价方程组,进行回代求解，可以得到：,首先写出增广矩阵,于是，采用Gauss消元法求解方程组,(2.1),然后进行消元，采用公式,最后进行回代得到方程组的解,得到相似增广矩阵,（ik+1，k+2，n）,在编程计算时，最后的增广矩阵存放的元素是：,算法.,顺序Gauss消元法可执行的前提,定理 1 给定线性方程组，如果n阶方阵 的所有顺序主子式都不为零，即 则按顺序Gauss消去法所形成的各主元素 均不为零，从而Gauss 消元法可顺利执行。,注：当线性方程组的系数矩阵为对称正定或严格对角占优阵时，按Gauss消元法计算是稳定的。,例2 用Gauss消元法求解方程组：,例3 用Gauss消元法求解线性方程组,a=1 2 1-2;2 5 3-2;-2-2 3 5;1 3 2 3;b=4 7-1 0;(ab)for k=1:3 d=a(k,k);c=a(k,:);c0=b(k);for i=k+1:4 l=a(i,k)/d;a(i,:)=a(i,:)-l*c;b(i)=b(i)-l*c0;endendb(4)=b(4)/a(4,4);for k=3:-1:1 b(k)=(b(k)-a(k,k+1:4)*b(k+1:4)/a(k,k);endb=b,MATLAB 程序,ans=2-1 2-1 b=2-1 2-1,3、列主元Gauss消元法,顺序Gauss消元法计算过程中的 akk(k)称为主元素，在第k步消元时要用它作除数,则可能会出现以下几种情况,若出现 akk(k)0，消元过程就不能进行下去。,akk(k)0，消元过程能够进行，但若|akk(k)|过小，也会造成舍入误差积累很大导致计算解的精度下降。,例4:在四位十进制的限制下，试用顺序Gauss消元法求解如下方程组,此方程组具有四位有效数字的精确解为,x117.46，x245.76，x35.546,解 用顺序Gauss消元法求解，消元过程如下,经回代求解得 x35.546，x2100.0，x1104.0,和此方程组的精确解相比,x35.546，x245.76，x117.46,有较大的误差。,对于此例，由于顺序Gauss消元法中的主元素绝对值非常小，使消元乘数绝对值非常大，计算过程中出现大数吃掉小数现象，产生较大的舍入误差，最终导致计算解 x1104.0 和 x2100.0 已完全失真。,为避免这种现象发生，可以对原方程组作等价变换，再利用顺序Gauss消元法求解。,写出原方程组的增广矩阵：,针对第一列找出绝对值最大的元素，进行等价变换：,求得方程的解为：x35.546，x245.76，x117.46,精确解为：x35.546，x245.76，x117.46,列主元Gauss消元法与顺序Gauss消元法的不同之处在于：,后者是按自然顺序取主元素进行消元,前者在每步消元之前先选取主元素然后再进行消元,下面将列主元Gauss消元法的计算步骤叙述如下：,给定线性方程组 Axb,记A（1）,b（1）A,b，列主元Gauss消去法的具体过程如下：,1.首先在增广矩阵A（1）,b（1）第一列的n个元素中选取绝对值最大的一个作为主元素，并把此主元素所在的行与第一行交换，即,2.其次进行第一步消元得到增广矩阵A（2）,b（2），在矩阵A（2）,b（2）第二列的后 n1个元素中选取绝对值最大的一个作为主元素，并把此主元素所在的行与第二行交换，即,3.再进行第二步消元得到增广矩阵A（3），b（3）。按此方法继续进行下去，经过n1步选主元和消元运算，得到增广矩阵A（n），b（n），它对应的方程组,A（n）xb（n）,是一个与原方程组等价的上三角形方程组，可进行回代求解。,容易证明，只要det（A）0，列主元Gauss消去法就可以顺利完成，即不会出现主元素为零或者绝对值太小的情形出现。,选列主元过程：,一、求主元 alk 使得|alk|=max|akk|,|ak+1,k|,|ank|;,二、判断：若|alk|，则输出错误信息并停机，否则 转三；,三、判断：若 l k，则交换增广矩阵中第 l 行与 k 行 的元素，否则不交换。,例5 用列主元法解,第一列中绝对值最大是10，取10为主元；,n 阶方程组第 k 轮消元时，选第 k 列的后(n-k+1)个元素中绝对值最大者作主元。,x1=0 x2=-1x3=1,第二列的后两个数中选出主元 2.5；,高斯消元法是一种顺序消元法。消元过程按方程和未知数的顺序依次作消元计算，缺乏一定的灵活性。计算中有除数为零时，算法便无法实现；即使所有的除数都不为零，可以计算出方程组的解，也不能保证计算结果有很高的精度。,m21=2/0.00001=2 105,=0.000003 106 0.40000 106=-0.40000 106,=0.000002 106 0.20000 106=-0.20000 106,x2=0.50000,x1=0.00000.,解法一：高斯消元法,m21=0.00001/2=0.5 10-5,=2-(0.5 10-5)3,x2=0.50000,x1=0.25000.,解法二：,列主元消元法,=0.20000 10 0.0000015 10=0.20000 10,=1-(0.5 10-5)2,=0.10000 10 0.000001 10=0.10000 10,2.2 矩阵的直接分解法,高斯消元法的矩阵形式,每一步消去过程相当于左乘初等下三角矩阵Lk,A 的 LU 分解(LU factorization),第一个方程组的系数矩阵为下三角矩阵，第二个方程组的系数矩阵为上三角矩阵，两个方程组都非常容易求解。,将A=LU 称为矩阵A的三角分解,这时线性方程组为：,令,则有,定理2：（矩阵的三角分解）设A为n n实矩阵，如果解AX=b用高斯消去法能够完成（限制不进行行的交换，即），则矩阵A可分解为单位下三角矩阵L与上三角知阵U的乘积。A=LU且这种分解是唯一的。,注：（1）L 为单位下三角阵而 U 为一般上三角阵的分解称为多利特尔（Doolittle）分解；（2）L 为一般下三角阵而 U 为单位上三角阵的分解称为库朗（Courant）分解。,Doolittle分解法：,根据 A=LU 有等式成立：,比较等式两端对应元素，有,可以解得：,当 i=1 时,当 j=1 时,当 i 1 时,当 j 1 时,下面，我们对具体矩阵进行Doolittle 三角分解。,为了表示和存储方便，可以将分解后的两个矩阵用一个矩阵表示,例7 利用Doolittle三角分解法分解矩阵,解：分解时用到如下公式,1,2,3,4,1,1,1,2,6,12,3,7,6,24,6,24,1,2,3,4,1,1,1,2,6,12,3,7,6,24,6,24,可以写成：,这时，矩阵的三角分解,如果我们要求解方程组,则由,得到,例8：利用Doolittle三角分解方法解线性方程组,解:进行三角分解ALU，可以对增广矩阵A,b作 三角分解:,1 2 3-2,-3,2,2,-3,-1,3,3,17,得到,1 2 3-2,-3,2,2,-3,-1,3,3,17,这时，相应的方程组为：,x135,x28，,例9：用矩阵分解方法解 AX=b,LY=bUx=Y,y1=6y2=-4y3=-4,x1=-13x2=8x3=2,A=LU,不选主元的LU分解,A=2，3，4；3，5，2；4，3，30 n=max(size(A);for k=1:n-1 d=A(k,k);for i=k+1:n m=A(i,k)/d;for j=k+1:n A(i,j)=A(i,j)-m*A(k,j);end A(i,k)=m;endend,A=2.00 3.00 4.00 1.50 0.50-4.00 2.00-6.00-2.00,平方根法：,在实际应用中，常见一类非常重要的线性方程组Axb，其中A为对称正定矩阵，即A是对称的且对任何非零向量 x 都有xTAx0。本节将对这类方程组导出更有效地三角分解求解方法，称之为平方根法。,设A为对称正定矩阵，那么A的所有顺序主子式均大于零，根据定理2.2，存在惟一三角分解 ALU，即,记 Ak（1kn）为A的 k 阶顺序主子阵，则det（Ak）为A的k阶顺序主子式。由上式，利用矩阵分块运算规则，容易验证,det（Ak）u11 u22 ukk,那么由 det（Ak）0，可知 ukk0，k1，2，n,这时，将上面的矩阵表示为：,即：A=LDM，其中 DM=U，M=D-1U。,当AAT为对称矩阵时，根据 ALDM,得到,ATMT D LT,再根据矩阵三角分解的唯一性,可知 MLT。于是,ALDLT,则有,令,如果对称正定矩阵A具有如下分解 AGGT，其中G为下三角矩阵，则称其为对称正定矩阵的 Cholesky（乔列斯基）分解。,为表示方便，可以记,给定对称正定方程组 Axb，对 A 进行 Cholesky分解ALLT，则原方程组等价于 LLTxb,Lyb LTx y,解此方程组即可得到原方程组的解x，这就是求解方程组的平方根法。,下面，我们通过比较矩阵的对应元素给出对称正定矩阵的平方根分解法。已知,即：LLTxb，等价于,比较对应元素：,当 i=j 时,当 j i 时,解得,于是，根据计算公式,可以对对称正定矩阵进行平方根分解,l11,l21,l22,l31,l32,l33,ln1,ln2,ln3,lnn,关于方程组 Ax=b,如果对系数矩阵进行了平方根分解 ALLT，则将方程组化为：Lyb,LTxy,解得,于是，关于系数矩阵是对称正定矩阵的线性方程组Ax=b的求解，分两步进行：,第一步：系数矩阵的平方根分解,第二步：解等价方程组,例10 用平方根法求解对称正定方程组,解:首先进行A 的Cholesky 分解,ALLT,2,-0.5,0.5,2,1.5,1,得 y12，y23.5，y31,得 x11，x21，x31,求解Lyb:,再求解 LTxy:,解三对角方程组的追赶法,考虑如下形式的三对角方程组问题,其中，系数矩阵满足条件:,该三对角方程组的增广矩阵 为:,第一步：第一行元素除以b1，取 1=c1/b1，y1=f1/b1，可得 增广矩阵：,第二步：从出发，用初等变换作 n-1 轮消元。作第 k 轮消元时，将矩阵中第 k 行元素乘以-ak+1 加到第 k+1 行元素上，然后将第 k+1 行主元单位化(k=1,2,n-1)。最后的增广矩阵：,在消元过程中所用除法次数为 2n-1，所用乘法次数为 2(n-1)。,根据矩阵初等变换的性质知，原三对角方程组,等价于如下方程组,对这一特殊的上三角方程组，在回代过程中只需用到 n 1 次乘法就可以求出方程组的解。,追赶法算法框图,例11 用追赶法求解方程组,最后得增广矩阵,所以，方程组的解为,a=0 1 1 1 1;b=4 4 4 4 4;c=1 1 1 1 0;f=2 1 1 1 2;d=b(1);f(1)=f(1)/d;for k=2:5,c(k-1)=c(k-1)/d;d=b(k)-a(k)*c(k-1);f(k)=(f(k)-a(k)*f(k-1)/d;endclear d a bfor k=4:-1:1,f(k)=f(k)-c(k)*f(k+1);endf,f=0.4808 0.0769 0.2115 0.0769 0.4808,MATLAB 程序,