1、第2章 线性方程组的直接解法,n阶线性方程组:,记:,则:,若,由Cramer法则,需计算n+1个行列式, 若采用逐行、逐列展开法, 需执行n! 次乘法. 若n较大, Cramer法则在实际计算时是不可行的.,解线性方程组的实用方法有两类：直接法, 迭代法.,本章介绍直接法.,一. 高斯(Gauss) 消去法,普通的Gauss消去法,先将方程组化为同解的上三角方程组 (称为消元过程);,第2章 线性方程组的直接解法,基本思想,再按 xn,xn-1,x2,x1的顺序求出各个解 (称为回代过程).,例, , ,消元过程:,例（续）,第2章 线性方程组的直接解法,回代过程:,n阶线性方程组的Gaus
2、s消去法求解步骤,消元过程: (采用增广矩阵的方法描述),初始元素用上标(1)记.,第1步: 消去该列元素.,设,消 时, 记,用 , 得,其中,消元过程（续-1）,第2章 线性方程组的直接解法,第2步: 消去该列元素.,设,消 时, 记,用 , 得,第k步: 消去该列元素.,消 时, 记,消元过程（续-2）,第2章 线性方程组的直接解法,第n-1步: 消元后,系数矩阵为上三角阵.,回代过程:,按变量的反序计算.,实际计算时, 元素 可以放在A矩阵的 中, 可放在b的 中, 即每次消元后, 新的矩阵元素可以替代原来的值.,输入n, A=(aij), b=(b1,bn)T. for k=1 to
3、 n step 1 以第k行对角元为主元进行消元 2.1 if akk=0, 输出计算失败信息, 停机. 2.2 for i=k+1 to n step 1 对k+1n行循环(消元aik)计算 aik/akkaik 存放公式中的 likfor j=k+1 to n step 1 对第i行的k+1n列循环(消元)计算 aij - aik*akjaijend for ( j ).计算 bi - aik*bkbiend for ( i ).end for ( k ). for k=n to 1 step -1 回代 3.1 计算 bk/akkbk 3.2 for i=1 to k-1 step 1
4、: bi-aik*bkbi end for (i).(移第k列)end for ( k ). 输出解x (=b),Gauss消去法的算法描述,第2章 线性方程组的直接解法,算法2.1 伪代码,这种做法能适用于分块求解.,Gauss消去法的计算量,第2章 线性方程组的直接解法,2. for k=1 to n step 1 n次循环 2.2 for i=k+1 to n step 1 n-k次循环aik/akkaik 每次k循环: n-k次除法for j=k+1 to n step 1 n-k次循环aij-aik*akjaij 每次k循环:(n-k)2次乘end for ( j ).bi-aik*
5、bkbi 每次k循环: n-k次乘法end for ( i ).end for ( k ).,只统计乘除法的次数.,消元过程:,除法次数:,乘法次数:,回代过程:,乘除法次数:,总的乘除法次数:,Gauss消去法时间复杂度: O(n3),for k=n to 1 step -1 n次循环 3.1 计算 bk/akkbk 每次k循环: 1次除法 3.2 for i=1 to k-1 step 1 k-1次循环bi-aik*bkbi 每次k循环: k-1次乘法end for (i).(移第k列)end for ( k ).,二. 主元素法,列主元素法,第2章 线性方程组的直接解法,消元过程:,Ga
6、uss消去法要求对角元(称为主元素) , 故只适用于1n-1阶顺序主子式均非0的矩阵. 若主元素非常小,则误差较大,因此Gauss消去法的数值稳定性较差 (参见书上p15的例子).,第1步,找出绝对值最大者,两行 交换,第1次 消元,第2步,找出绝对值最大者,第2次消元,列主元素法（续）,第2章 线性方程组的直接解法,第k步,在矩阵 的第k列中找,交换第k行与ik行, 进行第k次消元, 得,经过上述n-1步消元后, 矩阵A化为上三角阵.,回代过程:,与普通Gauss法相同, 回代得方程的解.,上述过程中, 主元是在列元素中选取的, 列主元素法由此得名. 另外, 选择列主元过程中需交换两行元素,
7、 这仅仅相当于交换了两个方程的次序, 解的次序并未改变.,全主元素法,第2章 线性方程组的直接解法,消元过程:,第1步,在所有元素中, 找出绝对值最大者,两行 交换,两列交换,第1次 消元,第2步,如此经过n-1步消元后, 将A化为上三角阵.,找出绝对值最大者,第2次消元,交换第2行与i2行; 交换第2列与j2列.,注意: 矩阵列元素的交换, 会改变解的排列次序!,回代过程:,同普通的Gauss消去法.,全主元素法例子,第2章 线性方程组的直接解法,交换 行列,x1,x2交换,x2 x1 x3,回代求解:,第1步,第2步,第3步,消元:,全主元素法程序实施的说明,第2章 线性方程组的直接解法,
8、全主元法计算精度高, 数值稳定性好. 但计算中要交换行列元素, 程序稍为复杂, 计算时间较长.,为避免交换元素, 将消元和回代时所用的行列号看成逻辑编号, 分别用i, j记, 而实际行列号分别用数组Lin(i)和Col(j)保存.,利用临时数组很快能得到按原来逻辑顺序排列的解. 首先将解b(i) 存放到临时数组c(i)中(i=1,2,n); 对于第i 个逻辑行, 实际行号ir=Lin(i), 其解存放于b(ir)中, 它所对应的解的实际序号为jr=Col(i). 令b(jr)=c(ir) 即可完成解的整理工作.,当需要交换逻辑行i与j时, 只要交换Lin(i) 与Lin(j)的值即可; 交换逻
9、辑列i与j时, 只要交换Col(i) 与Col(j) 的值即可. 不必交换系数矩阵A的元素以及右端项b的元素, 可节约交换元素的时间.,逻辑行号i由ir=Lin(i)给出, 逻辑列号j由jr=Col(j) 给出, 求解过程中用到aij元素时, 只要用A(ir, jr)即可, bi 则用b(ir) 代替.,按原来的逻辑顺序排列结果,全主元素法程序实施的例子,第2章 线性方程组的直接解法,全主元素法程序实施的例子(续),第2章 线性方程组的直接解法,结果次序整理, i=1:,b(3)=0, i=2:,b(2)=1, i=3:,b(1)=2,全主元法的Fortran程序,三. 直接三角分解法,Gau
10、ss消去法的矩阵形式,第2章 线性方程组的直接解法,消元第1步,消元过程相当于作一系列的初等行变换. 如取,其中,Gauss 消去法,只改变了第2行元素,第1列元素被消去,取,Gauss 消去法的矩阵形式 (续),第2章 线性方程组的直接解法,第k步,其中,第k列元素被消去,经过n-1步消元:,即,取,Doolittle分解式的导出,第2章 线性方程组的直接解法,记为,由,容易验证,记,则,结果,注意:,Doolittle分解时, 右端项并不参与, 因此与右端项无关; 有了Doolittle 分解式, 解线性方程组变得非常容易.,矩阵的三角分解,第2章 线性方程组的直接解法,定理2.1 设A为
11、n阶方阵, 若A的顺序主子式det(Aj) (j=1,2,n-1)均不为零, 则A存在唯一的Doolittle分解.,证明: 采用数学归纳法证明存在性.,使得,当j=1时, 若det(A1)=a110, 则存在,若结论对j=k-1成立, 即,其中,各阶主子式:,定理2.1 证明 (续-1) 存在性,第2章 线性方程组的直接解法,存在,使得,并且,于是结论对j=n-1成立, 即,存在 , 使得,即,存在性得证.,定理2.1 证明 (续-2) 唯一性,第2章 线性方程组的直接解法,用反证法: 若有两种Doolittle分解,唯一性得证. #,由于A非奇异, 因此 非奇异,如何进行三角分解,其中,
12、记,三角分解的求解步骤,第2章 线性方程组的直接解法,编程时可合起来写.,上面的计算必须按照一定的次序才能完成.,仍用A矩阵存放L,U矩阵的元素. 例如, 计算u3n时,a3n在计算完u3n后就不再使用, 因此可用来存放u3n.,同理,an3在计算完 ln3后就不再使用, 因此可用来存放 ln3.,三角分解的求解步骤（续）,第2章 线性方程组的直接解法,次序1: 逐行从左到右计算;,次序2: 逐列从上到下计算;,次序3: 实际编程时, 也可逐垮先计算u元素, 再计算L元素;,例. 求下面矩阵A的Doolittle分解,按逐垮计算方式的次序,直接三角分解法解线性方程组,第2章 线性方程组的直接解
13、法,记 则,该方法的时间复杂度与Gauss消去法相同, 都是O(n3/3).,前消,回代,解三对角方程的追赶法,第2章 线性方程组的直接解法,称A为三对角矩阵, 方程组为三对角方程组.,这样, 三角分解过程比满阵时大为简化.,于是, 三角分解时, L,U的形式为:,由Gauss消去法可知, 每一步消元过程只需消去1个元素. 如第1步只需消去a2, 第2步只需消去a3,三对角矩阵三角分解的存在性和唯一性由下面的定理给出.,三对角矩阵的三角分解定理,第2章 线性方程组的直接解法,定理2.2 设三对角矩阵A满足下列条件:,若,则称为严格对角占优.,则它可分解为:,元素ci都是矩阵A中的原始元素.,且
14、分解是唯一的.,证明: 首先, 由A=LU两边对应的下三角元素相等, 可得,下面采用数学归纳法证明存在性.,若ui0, 则分解式存在(下面证明).,故, 分解式唯一.,A的(i,i+1)元素, 即ci自然满足.,定理2.2 证明 (续-1),第2章 线性方程组的直接解法,由于,及,第1次消元前, u1=b10, 故存在初等矩阵,使得,其中,假设第k次消元前,故存在初等矩阵,使得,其中,定理2.2 证明 (续-2),第2章 线性方程组的直接解法,由,及,A可以消元成为上三角阵, 即 , 定理得证.,由数学归纳:,即 ,可继续消元.,从证明中还可看出, 最后一个条件可放松为:,因为, 若 , 取k
15、=n-1,追赶法解线性方程组,第2章 线性方程组的直接解法,第1步: 三角分解,第2步: 求解 Ly=d,消元过程 “追“的过程.,第3步: 回代过程 “赶“的过程,追赶法的具体实施,第2章 线性方程组的直接解法,算法2.2 追赶法计算步骤,输入n, A的系数(a2,an),(b1,bn),(c1,cn-1), 右端项(d1,dn). for i=2 to n step 1 消元过程ai/bi-1ai 用ai保存下三角元素libi-ci-1aibi 用bi保存上三角元素ui, 因此u1=b1不变di-aidi-1di 用di保存消元结果yiend for ( i ). 3. dn/bndn 得
16、到第n个解xn 4. for i=n-1 to 1 step -1 回代过程(di-cidi+1)/bidi 得到第i个解xiend for ( i ) 5. 输出方程的解d1,d2,dn,注:,计算过程中用向量a,b,c来保存系数矩阵A的元素, 并保存下三角和上三角元素, 这可以节约存储单元; 可证明, 当系数矩阵严格对角占优时, 此方法具有较好的数值稳定性.,追赶法的Fortran程序,四. 平方根法与改进的平方根法,平方根法 (Cholesky分解法),第2章 线性方程组的直接解法,若系数矩阵A对称正定, 则A的各阶顺序主子式全大于0, 特征值全是实数, 且大于0. 从而可导出下面的方法
17、.,定理2.3 对称正定矩阵A, 则有唯一的分解式:,证明: A对称正定 A的各阶顺序主子式全大于0. 由定理2.1 Doolitle分解:, P是单位上三角阵, 且有:,令,由A的对称性,由Doolitle分解的唯一性,或,对,由A的正定性,即D正定.,D的对角元大于0.,令,则,定理2.3 证明 (续),第2章 线性方程组的直接解法,记 为下三角阵, 对角元为 的元素, 因此都大于0.,于是,唯一性: 假定存在下三角阵GL, 对角元都大于0. 且有,于是,两端左乘L-1, 右乘G-T,的第i个对角元,的第i个对角元,唯一性得证.,Cholesky分解式的具体计算,第2章 线性方程组的直接解
18、法,对A的下三角元素(ij),这是按行进行计算的公式. 书上p30 (2-30)式是按列计算的公式.,关于求解L中各元素的顺序,求,求,用Cholesky分解法求线性方程组,第2章 线性方程组的直接解法,记,则,前消过程:,第i个方程:,回代过程:,第i个方程:,乘除法运算量n3/6, 约为Gauss消去法的一半, 但需作n次开方运算.,改进的平方根法 (LDLT法),第2章 线性方程组的直接解法,在上面定理2.3的证明中, 对称正定矩阵A, 有:,且对角元素大于0,计算A的下三角元素:,这是按列顺序的计算公式, 也可按行顺序计算.,如果仍用A矩阵存放L和D, 则计算 所涉及的元素如右图.,虽
19、然避免了开方运算, 但乘法次数比Cholesky分解增加约一倍, 总运算量为n3/3量级.,实际上, 有许多乘法是重复的, 可以用数组保存起来, 这样可减少计算工作量.,for i=2 to n step 1 按行顺序循环. i=1时,d1=a11for j=1 to i step 1 计算(i,j)位置的元素值.,改进的平方根法的具体计算,第2章 线性方程组的直接解法,定义一维数组uk (k=1,n). 仍用A来保存元素 .,end for ( j ). end for ( i ).,当j1时, 计算,if ji, 计算,采用上面的方法, 乘除法的运算次数与Cholesky分解相当, 同时避
20、免了开方运算.,按行的顺序计算L的第i 行元素时:,for k=1 to j-1 step 1 计算end for (k).,for k=1 to j-1 step 1 计算 End for (k).,按列的顺序计算L的第j列元素时:,第2章 线性方程组的直接解法,for j=2 to n step 1 按列的顺序循环, 从2n列.,定义一维数组uk (k=1,n). 仍用A来保存元素 .,for k=1 to j-1 step 1 计算end for (k).,end for ( i ). end for ( j ).,for i=2 to n step 1 计算第1列元素. i=1时,d1
21、=a11计算 End for ( i ).,if ij, 计算,for i=j to n step 1 计算(i,j)位置的元素值.,改进的平方根法解线性方程组,第2章 线性方程组的直接解法,由于利用了系数矩阵的对称性, 平方根法和改进的平方根法计算量只有Gauss消去法的一半, 而且其数值稳定性良好, 是求解中小型对称正定线性方程组的好方法.,则,前消过程:,第i个方程:,回代过程:,第i个方程:,规定:,该方法不仅可用于稠密系数矩阵, 也可用于稀疏的系数矩阵. 如, 在有限元方程的求解中, 经常用LDLT法来求解. 在求解过程中可保持L矩阵的稀疏性. 另外, 在以前计算机内存较小的情况下,
22、 还发展了对系数矩阵进行分块分解的求解方法.,改进平方根法的Fortran程序,五. 误差分析,向量和矩阵的范数,第2章 线性方程组的直接解法,向量的范数向量大小的一种度量,向量和矩阵的范数在线性方程的误差分析中起着十分重要的作用.,定义2.1 对任意向量 ,按某种规则对应于一实数 . 若 满足:,1非负性: 且 当且仅当,则称 为向量 的范数.,2齐次性: 对任意实数, 都有,3三角不等式:,范数,在Rn中引进函数 其中p1为正整数,容易验证, fp(x)满足条件 再根据Minkowski不等式:,立即推得满足条件 因此,是Rn中的范数, 称为 范数, 记为,常用的3种向量范数,第2章 线性
23、方程组的直接解法,(1) 2-范数,(2) 1-范数,(3) -范数, 向量的模长, 又称Euclid范数.,在数值分析中, 当p=1,2,三种情形时特别有用.,当n=3时: 2-范数的单位球为Euclid单位球; 1-范数单位球是以(1,0,0),(0,1,0),(0,0,1)为顶点的4-双棱锥; -范数的单位球是以原点为中心, 边长为2的立方体.,当 时, 称 为单位向量.,定义集合 为单位球.,向量范数的等价性,第2章 线性方程组的直接解法,Rn中的两个范数 和 , 若存在实数m,M0, 对 都有, 则称两种范数是等价的.,不难证明, 1-范数,2-范数和-范数都是等价的.,例如, 由C
24、auchy-Schwarz不等式,对任意xi,yi成立. 若取xi|xi|, yi=1, 则,另外,故1-范数和2-范数是等价的.,矩阵范数,第2章 线性方程组的直接解法,1 且 当且仅当,则称 为矩阵A的一种范数.,2,3,定义2.1 对任意n阶矩阵A, 在 中定义一个实值函数, 记作 , 若 满足:,4,矩阵范数与向量范数的相容性,定义 设在 中定义了一种矩阵范数 , 在 中定义了一种向量范数 . 若对任意n阶矩阵A及 中的向量x, 恒成立不等式:,则称矩阵范数 与向量范数 是相容的.,可以证明: 对任一种矩阵范数, 至少存在一种向量范数, 使得它们是 相容的.,矩阵范数 (续),第2章
25、线性方程组的直接解法,诱导范数,定理2.4 设A是n阶方阵, 是Rn中的向量范数, 则 是一种矩阵范数, 称为向量范数 诱导出的矩阵范数.,反之, 若 , 则 对任意 成立,. 对任意实数, 有,. 对任意方阵A和B,. 对任意n维非零向量y,即,于是,满足矩阵范数与向量范数的相容性.,常用的矩阵范数,第2章 线性方程组的直接解法,(3) 2-范数,(1) 1-范数,(2) -范数,上面3种都是由对应的向量范数诱导的矩阵范数., 又称矩阵的列范数., 又称矩阵的行范数., 又称谱范数. 为矩阵ATA的最大特征值.,(4) F-范数, 相当于向量2-范数的直接推广, 但它不是向量诱导的矩阵范数.
26、,因为, 如果F-范数是诱导范数, 则对于单位阵I, 有,而,(Frobenius),F-范数与向量的2-范数相容(利用Cauchy-Schwarz不等式可证). 即,若A是对称矩阵, 则 , 为A的按模最大特征值.,矩阵的误差表示,第2章 线性方程组的直接解法,方程组的状态与条件数,第2章 线性方程组的直接解法,方程组病态的概念,实际问题中数据是有误差的, AA+A, bb+b.,如果A,b 很小, 而x 很大, 则称方程组是“病态”的.,病态方程的例子,准确解,近似解,误差非常大, 方程组是“病态”的.,右端项扰动的影响,第2章 线性方程组的直接解法,取范数,另外,系数矩阵A扰动的影响,移
27、项,若 足够小, 使得,解的相对误差与 有关. 一般地, 越大, 解的扰动也越大.,因此, 控制了解的扰动程度.,条件数的概念,第2章 线性方程组的直接解法,系数矩阵A的条件数越大, 方程病态就越严重.,特别地, 若A是实对称矩阵,解的相对误差可改写成:,定义2.3 对非奇异矩阵A, 称 为矩阵A的条件数, 记为,系数矩阵A扰动时,常用的条件数,(3),(2),(1),其中, 是A的特征值, 且,条件数的性质,第2章 线性方程组的直接解法,证. 对正交阵A, ATA=1, 其所有特征值都是1.,证.,(1) 对任意n阶方阵A, 都有:,对任意由向量诱导的矩阵范数 (如1-,2-,-范数),(2
28、) 对任意实数0:,(3),(5) 对正交矩阵A:,(4),对任意矩阵范数 : 若A是任意n阶非零方阵,证.,病态矩阵的例子,第2章 线性方程组的直接解法,Hilbert矩阵是有名的病态矩阵,n=3时,当n=6时,求条件数需要求矩阵的逆, 其工作量较大. 根据数值经验, 下列情况下方程组通常是病态的.,病态方程的数值经验,(1) 在用主元法(列主元或全主元)消元时出现小主元;,(2) 系数矩阵有行(列)近似线性相关, 或行列式近似为零; 但有例外, 如A=I, 是小量, det(A)=n0, 但方程状态良好.,误差分析,第2章 线性方程组的直接解法,近似解的误差估计方法,例,取近似解,通常,
29、如果 很小, 就认为x*较精确.,但对于病态方程, 即使 很小, x*的误差仍然很大.,残量,但, 解的误差: 却非常大.,近似解精确程度的判断定理,第2章 线性方程组的直接解法,证明.,由此看出, 当方程病态时, cond(A)很大, 即使残量很小, 但解的误差仍会很大.,定理2.5 设x是方程Ax=b的精确解, x*是它的一个近似解, r=b-Ax*是近似解的残量, 则,关于算法数值稳定性的讨论,第2章 线性方程组的直接解法,不算大.,采用Gauss消去法求解, 消去A的(2,1)元素, 得,回代得,误差却非常大!,Gauss消去法数值稳定性较差, 列主元和全主元Gauss消去法的数值稳定
30、性较好.,定义: 对于条件数不太大的系数矩阵A, 采用某一给定的算法求解方程, 若初始数据的误差及计算过程中的舍入误差的积累、传播,对计算结果的影响较小, 则称该算法是数值稳定的. 反之, 称为数值不稳定的.,六. 超定线性方程组的最小二乘解,问题的提出,第2章 线性方程组的直接解法,目标: 寻找一组近似解, 使得每个方程的残量的平方和极小, 即:,定义 如果方程的个数多于变量个数, 则称为超定线性方程组.,通常, 不存在一组解: (x1,x2,xn), 使得它精确满足所有m个方程.,要使得F(x)取到极小. 得到的解称为最小二乘解.,这类问题称为最小二乘(LS)问题.,最小二乘问题的解,第2
31、章 线性方程组的直接解法,若记C=ATA, 则,即,如果A的秩为n (列满秩),ATA对称正定,如, 正交化方法、广义逆方法, 不作详细介绍. 下面仅给出其定义.,有唯一解.,其它方法,是ATb的第k行元素.,其中 是ATA矩阵的第(k,j)元素,即Cx 的第k行的值.,因此, 方程化为,广义逆的定义,第2章 线性方程组的直接解法,则称B是A的Moore-Penros广义逆, 记作A+, 即A+=B.,ABA=ABAB=B(AB)T=AB(BA)T=BA,定义2.4 设A为mn阶矩阵, B为nm阶矩阵. 如果满足:,可以证明: 当A为列满秩时, A+是由A唯一确定的, 且:,于是, 最小二乘问题 的解可表示为:,
