第三章 解线性方程组的直接法其中且解线性方程组3.1 高斯消元法 高斯消元法：思 路(1) 将A化为上三角阵， (2) 回代求解。=具体地,可先看书中P35,例3.1. 例3.1 用顺序Gauss消去法解线性方程组解 1、消元过程：2、回代过程消元记Step 1：设 ，计算因子 将增广矩阵 第 i 行 mi1 第1行，得到 其中Step k：设 ，计算因子且计算共进行 ？ 步n 1回代定理 若A的所有顺序主子式 均不为0，则高斯消元无需换行即可进行到底，得到唯一解。（例3.3，充要条件）注注：事实上，只要 A 非奇异，即 A1 存在，则可通过逐次 消元及行交换，将方程组化为三角形方程组，求出唯 一解。例3.3的证明：必要性 若顺序Gauss消去法可行，即 则可进行消去法的第k-1步 由于 是由A逐行实行初等(消法)变换得到的，这些运算不改变相应顺序主子式的值，故有充分性 用归纳法证明。当k=1时显然成立。设命题对k-1成立。现设 由归纳法假设有 因此，消去法可以进行第k-1步，A约化为其中 是对角元为 的上三角矩阵，而且A的k阶顺序主子式与 的相等.function x=magauss(A,b,flag) %用途：顺序Gauss消去法解线性方程组Ax=b %格式：x=magauss(A,b,flag), A为系数矩阵, b为右端项, 若flag=0, %则不显示中间过程，否则显示中间过程, 默认为0, x为解向量 if narginA=1 1 1 1;-1 2 -3 1;3 -3 6 -2;-4 5 2 -5;b=10 -2 7 0;x=magauss(A,b);x得到计算结果：x=1 2 3 4 运算量 由于计算机中乘除运算的时间远远超过加减运算的时间 ，故估计某种算法的运算量时，往往只估计乘除的次数，而 且通常以乘除次数的最高次幂为运算量的数量级。 Gaussian 消去法:Step k：设 ，计算因子且计算共进行n 1步(n k) 次(n k)2 次(n k) 次(n k) (n k + 2) 次消元乘除次数：1 次(n i +1) 次 回代乘除次数：Gaussian 消去法的总乘除次数为，运算量为 级。3.2 选主元消去法 例：单精度解方程组/* 精确解为 和 */8个8个用Gaussian 消去法计算：8个小主元可能导致 计算失败。 全主元Gauss消去法 每一步选绝对值最大的元素为主元素，保证 。Step k: 选取 If ik k then 交换第 k 行与第 ik 行;If jk k then 交换第 k 列与第 jk 列; 消元注注：列交换改变了 xi 的顺序，须记录交换次序，解完后再 换回来。 列主元Gauss消去法 省去换列的步骤，每次仅选一列中最大的元。例：注注：列主元法没有全主元法稳定。 完全选主元法比 Gaussian Elimination多出 ，保证稳定，但费时。 列主元法比 Gaussian Elimination只多出 ，略省时，但不保证稳 定。算法.(列主元高斯消去法)步1：输入系数矩阵，右端b，置k:=1;步2：对k=1,n-1进行如下操作：(1)选列主元，确定 ，使若 ，则停止计算，否则，进行下一步；(2)若 ，交换 的第 两行；(3)消元：对 计算步3：回代对 计算function x=magauss2(A,b,flag) %用途：列主元Gauss消去法解线性方程组Ax=b %格式：x=magauss(A,b,flag), A为系数矩阵, b为右端项, 若flag=0, % 则不显示中间过程，否则显示中间过程, 默认为0, x为解向量 if nargink t=A(k,:); A(k,:)=A(p,:); A(p,:)=t;t=b(k); b(k)=b(p); b(p)=t;end%消元m=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);b(k+1:n)=b(k+1:n)-m*b(k); A(k+1:n,k)=zeros(n-k,1);if flag=0, Ab=A,b, end End %回代 x=zeros(n,1); x(n)=b(n)/A(n,n); for k=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k); end例3.5 利用通用程序magauss2.m计算下列方程组解 在MATLAB命令窗口执行A=2 -1 -3 1;-1 1 2 1 3; 2 3 3 -1;-3 1 3 2 4;1 3 -1 4 4;b=11 14 4 16 18;x=magauss2(A,b);x得到计算结果：x=1.0000 2.0000 1.0000 -1.0000 4.00003.3 解三对角方程组的追赶法Step 1: 追：其中与Gauss消去法类似，一旦= 0 则算法中断，故并非任何 三对角阵都可以用此方法分解。Step 2: 赶即回代。（注：总共只有6n-5次乘除计算量 ）定理则追赶法是可行的。若方程组(3.7) 的系数矩阵元素满足条件证: 只需证 即可，显然 且 设 ，则对function x=machase(a,b,c,d) %用途：追赶法解三对角方程组Ax=d %格式：x= machase(a,b,c,d) a为次下对角线元素向量，b主对角元素 % 向量，c为次上对角线元素向量，d为右端向量，x返回解向量 n=length(a); for k=2:nb(k)=b(k)-a(k)/b(k-1)*c(k-1);d(k)=d(k)-a(k)/b(k-1)*d(k-1); end x(n)=d(n)/b(n); for k=n-1:-1:1x(k)=(d(k)-c(k)*x(k+1)/b(k); end例3.6 利用通用程序machase.m计算下列三对角方程组解 在MATLAB命令窗口执行a=ones(50,1); b=4*ones(50,1); c=ones(50,1);d=6*ones(50,1); d(1)=5; d(50)=5; x=machase(a,b,c,d)得到计算结果：x=(1,1,1,.,1)定理若 A 为对角占优 /* diagonally dominant */ 的三对角 阵，且满足 ，则追赶 法可解以 A 为系数矩阵的方程组。它意味着矩阵的对角元是非常大.多大算大呢?它们满足下面的不等式:注注： 如果 A 是严格对角占优阵，则不要求三对角线上的 所有元素非零。 根据不等式可知：分解过程中，矩阵元素不会过分增大，算法 保证稳定。 运算量为 O(6n)。3.4 矩阵的三角分解 高斯消元法的矩阵形式 (书P45例3.8,设系数矩阵顺序主子式不为零) Step 1:记 L1 =，则Step n 1:其中Lk =记为L单位下三角阵记 U =A 的 LU 分解定理若A的所有顺序主子式均不为0，则 A 的 LU 分解唯 一（其中 L 为单位下三角阵）。证明：由前面分析可知，LU 分解存在。下面证明唯一性。 若不唯一，则可设 A = L1U1 = L2U2 ，推出上三角阵对角元为1 的下三角阵注注： L 为一般下三角阵而 U 为单位上三角阵的分解称为 Crout 分解。实际上只要考虑 A* 的 LU 分解，即 ，则即是 A 的 Crout 分解。为什么要讨论三角分解？若在消元法进行前能实现三角分解 ，则 容易回代求解 ,共有O(n2) 的计算量要确定矩阵 中的元素,则可利用公式 来列式以确定各个元素,具体可参看书中 P46,(3.12)-(3.17).或下页由A=LU,可得即算法3.3步1 输入系数矩阵A，右端项b；步2 LU分解：对 计算步3 解下三角方程组Ly=b对 计算步4 用回代法解上三角方程组Ux=y: 对 计算function x,l,u=malu(A,b) %用途：用LU分解法解方程组Ax=b %格式：x,l,u=malu(A,b) A为系数矩阵，b为右端向量，x返回 % 解向量，l返回下三角矩阵，u返回上三角矩阵format short %LU分解 n=length(b);u=zeros(n,n); l=eye(n,n);u(1,:)=A(1,:); l(2:n,1)=A(2:n,1)/u(1,1); for k=2:nu(k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k)/u(k,k); end %解下三角方程组Ly=b y=zeros(n,1); y(1)=b(1); for k=2:ny(k)=b(k)-l(k,1:k-1)*y(1:k-1); end x=zeros(n,1); x(n)=y(n)/u(n,n); for k=n-1:-1:1x(k)=(y(k)-u(k,k+1:n)*x(k+1:n)/u(k,k); end例3.9 利用通用程序malu.m计算下列方程组解 在MATLAB命令窗口执行A=2 -1 -3 1;-1 1 2 1 3; 2 3 3 -1;-3 1 3 2 4;1 3 -1 4 4;b=11 14 4 16 18;x,l,u=malu(A,b);x得到计算结果：x=1.0000 2.0000 1.0000 -1.0000 4.0000LU 分解与 Gauss 消去法的关系则取记 Lk =由此也可看出高斯消元法的消元过程相当于 的 求解,而回代过程相当于 的求解. 3.5 平方根法 /* Choleskiy分解 */：对称正定矩阵的分解法定义一个矩阵 A = ( aij )nn 称为对称阵，如果 aij = aji 。定义一个矩阵 A 称为正定阵，如果 对任意非零 向量 都成立。 回顾：对称正定阵的几个重要性质 A1 亦对称正定，且 aii 0若不然，则存在非零解，即 存在非零解 。 对任意 , 存在 , 使得 , 即 。其中第 i 位 A 的顺序主子阵 Ak 亦对称正定对称性显然。对任意 有, 其中 。 A 的特征值 i 0 设对应特征值 的非零特征向量 为 ，则 。 A 的全部顺序主子式 det ( Ak ) 0因为将对称 正定阵 A 做 LU 分解U =uij=u11uij / uii111u22unn记为A 对称即记 D1/2 =为什么 uii 0?因为det(Ak) 0则 仍是下三角阵定理设矩阵A对称正定，则存在非奇异下三角阵 使得 。若限定 L 对角元为正，则分解唯一。注注： 对于对称正定阵 A ，从 可知对任意k i 有 。即 L 的元素不会增大，误差可控，不需选主元。现在利用 来求出 的元素 如下：假设L的第k-1列元素已求得，下面求L的第k列元素由于上