文档库 最新最全的文档下载
当前位置:文档库 › 第六讲解线性方程组

第六讲解线性方程组

绪 论

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

?? ??

?

???

= , 即

相关文档
相关文档 最新文档