绪 论
0.1数值计算方法与算法
1. 准确解和数值解
数学是在寻求准确解(解析解)中不断发展 数的定义
数的运算规则
2. 计算方法是用计算机求解数学问题数值解
理论上有解,无计算公式
例如,计算6阶矩阵的全部特征值。
有计算公式,无法承受计算量
解一个有100个未知量的线性方程组; 3.科学计算和计算方法
第三种科学方法:理论、实验、科学计算
计算方法是一门理论性和实践性都很强的学科,计算方法既有数学类课程中理论上的抽象性和严谨性,又有实用性和实验性的技术特征。计算方法的前提课程是微积分,线性代数,常微分方程和一门计算机语言。 0.2 误差与有效数字 绝对误差与绝对误差界
定义: 设*
x 为精确值(或准确值),x 是*
x 的一个近似值,
称 x x e -=*为近似值x 的绝对误差或误差。 绝对误差 = 精确值 - 近似值
定义:如果精确值*x 与近似值x 的误差的绝对值不超过某正数ε,即 ε≤-=||||*x x e
称ε 为绝对误差限或误差限。(例:迭代序列终止控制) 例:若经四舍五入得到 456.123=x ,
对于数 4564.123,4561.123,4555.123,4559.123的近似值都是456.123=x ,
它的误差限是:
34
*
102
1
510
||||--=?≤-=x x e
若0123456.0*
=x ,则它的误差限是: 78
*102
1
510||||--=?<-=x x e
相对误差与相对误差限
定义:设*
x 为精确值,x 是*
x 的一个近似值
称 *
**x x
x x e e r -== 为近似值x 的相对误差。
在实际计算中,有时得不到精确值*
x ,当r e 较小时*
x 可用近似值x 代替,即
x
x
x x e e r -==*
相对误差 =绝对误差精确值 或 相对误差 =近似值
绝对误差
相对误差r e 的值也可正可负,与绝对误差一样不易计算,常用相对误差限控制相对误差的范围。
定义:如果有正数r ε使得 r r x
e
e ε≤=||
* ,则称r ε为*x 的相对误差限。 产生误差的因素很多,产生误差的原因主要有 原始误差
由客观存在的模型抽象到物理模型产生的误差。包括模型误差和原始数据误差。 截断误差
用有限项近似无限项时,由截取函数的部分项产生的误差,称为截断误差。
例如:∑∞==+++++=02!!21n n k x
n x k x x x e ,在计算中用∑∑∞
==≈=00!
!n n N n n x
n x n x e
舍入误差
在数值计算中,通常都按有限位进行运算。例如,按照4舍5入的原则,2 / 3=0.666667 或2 / 3=0.667,由舍入产生的误差,称为舍入误差。
在实际计算中的数据通常是近似值,它们由观察、估计或一些计算而得到,这些数在计算机表示后也会带来进一步误差,即误差的积累和传播。关于误差的传播似乎没有多少统一的理论,通常积累误差的界是以通例分析为基础而建立的。
有效位数
定义:当x的误差限为某一位的半个单位,则这一位到第一个非零位的位数称为x的有效位数。
例如,x= 12.34 ,y= 0.004067均有4位有效数字,
3.00与3.0000分别有3位和5位有效位数。
0.3约束误差
选择收敛的稳定的方法
对同一问题选择不同的数值计算方法,可能得到不同的计算结果。在计算方法中,除了给出方法的数值计算公式,还要讨论计算公式的收敛性、稳定性和截断误差的特性。选择收敛性要求低、稳定性好的方法是约束误差扩张最重要的措施。例如,样条插值函数比高次多项式的效果好的多,是构造插值函数的首选方法。
提高数值计算精度(用优质计算机)
数值在计算机中存放的位数称为字长。有限位的字长是带来舍入误差和抑制数值计算精度的根源。对同一种方法,在字长大的计算机上的计算效果要比在字长小的计算机上优越。
在计算机上,用同一种数值计算方法对数据选用不同的数值类型,有时会直接影响到计算效果。例如,对病态的线性方程组,采用单精度数据使用消元方法,其数据解大大失真,而用双精度数据有时却可得到满意的数值解。
第6章解线性方程组
在数值计算和工程应用中,求解线性代数方程都具有中心的作用。在石油勘探,天气预报等问题中常常出现成百上千阶的方程组,也就产生了各种形式方程组数值解法的需求。研究大型方程组的解是目前计算数学中的一个重要方向和课题。
解方程组的方法可归纳为直接解法和迭代解法。从理论上来说,直接法经过有限次四则运算,假定每一步运算过程中没有舍入误差,那么,最后得到方程组的解就是精确解。但是,在计算过程中,完全杜绝舍入误差是不可能的,只能控制和约束舍入误差的增长和危害,这样直接法得到的解也不一定是绝对精确的。
迭代法是将方程组的解看作某种极限过程的向量极限的值。在用迭代算法时,我们不可能将极限过程算到底,只能将迭代进行有限多次,得到满足一定精度的方程组的近似解。
在数值计算历史上,直接解法和迭代解法交替生辉。一种解法的兴旺与计算机的硬件环境
和问题规模是密切相关的。一般说来,对同等规模的线性方程组,直接法对计算机的要求高于迭代法。对于中等规模的线性方程组(方程组200 未知量x x x n 12,,, 的数值, a x a x a x b a x a x a x b a x a x a x b n n n n n n nn n n 1111221121122222 1122+++=+++=+++=?????? ? (6.1) 其中a b ij i ,为常数。 写成矩阵形式 Ax b = 称 A a ij =()为系数矩阵,x x x x n T =(,,,)12 为解向量,b b b b n T =(,,,)12 为常数向量。 §6.1高斯(Gauss )列主元消元 §6.1.1 高斯(Gauss )消元法 高斯消元法是我们熟悉的古老、简单而有效的解方程组的方法。下列中学阶段解二元方程组(高斯消元法)的步骤: x x x x 121211328 2-=+=?? ?()() 方程(1)乘以 -3加到第(2)个方程,得到方程(3) 1221(1)055 (3) x x x -=?? +=? ∴=x 21 代入(1) 得x 12= 在线性代数中,对方程组的增广矩阵做行的初等变换 (,)(~,~ )()()A b A b =-?? ????→???-?? ?? ?=-?+131218101515312 ~ A 已是上三角阵, x x x 1221 55-==??? 为原方程组的等价方程组,再用回代方法,得到: 这就是高斯消元法解方程组的消元和回代过程。 在解方程组中施行以下三种初等变换: 对换某两个方程的次序; 对其中某个方程的两边同乘一个不为零的数; 把某一个方程两边同乘一个常数后加到另一个方程的两边。 我们以n = 4为例演示高斯消元过程。 方程组 AX b =,其增广矩阵为 (,)A b a a a a b a a a a b a a a a b a a a a b =?? ?? ?? ??11121314121 222324231 323334341 42 43 444 k=1,若 a 110≠, 将第一行乘以-a a 2111/加到第二行上; 将第一行乘以-a a 3111/加到第三行上; 将第一行乘以-a a 4111/加到第四行上, 即 11111111 24{,}i i i i i i a a for i to A A A b b b a a =- +→-+→ 得到 11 12 13 14 1(1)(1)(1)(1)2223332 (1)(1)(1)(1)(1)3233343(1)(1)(1)(1)42 43 44 4000 a a a a b a a a b A a a a b a a a b ?? ? ? = ? ??? 即 ) 1() 1(b X A =, k=2,若 a 2210() ≠ 将第二行乘以-a a 321221() () /加到第三行上;将第二行乘以-a a 421221() () /加到第四行上; 即 22222222 34{,}i i i i i i a a for i to A A A b b b a a =- +→-+→ 得到 11 12 13 141(1)(1)(1)(1)22 23242 (2) (2)(2)(2)33343(2)(2)(2)43 4440000 a a a a b a a a b A a a b a a b ?? ? ?= ? ??? k=3,若 a 3320() ≠ 将第三行乘以-a a 432332()() /加到第四行上; 即 33333333 44{,}i i i i i i a a for i to A A A b b b a a =- +→-+→ 得到 111213141(1)(1)(1)(1)2223242(2)(2) (2)33343(3)(3)4440(,)00000a a a a b a a a b A b a a b a b ?? ? ? = ? ??? (10) ~ A 已是上三角阵,于是得到了与原方程等价的上三角形式的方程组: a x a x a x a x b a x a x a x b a x a x b a x b 1111221331441221223132414213323342432443443+++=++=+==???? ? ? ?()()() ()()()() ()() (11) 再对上方程组(.11)依次回代解出x x x x 4321,,, 。 可以得到系数矩阵A 的行列式的值为~ A 的对角元素的乘积。即 Det A a a a a ()() () () =11221332443 这也正是在计算机上计算方阵A 的行列式的常规方法。 高斯消元法算法 在上述描述中用a ij k () 的上标是为了帮助理解在消元过程系数矩阵元素的变换。由于在消元过程中,将a ij k () 仍放在a ij k ()-1元素的位置上,在下列算法中将略去a ij k () 的上标。 13 14 {,} ik ik k i i k i i kk kk for k to for i k to a a A A A b b b a a ==+-+→-+→ 再对ik k i i kk a A A A a - +→细化为: 14 ik kj ij ij kk a for j k to a a a a =+- +→ 要将上列解方程步骤推广到n 阶方程组,只需将对控制量“4”的操作改成对控制量n 的操 作即可。 经过以上消元,我们已得到与b Ax =的等价方程组b x A =,其中A 已是一个上三角矩阵。为简单起见,记A U =的元素为ij u ,b ~ 的元素为i b ,回代 1,,1,,/)(1 -=- =∑+=n n i a x a b x n i j ii j ij i i 即已得到原方程组的解。 例6.1 求解上三角形方程组 u x u x u b u x u x b u x b n n n nn n n 11112211 22222+++=++==?? ?? ? ? ? 解 由第n 个方程 u x b nn n n = , 得 x b u n n nn =, 将n x 的值代入到第n -1方程 x b u x u n n n n n n n -----=-11111(),, 将11,,,+-i n n x x x 的值代入到第i 个方程 u x u x u x b ii i i i i in n i +++=++,11 得解的通式 x b u x u i n n i i ij j j i n ii =- =-=+∑(),,,,,1 121 计算x i 需要n i -+1次乘法或除法运算,i n n =-,,,,121 。因此,求解过程中的运算量为 ()() ()n i j n n O n j n i n -+== +===∑∑112 1 12 。 高斯消元法的基本思想就是通过对方程组做初等变换,把一般形式的方程组化为等价的上 (下)三角形或对角形方程组。 高斯消元算法 在算法中略去了变量、矩阵和向量的定义部分。在消元过程中,将a ij k () 仍放在a ij k () -1元素 的位置上。 step1 输入:方程组阶数n 、方程组系数矩阵A 和常数向量B step 2 for k = 1 to n -1 !消元过程 { for i = k+1 to n { t a a i k k k :/= for j = k+1 to n { a a t a i j i j k j :=-? } !ENDFORJ b b t b i i k :=-? } !ENDFORI } !ENDFORI step 3 for i = n to 1 ! 回代求解 ii n i j j ij i i a x a b x )(1 ∑+=- = { s b i = for j = i +1 to n s s a x i j j =-? x s a i i i :/= } step 4 输出方程组的解 n i x i ,,2,1, =. 高斯(Gauss )消元法的运算量 整个消元过程的乘法和除法的运算量为 ()()n k n k k n --+=-∑11 1+()n k k n -=-∑1 1 , 回代过程的乘除运算量为() n n +12 , 故Gauss 消元法的运算量为 ()n n n O n 32 3313 +-= §6.1.2 列主元消元法 在上面的消元法中,未知量是按照在方程组中的自然顺序消去的,也叫顺序消元法。 在消元过程中假定对角元素a kk k () ,-≠10 消元步骤才能顺利进行,由于顺序消元不改变A 的主子式值,故顺序高斯消元法可行的充分必要条件为A 的各阶主子式不为零。但是,实际上只要detA ≠0,方程组Ax = b 就有解。故顺序高斯消元法本身具有局限性。 另一方面,即使高斯消元法可行,如果a kk k () -1很小,在运算中用它作为除法的分母.会导致 其它元素数量级的严重增长和舍入误差扩散。这是高斯消元法的另一缺陷。例如: 例6.2 方程组000033000020001110000 100001000021212...() ...() x x x x +=+=?? ? (*) 的精确解为: x x 12132 3 ==,. 在高斯消元法计算中保留4位小数。 解: 方程()()/.()11000032?-+方程得: 0000330000 200019999066660122.....x x x +==?? ? ∴=x 206667.,代入方程(1)得: x 10= 由此得到的解完全失真,如果交换两个方程的顺序,得到等价方程组 10000 10000100000000330000200011212......x x x x +=+=?? ? 经高斯消元后有 10000100001000029997199980666703333 12221......,.x x x x x +==?? ?∴== 由此可看到,在有些情况下,调换方程组的次序,对方程组的解是有影响的。 如果不调换方程组的次序,取6位有效数字计算方程组(*)的解, 得到 33.0,666667 .012==x x ; 取9位有效数字计算方程组(*)的解, 得到 333333 .0,666667.012==x x ;由此可见有效数字的作用。 调换方程组的次序是依照在运算中做分母量的绝对值尽量的大,以减少舍入误差的影响。如果在一列中选取按模最大的元素,将其调到主干方程位置再做消元,则称为列主元消元法。 用列主元方法可以克服高斯消元法的额外限制,只要方程组有解,列主元消元法就能 畅通无阻地顺利求解,它在在消元法中有效的抑制了舍入误差的增长,同时又提高了解的精确度。 更具体地,第一步在第一列元素中选出绝对值最大的元素max{||}||111≤≤=i n i m a a ,交换第 一行和第m 行的所有元素,再做化简a a a n 21311,,, 为零的操作。 对于每个1,,2,1,-=n k k ; 在实现消元前,选出{||,||,,||}() ,()() a a a kk k k k k n k k -+--1111 中绝对值最大的元素(1) ,k m k a -,对k 行和m 行交换后,再做消元操作,这就是列主元消元法的 操作步骤。由于det A ≠0,可证a a a kk k k k k n k k () ,()() ,,,-+--1111 中至少有一个元素不为零,因此,列主元消元总是可行的。列主元消元法与高斯消元法相比,只增加了选列主元和交换两个方程组(即两行元素)的过程。 如果对于第k 步,从k 行至n 行和从k 列至n 列中选取按摸最大的||a i j ,k ≤≤i n , k ≤≤j n ,记||max () ,a uv k k i j n -≤≤=1{||() a i j k -1},对第k 行和第u 行交换,对第k 列和第v 列交换, 这就是全主元消元法。在k 列和第v 列交换时,还要记录下v 的序号,以便恢复未知量x k 和x v 的位置。 高斯消元法将系数矩阵化为上三角阵,再进行回代求解;Guass-Jordan 消元法是将系上三角阵继续做初等变换化为对角阵,再进行求解, (6)11 12 13 14 111 1(1)(1)(1)(1)(1) (5)22 23242 22 2 (2)(2)(2) (2)(4)33 34333 3(3)(3)(3)(3)44 44440 00 00a a a a b a b a a a b a b a a b a b a b a b ???? ? ? ? ?→ ? ? ? ? ?? ?? ? 那么 x b a i n i i ii ==/,,,,.12 用初等变换化系数矩阵为对角矩阵的方法称为Gauss-Jordan 消元法。 例6.1.2 解方程组与计算逆矩阵 在线性代数中逆矩阵是按其伴随矩阵定义的,若||A ≠0,则方阵A 可逆,且 A A A -=1 1|| * ,其中A * 为A 的伴随矩阵。要计算2 n 个1-n 阶的行列式才能得到一个伴随矩阵,在计算中因其计算工作量大而不被采用。通常对(,)A I 做行的初等变换,在将A 化成I 的过程中得到(,)I A -1 。在数值计算中,这仍然是一种行之有效的方法。 由逆矩阵的定义 A A I ?=-1 令 A a A x ij ij ==-(),()1 ,有 a a a a a a a a a n n n n nn 11121212221 2 ?? ??????x x x x x x x x x n n n n nn 11 121212221 21 000 100 01 ?? ??????= ?? ? ? ? ? ? ? 化为求解n 个方程组 j nj j j nn n n n n e x x x a a a a a a a a a =???? ?? ? ??=??????? ????????? ??010212 1 22221 11211 ,j n =12,,, j e 是第j 个分量为1,其余分量为零的n 维向量。 或记为: Ax e j n j j ==, ,,,.12 用直接法或迭代法算出 x j n j ,,,,=12 ,也就完成了逆矩阵A -1 的计算。 如果依次对(,),(,),,(,)A e A e A e n 12 作Gauss-Jordan 消元法,组合起来看有(当然也能组合起来做): (,,,,)(,)A e e e A I n 12 = 这正是在线性代数中的初等变换计算逆矩阵的方法。 由此可见,计算一个n 阶逆矩阵的工作量相当于解n 个线性方程组。在数值计算中,常常将计算矩阵逆的问题转化为解线性方程组的问题。 §6.2直接分解法 在高斯消元法中,每一个消元步骤即对方程组做一次初等变换,从矩阵乘法的角度看,每做一次初等变换相当于用一个初等矩阵左乘系数矩阵A 和常数项b 。 仍以n=4为例,k=1,记 1 111i i a a =- i =234,, 容易验证, 111213142122232421313233343141 42 43 44411 00010001000010001001000010010001000100010001a a a a a a a a l a a a a l a a a a ???? ?? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???????? 2131411000100010001?? ? ?= ? ??????? ?? ? ??4443 42 41 34333231 242322 2114131211 a a a a a a a a a a a a a a a a ===) 1(1A A T ???? ?? ? ? ?)1(44 )1(43 )1(42 )1(34)1(33)1(32)1(24)1(23)1(2214 13 12 11 000a a a a a a a a a a a a a (5.8) 由 T A x T b 11= ,得到 A x b ()() 11= ,其中 ?????? ? ??==)1(4)1(3)1(21) 1(1b b b b b b T 同理记 (1)(1) 2 222i i a a =- i =34, 11 12 13 14(1) (1)(1)(1)(1)(2)22 2324 2(2)(2)32 3334 (2)(2)42 43 44100001000100 01a a a a a a a T A A A l a a a a ???? ? ? ? ?=?= = ? ? ? ? ???? ? 记 () 4343233 2=-a () 11 12 13 141112 13 14(1) (1)(1)(1) (1)(1)22 2324 (2) 22 2324 3(2)(2)(2)(2)33 343334 (3)(2)(2)43 4443 441 0000100001 00 1a a a a a a a a a a a a a a T A A a a a a a a a ?????? ? ? ? ? ? ?=?== ? ? ? ? ? ? ?? ? ??? ? ??????? ??==)2(4)2(3 )1(21) 2()1(2b b b b b b T ,1(1)2 (2)(3) 3(2)3(3)4b b T b b b b ?? ? ?== ? ??? 即所有的消元步骤表示为:(,)A b 左乘一系列下三角初等矩阵;容易验证,这些下三角阵的乘积仍为下三角矩阵 T T T T ==3212131324142 43 10001 00101l l l l l l ?? ? ? ? ??? , T T T A A 321=~ 于是有 ~A TA = 或 A T A =-1~ 这里 T -1 仍为下三角阵,其对角元素为1,称为单位下三角阵。而~ A 已是上三角阵 记 T -1 =L ,~ A =U ,则有 A LU = 结果表明若消元过程可行,可以将A 分解为单位下三角L 与上三角阵U 的乘积。由此派生出解方程组的直接分解法。 §6.2.1 Dolittle 分解 由高斯消元法得到启发,对A 消元的过程相当于将A 分解为一个上三角阵和一个下三角阵的过程。如果直接分解A 得到L 和U ,A=LU 。这时方程AX=b 化为LUX=b , 令UX=Y ,由LY=b 解出Y ;再由LY=B ,解出X 。这就是直接分解法。 将方阵A 分解为A=LU ,当L 是单位下三角阵,U 是上三角阵时,称为Dolittle 分解;当L 是下三角阵,U 是单位上三角阵时,称为Courant 分解。分解的条件是若方阵A 的各阶顺序主子式不为零,则Dolittle 分解或Courant 是可行的。 Dolittle 分解步骤 在矩阵A = LU 的Courant 分解形式中,L 是下三角矩阵,U 是单位上三角矩阵, a a a a a a a a a a n n n n n n 11121212221 2 3 4 ?? ??????= 1112112l l l n n ?? ? ? ?? ??u u u u u u n n nn 11121222 ?? ? ? ? ? ? ? (5.16) 矩阵L 和U 共有n 2 个未知元素,按照一定的顺序,对每个a ij 列出(5 .16)两边对应的矩阵乘法关系,一个a ij 关系式解出一个L 或U 的元素。按下列步骤计算L 和U 的元素 计算U 的第一行元素n u u u 11211,,, 要计算 u 11,则列出(5.16)等号两边的第1行第1列元素的关系式 a 11=l u r r r n 111 =∑=() 100 u 1100 ?? ?? ?? ??=u 11, 故u 11=a 11, 一般地,由U 的第一行元素的关系式 a j 1=l u r r n rj 11 =∑=()100012 u u j j ?? ?? ????= u j 1 得到 ∴==u a j n j j 1112,,, ? 计算L 的第一列元素 13121,,,n l l l 要计算l 21,则列出(5.16)等号两边的第2行第1列的元素 a 21= ()l u l u l u r r r n 211 21 1121111000 0=∑=?? ?? ????= 故 112121/u a l = 一般地,由L 的第一列元素的关系式 ( ) a l u l l u l u i ir r r n i i i i 111 11 1111110000= =?? ?? ????==-∑ , 得到 ∴==l a u i n i i 1111 23/,,, , ? 计算U 的第二行元素 ? 计算L 的第二列元素 若已算出U 的前k -1行,L 的k -1列元素 ? 计算U 的第k 行元素 u u u kk k k k n ,,,,,+1 由U 的第k 行元素的关系式 ()????? ??? ? ??== -=∑0000111 ,211 jj j k k k k n r rj kr kj u u l l l u l a u 是上三角阵 ∴ (行标)k j ≤(列标) ∑∑∑-===+=== 1 1 1 1 k r kj rj kr k r rj kr n r rj kr kj u u l u l u l a 得到 ∴=- =+=-∑u a l u j k k n kj kj kr rj r k 1 1 1,,,, (5.17) 计算L 的第k 列元素l l l k k k k n k ++12,,,,,, 由L 的第k 列元素的关系式 a l u ik ir r n rk ===∑1(,,,,,,),l l u u i i i k kk 11110000 -?? ?? ?? ? ???? L 是下三角阵 ∴ (行标)i k ≥(列标) kk ik k r rk ir rk ir k r rk ir n r ik u l u l u l u l a +=== ∑∑∑ -===11 1 1 得到 ∴=- =-∑l a l u u ik ik ir rk r k kk (),1 1 i=k+1,k+2, ,n (5.18) 一直做到L 的第n -1列,U 的第n 行为止。 用LU 直接分解方法求解方程组所需要的计算量仍为13 032 n n +(),与用Gauss 消元 法的计算量基本相当。 可以看到在分解中A 的每个元素只在(5.17)或(5.18)中做而且仅做一次贡献,如果需要节省空间,可将U 以及L 的元素直接放在矩阵A 相应元素的位置上。 用直接分解法解方程Ax b =。首先作出分解A = LU ,令UX=Y ,解方程组LY b =,由于L 是单位下三角阵,容易得到 y b l x i n i i ij j i j =- ==-∑11 12,,,, (5.19) 再解方程组UX=Y ,其中 x y u x u i n n i i ij j j i n ii =- =-=+∑(),,,,,1 121 (5.20) 对A 进行LU 分解时,并不涉及常数项b 。因此,若需要解具有相同系数的一系列线性方程组时,使用直接分解法可以达到事半功倍的效果。 Dolittle 直接分解算法 例6.2.1用Doolittle 分解求解方程组 ??? ??=++=++=++5 2262342321 321321x x x x x x x x x 解: ???? ? ??????? ??=????? ??332322 1312 1132 3121000101 00 1221231112u u u u u u l l l 对 k =1,u a j i j 11123==,,, u u u 111213211===,, l a u i i i 111123==/,, l l 213112051205====/.,/. 对k=2 , u a l u 22222112305 25=-=-=.. u a l u 23232113205 15=-=-=.. ()6.0)5.02(5 .211 1231322232=-=-=u l a u l 对k=3 u a l u l u 33333113322306=--=. ∴ LU =10005100506121102515 0006......?? ??????? ?? ? ?? 解LY b = 100051005061465123...?? ??????? ?????=?? ?? ? ??y y y y y y 1234406===,,. 解UX = Y 21102515 00064406123....?? ??????? ?????=?? ?? ? ??x x x ∴===x x x 321111,, §6.2.2 Courant 分解 Courant 分解步骤 在矩阵A = LU 的Courant 分解形式中,L 是下三角矩阵,U 是单位上三角矩阵, 记 L =l l l l l l n n nn 11212212 ?? ? ?? ?? ? ,U u u u n n =?? ? ? ? ? ? ?1111212 (5.21) 矩阵 Courant 分解的次序与Doolittle 分解的次序不同,设A 的各阶主子式不为零,按计算L 的第一列,U 的第一行, ,L 的第k 列,U 的第k 行的顺序进行计算,类似于Doolittle 分解的手法,作出A=LU 两边矩阵元素的关系式 a l u i j i r r n r k =-∑1 ,可导出计算 公式 对k =1,2, ,n l a l u i k k n ik ik ir rk r k =- =+=-∑,,,,11 1 (5.22) u a l u l j k k n kj kj kr rj r k kk =- =++=-∑(), ,,,1 1 12 (5.23) 如果A 的各阶主子式不为零,则(5.23)中的l kk ≠0,Courant 分解才能顺利进行。在线性代中我们已知,系数矩阵A 的行列式不为零,方程组的解存在而且唯一,并没有对主子式有任何要求。在直接分解中也能采用与列主元消元法类似的方法进行处理,避免A 的某阶主子式为零或|l kk |很小。即每求一列L 的元素后,求出该列中按模最大的元素l uk ,然后将矩阵A 和L 的第k 行和第u 行元素逐个交换,并交换常数项b k 与b u 。 对于方程组AX b =,作Courant 分解后,方程组可写成LUX b =,仍记UX y = 由Ly b =,得到 y b l y l i n i i i j j j i i i =- ==-∑(),,,,1 1 12 (5.24) 由UX y =,得到 x y u x i i j i n ij j =-=+∑ 1 , j n n =-,,,,121 (5.25) 例6.2.2用Courant 分解求解方程组 121215016246350123----?? ??????? ?????=-?? ?????x x x 和 12121501611117123----?? ??????? ?????=--?? ?? ???x x x 解: 12121501600010100 11121 2231 32 33121323----?? ?????=?? ??????? ?? ???l l l l l l u u u k=1 ; 计算L 的第一列 l a i i i 11123==,,, ; l l l 112131120==-=,, 计算U 的第一行 u a l 121211==/2,u a l 131311==/ 1 k=2 ; 计算L 的第二列 l a l u 22222112=-= 3 ; l a l u 323231121=-=- ; 计算U 的第二行 u a l u l 2323211222=-()/= -1; k=3 ; l a l u l u 3333311332225=--=); L U =--?? ?????=-?? ?? ???100230015121011001 分别解方程: Ly b UX y 1111==, 和 Ly b UX y 2222==, ???? ? ??-=??? ?? ??-=????? ??--=9524, 506324510032001111y y Ly ; ???? ? ??==??? ?? ??-=947, 1001101211111X y X UX ????? ??--=5100320012Ly ????? ??-=??? ?? ??--=2711, 1711122y y ; ???? ? ??-==??? ? ? ??-=253, 1001101212222X y X UX 例6.1.3 推导三对角矩阵A LU =分解公式,求解方程组Ax f = 1111 22 222 1111 111n n n n n n n n a b u v c a b w u A v a b c w u c a ----?? ???? ? ?? ? ? ??? ?== ??? ? ??? ??? ?? ?? ? 解: 比较A L U =?两边元素,可得到 1 1111/,,,3,2,u b v a u n i c w i i ==== 2221222/,u b v v c a u =-= 若规定 c 10= ,可得到一般计算公式 n i u b v v c a u i i i i i i i ,,2,1,/1 =?? ?=-=- (5.27) 考虑三对角线系数矩阵的线性方程组 Ax f = 这里 x x x x f f f f n T n T ==(,,,),(,,,)1212 令 Ux y =,则有 Ly f = 于是有 n i u y c f y i i i i i ,,2,1, /)(1 =-=- (5.28) 规定0=n v ,由Ux y = 可得到 1,2,,1,, 1 -=-=+n n i x v y x i i i i 事实上,由(5.27),每计算出一个i u ,就可按(5.28)计算出y i 。若规定c 10=,由(5.27)和(5.28)可合并成 ??? ?? =-==-=--n i u y c f y u b v v c a u i i i i i i i i i i i i ,,2,1/)(/11 (5.29) 1,2,,1,, 1 -=-=+n n k x v y x k k k k (5.30) 称计算公式(5.29)和(5.30)为追赶法。 §6.2.3 对称正定矩阵的LDL T 分解 对称正定矩阵也是很多物理问题产生的一类矩阵,正定矩阵的各阶主子式大于零。由线性代数中的定理,若A 正定,则存在下三角矩阵U ,使A UU T =,直接分解A UU T =的分解方法,称为平方根法。对于对称正定矩阵,常用的是LDL T 分解。 对A 作Doolittle 分解A=LU A =1112112l l l n n ?? ?? ????u u u u u u n n nn 1112122 2 ?? ? ?? ? ? ? / 提出U 矩阵的对角元素 / =11121 1 2 n n ?? ?? ????u u u nn 1122 ?? ?? ????1111212u u u n n ?? ? ? ? ? ?? 由A 对称正定,可得u ii >0,令 D u u u d d d nn n =?? ??????=?? ?? ????1122 12 可证 1111212u u u L n n T ?? ?? ? ??? = , 即