线性方程组求解的数值方法.ppt

1,第二章 线性方程组求解的数值方法,第二章 线性方程组求解的数值方法,1.高斯消元法2.矩阵分解法3.向量范数与矩阵范数4.迭代法求解5.方程组的病态问题与误差分析,主要内容：,第二章 线性方程组求解的数值方法,理解各种线性方程组数值求解；掌握求解方法和解的误差分析方法；能编程实现求解算法。特别强调：遇到问题养成用计算机编程求解的习惯，不要习惯性的用笔算，而这是国内外学生的一个主要差距。,教学要求：,在自然科学和工程技术中，有很多问题的解决都需要用到线性方程组的求解。因此，求解线性方程组的问题是一个在科学技术中常见的普遍问题。,解线性方程组的数值解法：有直接法和迭代法两类。,直接法：计算过程没有舍入误差，经过有限次四则运算可求得方程组得精确解。（实际计算有舍入误差）高斯消元法，矩阵分解法,迭代法：核心是迭代求解的收敛条件和收敛速度。雅可比（Jacobi）迭代，高斯-赛德尔（Gauss-Seidel）迭代,基本思想方法：由行初等变换将系数矩阵约化为三角 矩阵；用回代的方法求解方程组。,例1 用消去法（高斯消元法）解方程组,高斯消元法是求解方程组的古典方法。,(2.1),3.1 高斯消元法,结论：,整个计算过程可分为两部分：（1）消元：把原方程组转化为系数矩阵为上三角矩阵的方程组；（2）回代：由系数矩阵为上三角矩阵的方程组求解,（2）回代求解，得：,解：（1）消元：,对于一般情形：n阶线性方程组的高斯消元法,若记,增广矩阵,（2.2）,（1）第1步(k=1)，,一般形式的高斯消元法：,设，首先对行计算乘数,对增广矩阵 进行行初等变换：,(即用 乘以2.2式的第1个方程，加到第i个方程上，消去2.2式的第二个方程直到第n个方程中的未知数),(代表第i行),得到与原方程组 等价的方程组。,记为,（2）一般第k步消元(),设已完成上述消元过程第1步，第2步，第k-1步，且,其中：,设，计算乘数,(即用 乘以2.2式的第k个方程，加到第i个方程上，消去2.2式的第k+1个方程直到第n个方程中的未知数),那么第k步消元操作即：,(3)继续这一过程，直到完成第n-1次消元，最后我们得到与原方程组等价的三角形方程组,(2.3)消元过程结束,求解三角形方程组2.3，得到求解公式,这个过程称为回代过程。,例题：考虑方程组,Gauss消去法中每步用来消去其他元素的 称为该步的主元素。Gauss消去法作为数值方法，主元素的选择是否会影响计算的结果呢？,则该方程的精确解为,而采用(，1)作为主元素，利用高斯消去法得到的解为,显然结果是错误的。,错误在哪个地方呢？原因是我们在消元时，利用了小主元，使得约化后的方程组元素数量级大大增长，再经舍入，而计算机的有效位数有限，造成消元后的三角形方程组就不准确了。,结论：在消元过程中可能出现主元素为零的情况，这时消去法将无法进行；即使不为零，在主元素很小时，用其做除数，也会导致其他元素数量级的严重增长和舍入误差的扩散，最后也使得计算解不可靠。,解决方法：对一般的矩阵来说，最好每一步选取系数矩阵（或消元后的低阶矩阵）的该列中绝对值最大的元素作为主元素，以使高斯消去法具有较好的数字稳定性。(高斯列主元素消去法),1.列主元法,第一列中绝对值最大是10，取10为主元,n阶方程组第k轮消元时，选第k列的后(n-k+1)个元素中绝对值最大作主元。,x3=6.2/6.2=1,x2=(2.5-5x3)/2.5=-1,x1=(7+7x2-0 x3)/10=0,x1=0 x2=-1x3=1,第二列的后两个数中选出主元 2.5,2 完全主元素消去法,整个矩阵中绝对值最大是10，取10为主元,n阶方程组第k轮消元时，选消元后元素中绝对值最大作主元。,x1=0 x2=-1x3=1,方框中6最大，交换行列，交换列的时候要做记录（即x3和x2交换了位置）：,完全主元素消除法与列主元素消除法的优缺点比较：,优点：数值更加稳定；缺点：计算量大；,对矩阵A实行初等行变换相当于用初等矩阵左乘A，于是对(2.2)做第一次消元后，化为 化为，即，其中,3.1 矩阵的三角分解 LU分解,第k步的初等矩阵为,并且,重复这一过程，最后得到,将上三角矩阵 记为U，则,将上三角矩阵 记为U，则，其中,则，L为单位下三角矩阵。,高斯消去法实质上是产生了一个将A分解为两个三角形矩阵相乘的因式分解。如果A是非奇异阵，则LU分解是唯一的，否则分解不唯一。,消元法：,消元法与三角分解法间的关系：,三角分解法：,讨论,直接三角分解法解线性方程组（）的具体流程：,1.,2.计算U的第r行，L的第r列元素 r=2,3,n,3.,(一)LU分解,再求解Ly=b,Ux=y计算公式：,(二)x的计算,例 用直接三角分解法解方程组,解:(一)矩阵LU分解,(1),故：,(2),经计算：,(二)求解x：,从而,3.2 矩阵的Cholesky分解,对称正定矩阵A：AT=A,A的各阶顺序主子式大于零.,前面指出非奇异的矩阵可以有三角分解，当A是某些特殊矩阵时，它的LU分解会有更加简洁的形式。,A的LU分解,(2.4),uii 0(i=1,n),由于A是对称正定的，则有,将U进一步分解：,根据A的对称性：,故：,由分解的唯一性知：,故：,其中P为上三角矩阵，这种对称正定矩阵的分解称为Cholesky分解。,在Matlab中函数“chol”给出对称正定矩阵的Cholesky分解。,经过n步可直接求得,思路:(1)分解对称正定矩阵,设n阶对称正定矩阵A有分解,先用待定系数法求L的元素,Cholesky分解的具体步骤(平方根法)：,(2)分解,求解方程组,Cholesky方法具体计算公式：,分解计算：,求解：,3.3 向量范数与矩阵范数,向量、矩阵与线性方程组有着密切的关系，向量、矩阵范数是解方程组以及研究与探讨方程组本身性质的工具。,一、向量范数的定义,定义 范数为,其中的，最常用的范数是，以及,容易证明向量范数满足如下性质：(1)如果，则；而且，当且仅当(2)对任何的数c，都有(3),范数也称为2-范数或Euclidean函数；范数也称为-范数或一致范数。在Matlab中用函数 用来计算向量x的 范数。,由性质(3),还容易得到：,二 矩阵范数的定义,在向量范数的基础上可以进一步定义矩阵范数,如果上式中的向量范数取为 范数，则可以证明定义出的矩阵范数（称为矩阵 范数或者列范数）为,如果向量范数取为2-范数，则,其中（模），称为矩阵B的谱半径。(为B的任意特征值。),类似地，Matlab函数norm(A,p)给出矩阵p=1,2或 范数,如果向量范数取为，则可以证明定义出的矩阵范数（称为矩阵的 范数或者行范数）为,容易证明矩阵范数满足如下性质：(1)如果，则；而且，而且仅当(2)对任何的数，(3)(4),而且由矩阵范数的定义还有称之为矩阵范数与相应的向量范数是相容的。,39,3.4 古典迭代法的构造,求解线性代数方程组的方法,中小规模问题,直接法,迭代法,大规模，超大规模问题,古典方法,现代方法,40,线性方程组的一般形式为：,如果 是非奇异的，则上式有唯一解。,我们将通过一个具体线性方程组的例子来讲解迭代法。取：,即线性方程组为：,方程的精确解(直接法计算),41,如果对该方程组进行变形，将变量 分别从三个方程中分离出来：,(1),令初值,所以(1)式可以表示为迭代形式：,所以我们可以得到一个向量的序列，只要该序列当 时有极限，那么这个极限就是该线性方程组的解。,上面这种迭代求解线性方程组的方法称为Jacobi迭代法。,考虑线性方程组的一般形式：,其中矩阵A为nn阶方阵，则方程的分量形式为：,改写成：,从而得到迭代公式：,上面式子就是一般形式的Jacobi迭代法，也可也记做：,在Jacobi迭代法中充分利用新值，可得到下面迭代：,上面这种迭代方法也叫做Gauss-Seidel迭代，也可以记为：,方程组的精确解为x*=(1,1,1)T.,解 Jacobi迭代法计算公式为,取初始向量x(0)=(0,0,0)T,迭代可得,例1：利用Jacobi法和Gauss-Seidel法求解线性方程组,可见,迭代序列逐次收敛于方程组的解,而且迭代7次得到精确到小数点后两位的近似解.,计算结果列表如下:,Gauss-Seidel迭代法的计算公式为:,同样取初始向量x(0)=(0,0,0)T,计算结果为,由计算结果可见,Gauss-Seidel迭代法收敛较快.取精确到小数点后两位的近似解,Gauss-Seidel迭代法只需迭代3次,而Jacobi迭代法需要迭代7次.,为了进一步研究,从矩阵角度来讨论上述迭代法.,对线性方程组Ax=b,记,D=diag(a11,a22,ann),则有 A=D-L-U,于是线性方程组 Ax=b 可写成(D-L-U)x=b,等价于 Dx=(L+U)x+b 或 x=D-1(L+U)x+D-1b,由此建立Jacobi迭代法迭代公式 x(k+1)=D-1(L+U)x(k)+D-1b k=0,1,2,或写成 x(k+1)=Bx(k)+f k=0,1,2,其中,所以Gauss-Seidel迭代法可以写成,其中,Gauss-Seidel迭代法迭代公式为,写成矩阵形式为：,3.5 迭代法的分析,构造迭代法的途径可以有很多，那么是不是所有的方法都能收敛呢？收敛速度的快慢与什么有关系呢？,它的精确解为,例：,可见,并不是所有的方程组都收敛，即使收敛对于不同的方法(例如Jacobi与Gauss-Seidel)收敛速度也是不同的.,计算结果列表如下:,设 是矩阵B任意一个特征值，是相应的特征向量，即,再假设 为任意的向量范数及与之相融的矩阵范数，则有：,首先引入几个定义、引理：,所以对矩阵B任意一个特征值，都有,记，称之为矩阵B的谱半径，则,先引入两个引理（详细证明见附录）：,引理3.1 矩阵，则 的充分必要条件是,引理3.2 矩阵，若，则 是非奇异的。,定理3.1 对任意的 和任意的初始向量，迭代法收敛的充分必要条件是,证明：必要性 设迭代法产生的序列 收敛，记 是该序列的极限点，所以 满足,又由迭代关系，可以得到,由 的任意性，知：,由引理3.1知：,充分性 由 及引理2知（即）有唯一解，记为,类似于必要性的证明，得到,再由第一步知，故,定理3.1给出了迭代收敛的充要条件，但 不宜计算，所以在实际使用中通常并不好用。由 知，只要对某种相容的矩阵范数有,那么当然，所以这常常是实际中很有效的收敛判别准则。,严格对角占优矩阵：,考虑，设，当 时的矩阵。,定理3.2 如果A是严格对角占优的矩阵，则它一定非奇异。,由于A是对角占优矩阵，所以，故矩阵,是可逆的。考虑,矩阵的 范数为,由A是对角占优矩阵，得,由引理3.2知 是非奇异的，又由于 是非奇异的，所以A是非奇异的。,得证！,引理3.2 矩阵，若，则 是非奇异的。,定理3.3 系数矩阵为严格对角占优的线性代数方程组，它的Jacobi迭代法和Gauss-Seidel迭代都是收敛的。,证明：Jacobi方法的迭代矩阵为,由定理3.2中的证明知，严格对角占优矩阵满足：,再由定理3.1，得Jacobi迭代法收敛。,(a)Jocobi的证明,1,再考虑Gauss-Seidel方法，其迭代矩阵为,用反证法，假设Gauss-Seidel迭代不收敛，则由定理3.1知：,所以存在模大于1的特征值.设有，使得行列式：,(b)Gauss-Seidel的证明,故，只有才能使得 成立。,由于A是对角占优的，所以矩阵,也是对角占优的，那么该矩阵一定非奇异，这与上面该矩阵 的行列式等于0矛盾。,由于 可逆，所以,得证！,定理3.4 若迭代矩阵B满足，则迭代生成的序列 满足,收敛速度问题,其中 表示精确解。,证明：,先证明(b),然后证明(a),而,所以,故，本页第一个式子可写为,由于，故,(b)得证,显然对于任意的正整数p，都有,那么在本页第一个式子中反复利用上式就可以得到(a),(a)得证,给出的关系可以作为计算终止的判别准则。即只要 和 已足够接近，表明 与 便已足够靠近。,定理3.4(a),给出的是迭代收敛速度的一个估计。显然 越接近0，生成序列 收敛得越快。,定理3.14(b),定理3.4的讨论,3.6 超松弛迭代(SOR)及分块迭代方法,对于给定的迭代法，每步迭代所需的工作量是确定的。如果迭代法收敛速度缓慢，则需要比较多的迭代次数，由此导致算法工作量太大而失去使用价值，因此各种迭代法的加速技术具有重要意义。,假设 是已经得到的迭代值，以 为初值进行一步Gauss-Seidel迭代得：,将这一迭代值与 的值组合作为新一步的迭代值,其中 为待定的参数，写成矩阵形式为：,这时迭代矩阵：,新的迭代每步的计算代价与Gauss-Seidel方法相差无几，如果能取得恰当的 值，使新构造的方法比Gauss-Seidel方法收敛得更快，松弛就起到了加速的作用。在实际上真正使用的 值通常的范围为，被统称为超松弛方法，简称SOR(Successive Over-Relaxation)方法。,上面的这个方法就叫做松弛方法，它可以视为Gauss-Seidel方法的加速。显然 时，上面这个方法就是Gauss-Seidel方法。,定理3.5 SOR方法收敛的必要条件是：,证明：,假设SOR方法收敛，所以,设 为 的特征值，则,另一方面，经过计算不难得到,故，结论得证。,分块矩阵以及分块迭代方法,在许多实际应用中，矩阵可以写成如下分块的形式：,其中 是 阶的方阵，而且 据此可以构造分块形式的迭代法，考虑,其中,求解方程组，分块的Jacobi迭代法为：,分块的Gauss-Seidel迭代法为：,分块的超松弛迭代法(BSOR)为,3.7 线性方程组的条件数,对于线性方程组 如果解x关于问题（即矩阵A和向量b）的“微小”变化(来源可能是模型误差或数据上的误差，也可能是数值计算中的舍入误差)不敏感，那么这个问题即是一个“好”问题，反之就是“坏”问题或者病态的问题。,关于病态的例子，我们在第一章中已经包含了一个，这一章主要用范数作为工具来分析线性方程组的好坏。,例 设有方程组,已知系数,测量数据,待求数据,其中A为，若测量数据为准确数据，即，那么很容易计算,如果测量数据不够精确，即测量数据仅保留了2位有效数字，即，计算结果为,如果测量数据和系数矩阵的数据都仅保留了2位有效数字，即,，则计算结果为,以上问题称为病态的问题，病态是问题本身固有的。,(1)考虑由于右端的扰动所引起的解的变化，,已知,两式相减得,利用范数关系可以得到：,综合上述两式有：,这个式子给出了右端的扰动可能引起解扰动的上界，显然 越小，方程右端的变化引起解的变化就越小。,(2)考虑由于矩阵A有扰动时所引起的解的变化,可以得到：,其中利用了,经过整理可以得到,(3)考虑更一般的情形所引起的解的变化,假设满足如下的特殊形式：,又Tylor展开，可以写成,对 求导得,令，并考虑到，则,记，从而有,从上面三种情况可见，问题扰动所引起解的扰动是被同一个因子 控制的,定义3.1 设A为可逆矩阵，是与某种向量范数相容的矩阵范数，称,为线性方程组 相应于该矩阵范数的条件数。,由前面的分析可见，条件数大的就是病态的矩阵，反之是良态的矩阵。求解线性方程时了解它的条件数大小是必要的，它可以帮助你判断所得数值解的可信性及模型的合理性。,在Matlab中可以用cond（A）来计算条件数或者用rcond(A)来估计其量级的倒数。,设，则，即对于该函数，误差会被放大n倍。,二、病态问题与条件数（针对问题本身）,定义：输入数据的微小变动导致输出数据的较大误差，就被称为病态问题。,衡量是否病态的标准：条件数,对于函数值计算问题，条件数定义为：,不同的问题，条件数具体定义不同。,一般情况下，条件数大于10，就认为问题病态。,通过构造特殊算法来解决,条件数的重要性质：,定理3.6(a)对任意的非奇异矩阵A，cond(A)是由任意的矩阵 范数定义的条件数，则,(b)如果Q是正交矩阵，则对于任意的非奇异矩阵A，都有,其中,证明：,(a)中的各条性质是显然的，例如，对任意的非奇异矩阵A，有,(b)由于Q是正交矩阵，所以,另外，,