线性代数方程组.ppt

第二章 线性方程组的数值解法 确定小行星轨道 以太阳为原点在轨道平面 内建立直角坐标系,取天文测量单位,在五个不同时间 观察小行星,测得坐标数据： x 4.5596 5.0816 5.5546 5.9636 6.2756 y 0.8145 1.3685 1.9895 2.6925 3.5265 通过计算确定椭圆方程 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 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; end end b(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); end b=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消元法中的主元素绝对值非常 小，使消元乘数绝对值非常大，计算过程中出现大数吃掉小 数现象，产生较大的舍入误差，最终导致计算解 x1 104.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| 1 时 当 j 1 时 下面，我们对具体矩阵进行Doolittle 三角分解。 为了表示和存储方便，可以将分解后的两个矩阵用 一个矩阵表示 例7 利用Doolittle三角分解法分解矩阵 解：分解时用到如下公式 1234 1 1 1 2612 3 7 6 24 624 1234 1 1 1 2612 3 7 6 24 6 24 可以写成： 这时，矩阵的三角分解 如果我们要求解方程组 则由 得到 例8：利用Doolittle三角分解方法解线性 方程组 解: 进行三角分解ALU，可以对增广矩阵A,b作 三角分解: 1 2 3 -2 -3 2 2 -3-1 3317 得到 1 2 3 -2 -3 2 2 -3-1 3317 这时，相应的方程组为： x135x28， 例9： 用矩阵分解方法解 AX = b LY = b Ux = Y y1 = 6 y2 = -4 y3 = -4 x1 = -13 x2 = 8 x3 = 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; end end 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 , 如果对系数矩阵进行了平方根分解 A LLT，则将方程组化为： Lyb , LTxy 解得 于是，关于系数矩阵是对称正定矩阵的线性方程组Ax=b 的求解，分两步进行： 第一步：系数矩阵的平方根分解 第二步：解等价方程组 例10 用平方根法求解对称正定方程组 解: 首先进行A 的Cholesky 分解 ALLT 2 -0.5 0.5 2 1.5 1 2-0.50.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 )。最后的增广矩阵 ： 其中， i = ci / (bi - aii-1 ) ( i = 2,3, , n-1 ) yi = ( fi aiyi-1) / ( bi - aii-1 ) ( i = 2, 3, , n) 在消元过程中所用除法次数为 2n - 1，所用乘法次数为 2(n - 1)。 根据矩阵初等变换的性质知，原三对角方程组 等价于如下方程组 对这一特殊的上三角方程组，在回代过程中只需用到 n 1 次乘法就可以求出方程组的解。 xn= yn k 从 (n-1) 到 1 循环 输出 x1，x2，xn ；结束 输入已知数据 ai，bi，ci，fi i 从 2 到 n 循环 y1 = f1 /b1，d = b1 i-1 = ci-1 /d; d = bi - aii-1; yi = (fi - aiyi-1)/d. xi = yi - i xi+1 追赶法算法框图 例11 用追赶法求解方程组 最后得增广矩阵 所以，方程组的解为 k12345 k0.25000.26670.26790.2679 yk0.50000.13330.23210.20570.4808 xk0.48080.07690.