数值分析第5章 解线性方程组的直接法.doc_第1页
数值分析第5章 解线性方程组的直接法.doc_第2页
数值分析第5章 解线性方程组的直接法.doc_第3页
数值分析第5章 解线性方程组的直接法.doc_第4页
数值分析第5章 解线性方程组的直接法.doc_第5页
已阅读5页,还剩45页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

第五章 线性代数方程组的数值解法线性方程求解问题是科学研究和工程计算中最常见的问题。如电学中的网络问题、工程力学中求解连续力学体(微分方程)问题的差分方法、有限元法、边界元法及函数的样条插值、最小二乘拟合等,都包含了解线性方程组问题。因此,线性方程组的解法在数值计算中占有极其重要的地位。对于阶线性方程组,若,则方程组有惟一解。由克莱姆(Cramer)法则,其解为其中为用向量代替中第列向量所得矩阵。每个阶行列式共有项,每项都有个因子,所以计算一个阶行列式需做次乘法,我们共需要计算个行列式,要计算出,还要做次除法,因此用Cramer法则求解要做次乘除法(不计加减法),计算量十分惊人。如时,就需作约次乘法。可见Cramer法则在理论上是绝对正确的,但当较大时,在实际计算中却是不可行的。因此寻求有效的数值计算方法就成为非常必要的课题。线性方程组的类型很多,若按其系数矩阵阶数的高低和含零元素多少,大致可分为两类:一类是低阶稠密线性方程组,即系数矩阵阶数不高(例如,),含零元素很少。另一类是高阶稀疏线性方程组,即系数矩阵阶数高,零元素占绝对优势(比如占以上)。线性方程组的数值解法也可分为两大类:直接法和迭代法。直接法是在没有舍入误差的情况下,通过有限步运算可以得到方程组精确解的方法。但是,在实际计算时,由于初始数据变为机器数而产生的误差以及计算过程中所产生的舍入误差等都要对解的精确度产生影响,因此直接法实际上也只能算出方程真解的近似值。目前常用的有效算法是Gauss消去法和矩阵的三角分解法。迭代法是用某种极限过程去逼近准确解的方法。如对任意给定的初始近似解向量,按照某种方法逐步生成近似解序列使极限为方程组的解。因为在实际计算时,只能做到有限步,所以得到的也是近似解。常用的迭代法主要有Jacobi迭代法、Gauss-Seidel迭代法、逐次超松弛法及共轭斜量法等。直接法的优点是运算次数固定,并且可以事先估计。它的缺点是通常需要存储系数矩阵和常数项的所有元素,因而所需存贮单元较多,编写程序较复杂,一般适用于求解低阶线性方程组;迭代法的优点是原始系数矩阵始终不变,且只需存储原方程组系数矩阵中的非零元素,因而算法简单,编写程序较方便,且所需存贮单元也较少,所以被广泛用于求解高阶稀疏线性方程组。缺点是存在收敛性和收敛速度问题,因而只能对具有某些性质的方程组才适用。1 消元法1。1 三角方程组的解法形如 (5-1)的方程组称为上三角形方程组。写成矩阵形式为。若,则(5-1)有惟一解我们称求解上三角方程组(5-1)的过程为回代过程。同时也看到求解这类方程组是件容易的事,所以对一般形式的方程组,应设法将它化为(5-1)式的形式,然后再求解。1。2 Gauss消去法考虑方程组 (5-2)其矩阵形式为其中化线性方程组(5-2)为等价的三角形方程组的方法有多种,由此导出不同的直接方法,其中Gauss消去法是最基本的一种方法。Gauss消去法分消元计算和回代求解两个过程。为后面的符号统一起见,记方程组(5-2)为其中。消元计算第一步:就是要将的第一列主对角元以下的元素全约化为0。设,计算用乘(5-2)的第1个方程,加到第个方程上,完成第一步消元,得(5-2)的同解方程组 (5-3)简记为,其中的元素的计算公式为假设前步消元完成后,得(5-2)的同解方程组为 (5-4)简记为。第步:就是要将的第列主对角元以下的元素全约化为0。设,计算用乘(5-4)的第个方程加到第个方程上,完成第步消元。得同解方程组其中元素的计算公式为完成步消元后,(5-2)化成同解的上三角方程组 (5-5)简记为回代求解因,故,于是 (5-6)Gauss消去步骤能顺利进行的条件是,现在的问题是矩阵应具有什么性质,才能保证此条件成立。若用表示的顺序主子式,即有下面定理定理1 约化的主元素的充要条件是矩阵的顺序主子式。证明 必要性。因主元素,可进行步消元,每步消元过程不改变顺序主子式的值,于是必要性得证。用归纳法证明充分性。时命题显然成立。假设命题对成立。设。由归纳法假设有,Gauss消去法可以进行步,约化为其中是对角元为的上三角阵。因是通过步消去法得到的,每步消元过程不改变顺序主子式的值,所以的阶顺序主子式等于的,即由知,充分性得证。1。3 Gauss消去法的计算量1)消元过程的工作量第步消元过程:计算乘子需作次除法运算;第行乘加到第行,需要乘法和减法各次(第列元素无需计算),共行(),所以需要乘法次,故消元过程中的乘、除和减法运算量为乘除法次数:加减法次数:2)回代过程的工作量计算需要次乘除法和次加减法运算,整个回代过程的运算量为:乘除法:减法:所以Gauss 消去法的总运算量为乘除法:(当较大时)加减法:(当较大时)当时,用Gauss消去法约需430次乘除法运算,而用Gramer法则约需次乘除法。1。4 Gauss消去法的矩阵形式如果用矩阵形式表示,Gauss消去法的消元过程是对方程组(5-2)的增广矩阵进行一系列初等变换,将系数矩阵化成上三角形矩阵的过程,也等价于用一串初等矩阵去左乘增广矩阵,因此消元过程可以通过矩阵运算来表示。第一次消元等价于用初等矩阵左乘矩阵,其中,即一般地,第次消元等价于用初等矩阵左乘矩阵,其中,经过次消元后得到即将上三角矩阵记为,得到其中为单位下三角阵。这说明,消元过程实际上是把系数矩阵分解为单位下三角阵与上三角矩阵的乘积的过程。这种分解称为杜利特尔(Doolittle)分解,也称为分解。定理2 设为阶方阵,如果的顺序主子式,则存在惟一的单位下三角阵和上三角阵,使。证明 以上的分析已证明了可作分解。下面证明分解的惟一性。如果非奇异,设有两个分解式其中都是单位下三角阵,都是上三角阵。因非奇异,所以都可逆。于是上式左边为单位下三角阵,右边为上三角阵,所以即有,惟一性得证。若为奇异阵,因其阶顺序主子式不等于零,故可记其中均为阶方阵,且非奇异。由,得由已证明的结果知的分解是惟一。所以亦是惟一确定的。进而,也是惟一确定的。定毕。定理2的逆命题是定理3 设阶方阵非奇异,若有惟一的分解(其中为单位下三角阵,为上三角阵),试证的顺序主子式。(证明留给读者)2 主元素法前面已指出Gauss消去法必须在各个约化主元素下才能进行。现在还需指出的是:即使,但当很小时,也不能做主元素的,因为作第次消元时,需将第个方程乘以,因此当很小时,乘数很大,用去乘第行的元素,将导致的数量级迅速增长,这样在消元计算 时,会出现大数吃掉小数的现象,因而导致最后计算结果的精度很差甚至失真。例1 解方程组(1) 求精确解;(2) 在的浮点机上用Gauss消去法求解。解 (1) 容易求得方程组的精确解为。(2) 在所给浮点机上原方程组为由于消去第二式的,得对阶,出现大数“吃掉”小数,结果有解得,回代第一式得。与精确解相比较,已无精确度可言。产生不准确的原因是主元素太小,致使很大。改变上述状况的办法是将方程组的一、二两式对换,得然后再使用Gauss消去法,此时于是得近似解。此例表明,高斯消去法解方程组时,小主元可能带来严重的后果,因此应尽量避免小主元的出现;另一方面,通过对换方程组中方程的次序或改动变量次序,选择绝对值大的元素做主元,可减少舍入误差,提高计算精度。1 列主元与全主元消去法1列主元素消去法。在第步消元时,在的第列元素中选取绝对值最大者作为主元,并将其对换到位置上,然后再进行消元计算。用方程组(5-2)的增广矩阵 (5-7)表示它,并直接在增广矩阵上进行计算。具体步骤如下:第一步:首先在上述矩阵的第一列中选取绝对值最大的,比如,则。将(5-7)中第一行与第行互换。为方便起见,记行互换后的增广矩阵为,然后进行第一次消元,得矩阵假设已完成步的主元素消去法,约化为第步:在矩阵的第列中选主元,如,使。将的第行与第行互换,进行第次消元。经过步,增广矩阵(5-7)被化成上三角形算法1(Gauss列主元素消去法)说明:消元结果冲掉,行乘数冲掉,det存放行列式值。输入,置,1。 消元过程对(1) 选主元:(a) 按列选主元,即确定,使(b) 若,停机;(c) 若(进行行交换)(2) 对()对(3) 回代过程(a) 若,输出失败信息,停机(b) (c) 对注:计算程序中对的判断可以用或(是预选的很小正数)。2完全主元素法。在第步消元时,在的右下方阶矩阵的所有元素中,选取绝对值最大者作为主元,并将其对换到位置上,再做消元计算。全主元消去法和列主元消去法相比,每步消元过程所选主元的范围更大,故它对控制舍入误差更有效,求解结果更加可靠。但全主元法在计算过程中,需同时作行与列的互换,因而程序比较复杂,计算时间较长。列主元法的精度虽稍低于全主元法,但其计算简单,工作量大为减少,且计算经验与理论分析均表明,它与全主元法同样具有良好的数值稳定性,故列主元法是求解中小型稠密线性方程组的最好方法之一。2 列主元消去法的矩阵形式第步消元时,需先交换的第行与第行,然后再做消元计算,相当于对进行如下的矩阵计算:(1) 用初等排列矩阵左乘,即其中由交换单位矩阵的第与两行所得。显然。(2) 用初等排列矩阵左乘,即于是对施行带行交换的步消元过程,可用矩阵表示为从而有因为,所以上式可改写为可以证明,若是指标为的单位下三角阵,则仍是一个指标为的单位下三角阵。若令则是排列阵。又设于是即注意到仍为单位下三角阵,若令则有上式表明,带行交换的Gauss消去过程所产生的矩阵分解,相当于对系数矩阵先施行每步消元时所做行交换后,再将所得矩阵进行分解。以上讨论可叙述为定理4(列主元三角分解定理)如果为非奇异矩阵,则存在排列矩阵,使其中为单位下三角矩阵,为上三角矩阵。3 Gauss-Jordan消去法高斯消去法始终是消去对角线下方的元素,现考虑高斯消去法的一种修正,即消去对角线下方和上方的元素,这种方法称为Gauss-Jordan消去法。设用高斯-若当消去法已完成步,于是化为等价方程组,其中在第步计算时,考虑对上述矩阵的第行上、下都进行消元计算。1。 按列选主元元素,即确定使换行(时)交换第行与第行元素3。 计算乘数(,且),(可保存在存放的单元中)。4。 消元计算(,且,)(,且)计算主元()经过步消元后,可得说明用高斯-若当方法将约化为单位矩阵,常数项位置就是计算解,无需回代过程。用高斯-若当方法解方程组计算量大约需要次乘除法,要比高斯消去法大,但用高斯-若当方法求一个矩阵的逆矩阵还是比较合适的。定理2(高斯-若当法求逆矩阵)设为非奇异矩阵,方程组的增广矩阵为。如果对应用高斯-若当方法化为,则。事实上,为非奇异矩阵,若我们已通过左乘一系列初等矩阵把它化成了一个单位矩阵,即则说明用同样的初等阵左乘单位矩阵,即得到的逆矩阵。3 直接三角分解法3。1 不选主元的直接三角分解法设为如下形式由矩阵的乘法规则,得由此可得计算和的公式 (5-8)具体步骤如下:1 计算的第1行,的第1列2 计算的第行,的第列例2 求矩阵的三角分解。解 按式(5-8)所以紧凑格式:例2中矩阵的三角分解按紧凑格式计算,结果见下表如果线性方程组的系数矩阵已进行三角分解,则解方程组等价于求解两个三角形方程组。第一步:先求解下三角方程组得解第二步:再求解上三角方程组得解3。2选主元的直角三角分解法不选主元的直接三角分解过程能进行到底的条件是,实际上,即使非奇异,也可能出现某个的情况,这时分解过程将无法进行下去。另外,如果但很小,会使计算过程中的舍入误差急剧增大,导致解的精度很差。但如果非奇异,我们可通过交换的行实现矩阵的分解,实际上是采用与列主元消去法等价的选主元的三角分解法,即只要在直接三角分解法的每一步引进选主元的技术即可。设步分解已经完成,这时有考虑步分解,对于做到(4)(1) 计算辅助量,冲掉(2) 选主元。即确定行号,使,并记录主行号:(为一维数组)(3) 交换的行与行元素(4) 计算的行与的第列元素(这里)上述计算过程完成后,就实现了的分解,保存在的上三角部分,保存在的严格下三角部分,排列阵由的最后记录可得。例如,设,则于是方程组的求解等价于解,即,可通过及完成。例3 用选主元的直接三角分解法解方程组,其中,解(1) 对作选主元的分解,中间结果冲掉相应位置元素,数组记录主行号。第一步分解:,因为,于是有第二步分解:,因为,故有于是由选主元过程,。,(2) 解,得;解,得3。3 解三对角方程组的追赶法在数值计算中,如三次样条插值或用差分方法解常微分方程边值问题,常常会遇到求解以下形式的方程组其中,系数矩阵 (5-9)是一种特殊的稀疏矩阵,它的非零元素集中分布在主对角线及其相邻两条次对角线上,称为三对角矩阵。方程组称为三对角方程组。Gauss消去法用于三对角方程组时过程可以大大简化。由于第次消元时,只需消去,所以消元过程就是用单位下二对角矩阵左乘矩阵。由于每个以上的列元素为零,不难看出具有形式定理5 设矩阵(5-9)满足下列条件 (5-10)则它可分解为 (5-11)其中为矩阵(5-9)中所给出,且分解是惟一的。证明 将式(5-11)右端按乘法法则展开,并与的元素进行比较,得如果,那么我们可计算,即 (5-12)从以上的公式和消元过程我们可看出,要使得分解得以实施,必须满足现在我们用归纳法证明:。当时,显然成立。假定成立,即,我们将证明时也成立。有当矩阵(5-9)按式(5-12)计算进行分解后,求解方程组可化为求解方程组及。解得 (5-13)再解,得 (5-14)按上述过程求解三对角方程组成为追赶法。式(5-12),(5-13)结合称为“追”的过程,相当于Gauss消去法中的消元过程。式(5-14)称为“赶”的过程,相当于回代过程。算法21)输入2)对 3)4)对 5)输出,停机。追赶法的基本思想与Gauss消去法及三角分解法相同,只是由于系数中出现了大量的零,计算中可将它们撇开,从而使得计算公式简化,也大大地减少了计算量。为节省计算机存贮单元,计算得到的分别存放在的存贮单元内,而存放在的存贮单元内。可以证明,当系数矩阵为严格对角占优时,此方法具有良好的数值稳定性。4 平方根法与改进的平方根法工程技术中的许多实际问题所归结的线性方程组,其系数矩阵常具有对称正定性,对于具有此类特性的系数矩阵,利用矩阵的三角分解法求解是一种较好的有效方法,这就是对称正定矩阵方程组的平方根法及改进的平方根法,这种方法目前在计算机上已被广泛应用。4。1平方根法(Cholesky分解法)定理6(对称阵的三角分解定理)设是对称阵,且的所有顺序主子式均不为零,则存在惟一的单位下三角阵和对角阵,使 (5-15)证明 因为的各阶顺序主子不为零,由定理2,存在惟一的Doolittle分解其中为单位下三角矩阵,为上三角矩阵。令,将再分解为其中为单位上三角矩阵。于是又由分解的惟一性即得,式(5-15)得证。定理7(对称正定阵的Cholesky分解)设是对称正定矩阵,则存在惟一的对角元素为正的下三角阵,使 (5-16)证明 因为对称性,由定理6知,其中为单位下三角阵,。若令,则为的Doolittle分解,的对角元即的对角元。不难验证,的 ()阶顺序主子式为对应的与的阶顺序主子阵的乘积,因此的顺序主子式。因为正定,有,由此可推出。记则有其中,它为对角元为正的下三角阵,所以(5-16)成立。由分解的惟一性,可得分解(5-16)的惟一性。证毕。分解式称为正定矩阵的乔列斯基(Cholesky)分解。利用Cholesky分解来求系数矩阵为对称正定矩阵的方程组的方法称为平方根法。用比较法可以导出的计算公式。设比较与的对应元素,可得 (5-17)这里规定。计算顺序是按列进行的,即。当矩阵完成Cholesky分解后,求解方程组就转化为依次求解方程组它们的解分别为 (5-18)当的元素求出后,的元素也就求出,所以平方根法约需次乘除法,大约为一般分解法计算量的一半。同时还节省存贮,只需存系数矩阵的下三角部分和右端项,中间结果冲掉,计算解冲掉。此外,由式(5-17)的第一式得所以上式表明,在矩阵的乔列斯基分解过程中的平方不会超过的最大对角元,舍入误差的放大受到了控制,且对角元素恒为正数,于是不选主元素的平方根法是数值稳定的,计算实践也表明了不选主元已有足够的精度,所以平方根法是目前计算机上解决这类问题的最有效的方法之一。4。2改进的平方根法(法)利用平方根法解对称正定线性方程组时,计算矩阵的元素时需要用到开方运算。另外,当我们解决工程问题时,有时得到的是一个系数矩阵为对称但不一定是正定的线性方程组,为了避免开方运算和求解对称(未必正定)方程组,我们可以引入下面改进平方根法。由定理6我们知道,为对称矩阵时,它可分解成其中为单位下三角矩阵,为对角矩阵。记,由矩阵乘法运算,并注意到,得于是导出分解的计算公式:对 (5-19)计算顺序为:。按式(5-19)进行分解,虽然避免了开方运算,但在计算每个元时多了相乘的因子,故乘法运算次数比Cholesky分解约增多一倍,乘法总运算量又变成数量级。仔细分析式(5-19)可以看出,式中有许多计算是重复的,如为此我们将引进辅助量,于是式(5-19)可改写成 (5-20)具体计算过程是:对矩阵作分解后,解方程组可分两步进行:首先解方程组,再由求出。具体公式为 (5-21)例4 用改进平方根法求解方程组解 容易验证,系数矩阵为对称正定阵。按式(5-20)计算分解式,得按式(5-21)计算方程组的解,得所以方程组的解为。5 向量和矩阵的范数为了研究线性代数方程组近似解的误差估计和迭代法的收敛性,我们需要引入衡量向量和矩阵“大小”的度量概念向量和矩阵的范数概念。1 向量范数定义1 设对任意向量,按一定的规则有一实数与之对应,记为,若满足正定性:,而且当且仅当;齐次性:对任意实数,都有;三角不等式:对任意,都有则称为向量的范数。以上三个条件刻划了“长度”、“大小”及“距离”的本质,因此称为范数公理。对上的任一种范数,显然有。向量空间上可以定义多种范数,常用的几种范数:1)向量的1-范数:;2)向量的5-范数:;3)向量的-范数:;4)更一般的-范数:。容易证明,及确实满足向量范数的三个条件,因此它们都是上的向量范数。此外,前三种范数是-范数的特殊情况()。事实上,我们只需表明3)。事实上,由于及,故由数学分析的夹逼定理有。定理8 设给定,则对上每一种向量范数都是的元连续函数。证明 对于给定的,设为的列向量,将写成分块形式对,由三角不等式有其中,所以对任意的,当时,有这就证明了的连续性。推论 是的元连续函数。下面讨论范数的等价性问题定义2线性空间上定义了两种范数和,如果存在常数,使 (5-22)则称和是上等价的范数。显然,范数的等价性具有传递性,即若与等价,与等价,则有与。定理9 上所有范数是彼此等价的。证明 只要证明上任一种范数都与等价。设是上定义的任一范数,记它是上的有界闭集。由定理8的推论,是上的元连续函数,故在上有最大值和最小值,且。现在考虑,且,则有,所以有从而而当时上式自然成立,这就证明了与的等价性。2 矩阵的范数这里主要讨论中的范数及其性质(的可类似讨论),其范数要符合一般线性空间范数的定义1,为了考虑矩阵乘法运算的性质,我们在矩阵范数的条件中多加一个条件。定义3 如果对上任一矩阵,按一定的规则有一实数与之对应,记为。若满足 ,且当且仅当; 对任意实数,都有; 对任意的两个阶方阵,都有; (相容性条件)。则称为矩阵的范数。这里与向量范数是一致的,条件则使矩阵范数在数值计算中使用更为方便。例5(矩阵的Frobenius范数)证明满足矩阵范数定义3。它可以看成维向量的5-范数,所以满足定义3的前三个条件,可验证条件也成立。记,利用矩阵的乘积及Cauchy不等式有即得证。由于在实际中,经常用到矩阵与向量的乘积运算,为了估计矩阵与向量相乘积的范数,那么在矩阵范数与向量范数之间具有某种协调关系是有好处的。为此,我们先引入相容性的概念,再由一种向量范数导出对应的矩阵范数。定义4 对于给定的上一种向量范数和上一种矩阵范数,若有 (5-23)则称上述矩阵范数与向量范数相容。不难验证如下不等式:对于上的一种向量范数,对任一,对应一个实数,下面我们定理表明它定义了上的一种矩阵范数。不难验证它有等价的形式 (5-24)定理10 设是上任一种向量范数,则对一切,由式(5-23)确定的实数定义了上的一种范数,把它记为,且有 (5-25)证明 首先证明(5-23)中的可以写成。因为由定理8,是中有界闭集上的连续函数,故在上有最大值,即存在,使,所以式(5-24)中的上确界就是最大值。记为,有式(5-25),并且由此可以得到 (5-26)余下将验证以上确定的满足范数定义3。条件,是明显的。对任意,有即得证,同理可证。定义5 对于上任意一种向量范数,由(5-25)所确定的矩阵范数,称为从属于给定向量范数的矩阵范数,简称从属范数(也称由向量范数诱导出的矩阵范数,自然范数或算子范数)。显然,单位矩阵的任一种从属范数。所以当时,不是从属范数。显然从属范数与所给定的向量范数相容,但是矩阵范数与向量范数相容,却未必有从属关系。例如可以证明与相容,即但不从属于。我们把从属于向量-范数、1-范数及5-范数诱导的矩阵范数,并分别称为矩阵的-范数、1-范数及5-范数。定理11 设,则1(称为的行范数)2(称为的列范数)3(称为的5-范数)其中是矩阵的最大特征值。证明 只就1,3给出证明,2同理。设,则所以有另一方面,设整数满足令,其中显然有,而且3对任意,从而为非负定的对称阵,其特征值为非负实数,依次排列为,对应一组规范正交的特征向量,对任意的,可表示为如果,则有,特别地,取,上式等号成立,故。定义6 设的特征值为,称为的谱半径。引理1 设为上任一种向量范数,从属于它的矩阵范数记为,为非奇异阶方阵,证明(1) 定义了上一种范数。(2) 从属于的矩阵范数,记为,有(证明留给读者)定理12 (1)设为上任一种(从属或非从属)矩阵范数,则对任意的,有 (5-27)(2)对任意的及实数,至少存在一种从属范数,使 (5-28)证明 (1)设满足,且。必存在向量,使不是零矩阵。对于任意一种矩阵范数,由矩阵范数定义可得即可推出(5-27)。(2) 对任意的,总存在非奇异的,使为Jordan标准形,即是对角块形式的矩阵,其中的Jordan块为对于实数,定义对角阵容易验证仍为块对角形式,其分块与J相同,即其中它的阶数与相同。取的范数,可得而为非奇异阵,由引理1知,定义了上一种向量范数,从属于此向量范数的矩阵范数为定理13 设是上一种从属范数,矩阵,满足,则为非奇异矩阵,且 (5-29)证明 若奇异,存在非零向量,使,所以有一个特征值1,故有。根据定理12,有,这和定理假设条件矛盾,所以非奇异。记,则于是。同理可证的情形。6 误差分析6。1 方程组的性态与条件数在用数值计算方法解线性方程组时,计算结果有时不准确,这可能有两种原因:第一是计算方法不合理;另外一种情况可能是线性方程组本身的问题。判断一个计算方法的好坏,可用方法是否稳定、解的精确度高低以及计算量、存贮量大小等来衡量。然而,对于不同的问题,同一方法却可以产生完全不同的效果,这就涉及所提问题的性态,即“好、坏”。例6 容易看出,方程组的解为。而方程组的解为。比较这两个方程组可以看出,只是右端项有微小的差别,最大相对误差为,竟使解的结果面目全非。如果我们把第二个方程组看作是第一个方程组的常数项经微小扰动得到的,所以也可称第二个是第一个的扰动方程组。可见方程组的解对方程组的初始数据的扰动十分敏感,这种性质与求解方法无关,是方程组本身固有的性质,即方程组的性态决定的。定义7 如果矩阵或常数项的微小变化,引起方程组解的巨大变化,则称此方程组为“病态”方程组,矩阵称为“病态”矩阵(相对于方程组而言),否则称方程组为“良态”方程组,称为“良态”矩阵。一般地说,在计算机上解方程组,实际上都是解所给方程的扰动方程,因为至少在对输入的数据作十进制与二进制转换时就要产生舍入误差。对于良态问题,由于它的解与其扰动方程的解差别不大,所以只要算法是数值稳定的,就可以得到较好的结果,而对于病态问题,即使算法是稳定的,其计算结果依然会很坏。以下我们研究方程组的系数矩阵和向量的微小误差(扰动)对解的影响。设为任何一种向量范数,矩阵范数是从属范数。设为非奇异且有微小扰动,有微小扰动,则的扰动方程为 (5-30)由,得由于非奇异,所以矩阵存在,因此为了估计,对上式两端取范数,则有移项得假定足够下,使得,则相对误差为利用,得到 (5-31)式(5-31)表示了解的相对误差与系数矩阵的相对误差与右端向量的相对误差之间的关系。由此得到:定理14 设是非奇异阵,。若系数矩阵及右端项分别有微小误差及,引起解的误差,当时,它们的相对误差间有如下关系推论1 若,则 (5-32)推论2 若,且充分小,则 (5-33)事实上,记,对于给定的,为确定的数,扰动充分小,可使,于是,所以(5-31)可以写成不等式(5-31)、(5-32)及(5-33)式给出的都是解的相对误差的上界。推论1、推论2分别指出了当只有或的误差时,解的相对误差都不超过它们的相对误差的倍数。刻划了线性方程组的解对初始数据误差的敏感度,此数越大则在很小的或下可能使解的相对误差很大,从而大大破坏了解的准确性,而且是方程组本身一个固有的属性,它与如何求解方程组的方法无关。因此可以用它来表示方程组的性态。定义8 设是非奇异阵,称数为矩阵的条件数,用表示,即矩阵条件数与范数有关,通常使用的条件数有:(1) ;(2) 的谱条件数当是对称矩阵时其中与为的绝对值最大和最小的特征值下面定理给出条件数的一些性质,证明留给读者。定理15 设,非奇异,则有(1),。(2)为正交阵,则;若为正交阵,则(3)设与为按绝对值最大和最小的特征值,则若对称,则。定理15表明,矩阵的条件数不小于1,当条件数接近1时,矩阵称为“良态”的;当条件数比1大得多,矩阵是“病态”的。正交矩阵是最稳定的一类矩阵。例7 计算例6方程组系数矩阵的条件数。解 系数矩阵为其逆矩阵为于是有条件数很大,所以方程组是“病态”的。例8 已知希尔伯特(Hilbert)矩阵计算的条件数。解 的逆矩阵的元素是所以(1)计算条件数,所以。同样,可计算,。当愈大时,病态愈严重。(2)考虑方程组的扰动方程(及的元素取3位有效数字)有简记为。方程与它的扰动方程的解分别为:,。于是,而,有,所以,这表明与的相对误差不超过,而引起解的相对误差超过50%。Turing在1948年第一个使用了“条件数”这一术语,接着不少数学工作者相继提出用“矩阵条件数”来度量方程组求解等问题的病态程度,但是矩阵条件数大到怎样才算大,并没有一个客观的衡量标准。另外,当较大时, 的计算也并非易事。在实际计算中,如果遇到下列几种情况,就应当考虑方程组可能是病态的。(1) 用主元素消去法解方程组时,出现小主元;(2) 系数矩阵某两行(列)几乎线性相关;(3) 系数矩阵的元素间数量级相差很大,且无规律;(4) 近似解向量已使剩余向量的范数很小,但解仍不符合客观规律。在求得方程组的一个近似解后,检验精度的一个简单方法是将代入方程组求得残量(余量)。如果很小,就认为解比较准确。但对“病态”严重的方程组,可能即使残量很小,近似解与准确解的差仍很大。如例6中,方程组的准确解为。若以作为它的近似解,其残量很小,但解的误差却不小。事实上,由和得于是可推出上式说明,尽管已很小,如果很大,近似解得相对误差也会很大,所以方程组还是病态的。6.2 病态方程组的解法病态方程组的求解问题是计算方法的一个重要的课题。在这里我们仅就不是过分病态的方程组,讨论改善解的精度,有如下几种方法:(1)采用高精度的算术运算(如采用双倍字长进行运算),以便改善和减轻病态矩阵的影响;(2)当方程组的系数矩阵元素的数量级相差悬殊时,计算经验表明,即使使用主元素消去法,求解过程中有效数字也会严重损失,影响解的精度。为避免这一现象发生,可将系数矩阵元素的数量级大体均衡一下,比如先用系数矩阵每行元素的最大模除遍行或列各元素,称为行或列标度化过程,然后施行主元素消去法。适当选择非奇异对角阵,使求解的问题转化为求解等价方程组其中的条件数得到改善。记。多数场合,选择的原则是使的元素满足和常常,为了简单,只对方程组进行行标度化,即只选择。例9 设简记为,计算。解 ,有当用列主元消去法求解时(计算到三位有效数字)于是得到很坏的结果:。取,考察等价方程组的系数矩阵的条件数。则,大大改善了系数矩阵的条件数。再用列主元消去法求解,得到从而得到较好的计算解:。以上对方程组的系数矩阵进行预处理的过程,称为显式标度化方法。它在计算时,会引进新的误差,实践中使用如下的隐式标度化法比较普遍。隐式法与显式法的不同之处在于它不算出,列主元素消去法在本质上依旧对,而不是,只不过在第高斯消去法时,用选择列主元。例10 解方程组用舍入的三位浮点数计算。解 方程组的系数矩阵元素之间的数量级相差很大,我们将用隐式标度化法求解。方程组的增广矩阵为1消元计算(1) 计算(2) 按列选主元,即确定,使于是,即应为主元素,所以对换的两行:(3) 消元计算:2回代求解通过此例看到,标度化的目的是找到真正的主元,它只影响主元素的选取。本例的准确解是。(3) 迭代改善法。设有方程组,其中为非奇异矩阵,且若用Gauss消去法(或直接三角分解法或)求得计算解(精度不高),为获得方程组高精度的解,一般可采用下述方法改善解的精度。(1) 计算剩余(用双倍字长计算,再舍入成单字长的数)(2) 求解(用单字长运算)。求出修正量。(3) 计算改善解。(4) 如果(指定的精度),则就是满足精度要求的近似解;否则对重复上述过程(1)(4),就得及,如此做下去,即可求得方程组的一个近似解序列。当方程组不是过分病态时,通常会很快收敛到准确解。称上述方法为迭代改善法。例1 用直接三角分解法解方程组解 因方程组的系数矩阵的两行(列)几乎线性相关,所以方程组可能是病态的。事实上方程组的确是病态的。利用迭代改善法:(1)实现分解得(2)解及得(3)计算(用双倍字长计算)(4)求解得(5)计算改善解(6)计算。重复(3)(6)步做法,得上述迭代改善过程还可以继续下去,如果给定,则就是满足精度要求的近似解。方程组的准确解。7 矩阵的正交三角化及应用7。1 初等反射阵定义9 设满足,则称为Housholder矩阵(变换),或初等反射矩阵。定理16 Housholder矩阵具有以下性质:(1)是实对称的正交矩阵,即。(2)。(3)仅有两个不等的特征值,其中是重特征值,是单重特征值,为其相应的特征向量。图5.1(4)对任意的,有从图5。1可以看出,与关于超平面对称,反射变换的名称由此而来。证明 (1)显然是实对称矩阵,又所以是实对称的正交矩阵,且。(2)(3) (4)定理17 设为两个不相等的维向量,且,则存在一个初等反射矩阵,使。证明 令,则得到一个初等反射阵而且因所以定理17表明,一个向量经某一个反射变换后可变成任何方向的等长向量。定理18(约化定理)设,则存在初等反射阵使,其中证明 记,设,取,则有,于是由定理17知,存在变换其中,使。记,于是其中,。显然如果和异号,哪么计算时有效数字可能损失,我们取和有相同的符号,即取在计算时,可能上溢或下溢,为了避免溢出,将规范化(设)则有使,其中算法6 设,本算法计算及使,分量冲掉的分量。计算 关于的计算,设其中为的第列向量,则,因此计算就是要计算于是计算只需计算两向量的数量积和两向量的加法。计算只需作次乘法运算。7.2 平面旋转矩阵设,则变换其中为正交矩阵。是平面上向量的一个旋转变换。中变换:,其中 ,及称为中平面的旋转变换(或Givens变换),称为平面旋转矩阵。显然,具有如下性质:(1) 与单位阵只是在位置元素不一样,其他相同。(2) 为正交矩阵()。(3) (左乘)只需计算第行与第行元素,即对有 其中。(4)(右乘)只需计算第列与第列元素 利用平面旋转变换,可使向量中的指定元素变为零。定理19(约化定理)设,其中不全为零,则可选择平面旋转阵,使其中, 。证明 由,知只要选择满足可使。评注:实际应用中,我们只关心之值,并不需要算得。保证是目的,因此可同时取以上值的相反数,相应得。7.3 矩阵的分解以下我们讨论如何用正交矩阵来约化矩阵。设为非零矩阵,将矩阵按列分块(1) 第1步约化:若,取,即这一步不需要约化。不妨设,于是可选取初等反射阵,使,于是其中是维零列向量,。(2) 第步约化:假设已对完成上述前步的约化,即存在初等反射阵使其中为阶上三角阵,是阶零矩阵,。将写成按列分块的形式,则有 不妨设,否则这一步不需要约化(如果列满秩,则)。于是,可选取阶初等反射阵,使根据构造阶矩阵就有而其中,是行的零列向量, 。因此其中为阶上三角阵。令,继续上述过程,最后有 (上梯形)第步约化大约需要次乘法运算。以上讨论可归结为如下结果:定理20(矩阵的正交约化)设且,则存在初等反射阵,使为上梯形阵。约化过程大约需要()次乘法运算。定理21(矩阵的分解)(1) 设且的秩为,则存在初等反射阵使其中为阶非奇异上三角阵。(2) 设为非奇异矩阵,则有分解其中为正交矩阵,为上三角矩阵。且当具有正对角元时,分解唯一。证明 (1) 由定理20可得。(2) 由假设及定理20存在初等反射阵使记,则上式为,即其中为正交矩阵,为上三角矩阵。唯一性。设有,其中为正交阵,为非奇异上三角阵,且具有正对角元素,则由假设及对称正定矩阵的Cholesky分解的唯一性,则,从而可得。下面考虑用平面旋转变换来约化矩阵。定理22 (用Givens变换计算矩阵的分解)设为非奇异矩阵,则(1) 存在正交矩阵,使为非奇异上三角矩阵。(

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

最新文档

评论

0/150

提交评论