数值分析讲义--线性方程组的解法

第三章 线性方程组的解法

3.0 引 言

重要性:解线性代数方程组的有效方法在计算数学和科学计算中具有特

殊的地位和作用. 如弹性力学、电路分析、热传导和振动、以及社会科学及定量分析商业经济中的各种问题.

分类:线性方程组的解法可分为直接法和迭代法两种方法.

(a) 直接法:对于给定的方程组,在没有舍入误差的假设下,能在预定的运算次数内求得精确解. 最基本的直接法是Gauss 消去法,重要的直接法全都受到Gauss 消去法的启发. 计算代价高.

(b ) 迭代法:基于一定的递推格式,产生逼近方程组精确解的近似序列. 收敛性是其为迭代法的前提,此外,存在收敛速度与误差估计问题. 简单实用, 诱人.

(AX =b )

1

基本思想:

与解f (x )=0 的不动点迭代相类似, 将AX =b 改写为X =BX +f 的形式, 建立雅可比方法的迭代格式:X k +1=BX (k ) +f , 其中,B 称为迭代矩阵. 其计算精度可控, 特别适用于求解系数为大型稀疏矩阵(sparse matrices)的方程组. 2

问题:

(a) 如何建立迭代格式?

(b) 向量序列{X k }是否收敛以及收敛条件? 3 例题分析:

⎧10x 1-x 2-2x 3=7. 2

考虑解方程组⎨-x 1+10x 2-2x 3=8. 3 (1)

⎪-x -x +5x =4. 2⎩123

其准确解为X *={1, 1.2, 1.3}. 建立与式(1)相等价的形式:

⎧x 1=0. 1x 2+0. 2x 3+0. 72

⎨x 2=0. 1x 1+0. 2x 3+0. 83 (2) ⎪x =0. 1x +0. 2x +0. 84⎩312

据此建立迭代公式:

(k ) k

⎧x 1(k +1) =0. 1x 2+0. 2x 3+0. 72⎪(k +1) (k ) (k ) x =0. 1x +0. 2x +0. 83 (3) ⎨213

⎪x (k +1) =0. 1x (k ) +0. 2x (k ) +0. 84

12⎩3

取迭代初值x

(0)

1

=x

(0) 2

=x

(0) 3

=0, 迭代结果如下表.

x1 x2 x3

0 0 0 0 1 0.72 0.83 0.84 2 0.971 1.07 1.15 3 1.057 1.1571 1.2482 4 1.08535 1.18534 1.28282 5 1.095098 1.195099 1.294138 6 1.098338 1.198337 1.298039 7 1.099442 1.199442 1.299335 8 1.099811 1.199811 1.299777 9 1.099936 1.199936 1.299924 10 1.099979 1.199979 1.299975 11 1.099993 1.199993 1.299991 12 1.099998 1.199998 1.299997 13 1.099999 1.199999 1.299999 14 1.1 1.2 1.3 15 1.1 1.2 1.3

4 Jocobi 迭代公式:

设方程组AX=b, 通过分离变量的过程建立 Jocobi 迭代公式,即

i =1

∑a ij x j =b i , a ii ≠0(i =1, 2, , n )

n 1

x i =(b i -∑a ij x j ) (i =1, 2, , n )

j =1a ii

j ≠i

n

由此我们可以得到Jacobi 迭代公式:

x

(k +1)

i

n 1

=(b i -∑a ij x i k ) (i =1, 2, , n )

j =1a ii

j ≠i

[Jacobi迭代公式的算法] 1: 初始化. n , (a ij ), (b j ), (x 1) , M . 2: 执行k =1直到M 为止. ①

执行i =1直到n 为止.

u i ←(b i -∑a ij x j ) /a ii ;

j =1j ≠i

n

② ③

执行i =1直到n 为止.

x i ←u i ;

输出k , (x i ).

另外,我们也可以建立Jacobi 迭代公式的矩阵形式. 设方程组AX =b ,

其中,A =(a ij ) n 为非奇异阵,

X =(x 1, x 2,…,x n ) T , b =(b 1, b 2,…,b n ) T

将系数阵A 分解为: A =U +D +L ,

U 为上三角矩阵,D 为对角矩阵,L 为下三角矩阵.

AX =b 可改写为 (U +D +L ) X =b

⇔ X =D -1b -D -1(U +L ) X

由此可得矩阵形式的Jocobi 迭代公式: X k +1=BX (k ) +f

§3.2

注意到利用Jocobi 迭代公式计算x i

(k +1)

时,已经计算好

(k ) k ) x 1(k ) , x 2, , x i (-1的值,而Jocobi 迭代公式并不利用这些最新的近似值计

算,仍用x 1

(k )

(k ) k )

, x 2, , x i (-1. 这启发我们可以对其加以改进, 即在每个分量的

计算中尽量利用最新的迭代值,得到

x

(k +1)

i

i -1n 1(k +1)

=(b i -∑a ij x j -∑a ij x k ) (i =1, 2, , n ) j

j =1j =i +1a ii

上式称为Gauss-Seidel 迭代法. 其矩阵形式是

X =-(D +L ) -1UX +(D +L ) -1b , X k +1=BX (k ) +f .

迭代次数 x1 x2 x3 0 0 0 0 1 0.72 0.902 1.1644 2 1.04308 1.167188 1.282054 3 1.09313 1.195724 1.297771

4 1.099126 1.199467 1.299719 5 1.09989 1.199933 1.299965 6 1.099986 1.199992 1.299996 7 1.099998 1.199999 1.299999 8 1.1 1.2 1.3

.3

1

基本思想:

逐次超松弛迭代法(Successive Over Relaxation Method,简写为SOR) 可以看作带参数ω的高斯-塞德尔迭代法,是G-S 方法的一种修正或加速. 是求解大型稀疏矩阵方程组的有效方法之一. 2 SOR 算法的构造:

设方程组AX =b , 其中,A =(a ij ) n 为非奇异阵,X =(x 1, x 2,…,x n ) T , b =(b 1, b 2,…,b n ) T . 假设已算出x (k ) ,

(k +1) i

i -1n 1(k +1)

=(b i -∑a ij x j -∑a ij x k j ) (i =1, 2, , n ) (1)

j =1j =i +1a ii

相当于用高斯-塞德尔方法计算一个分量的公式. 若对某个参数ω,作(k +1)

与i

x i (k ) 加权的平均, 即

(k +1)

i

x

k +1i

=(1-ω) x

(k ) i

+ω=x i (k ) +ω((k +1) i

-x i (k ) )

(2)

其中,ω称为松弛因子.

用(1)式代入(2)式,就得到解方程组AX =b 的逐次超松弛迭代公式:

⎧x i (k +1) =x i (k ) +∆x i ⎪

i -1n ⎨ (3) ω(k +1) (k )

a ij x j -∑a ij x j ) (i =1, 2, , n ) ⎪∆x i =a (b i -∑j =1j =i

⎩ii

显然,当取ω=1时,式(3)就是高斯-塞德尔迭代公式. 3 例题分析:

SOR 方法解方程组

⎧4x 1-2x 2-x 3=0⎪

⎨-2x 1+4x 2-2x 3=-2 (1) ⎪-x -2x +3x =3⎩123

其准确解为X *={1, 1, 2}. 建立与式(1)相等价的形式:

⎪x 1=0. 5x 2+0. 25x 3⎪

⎨x 2=0. 5x 1+0. 5x 3-0. 5 (2) ⎪12

⎪x 3=x 1+x 2+1⎩33

据此建立迭代公式:

⎧(k +1)

(k ) k

=0. 5x 2+0. 25x 3⎪x 1

⎪(k +1) (k )

⎨x 2=0. 5x 1(k ) +0. 5x 3-0. 5 (3)

⎪1(k ) 2(k ) (k +1)

⎪x 3=x 1+x 2+1⎩33

(0) (0) (0)

利用SOR 算法,取迭代初值x 1=x 2=x 3=1,

ω=1.5,迭代结果如下表.

逐次超松弛迭代法

次数 x1 x2 x3 1 0.625000 0.062500 1.750000 2 0.390625 0.882813 1.468750 3 1.017578 0.516602 1.808594

4 0.556885 0.880981 1.710449 5 1.023712 0.743423 1.868103 6 0.746250 0.908419 1.838737 7 0.997715 0.860264 1.913894 8 0.864050 0.936742 1.908605 9 0.986259 0.922225 1.945523 10 0.928110 0.958649 1.947493 11 0.985242 0.955944 1.966198 12 0.961661 0.973818 1.969521 13 0.988103 0.974699 1.979289 14 0.979206 0.983746 1.982172 15 0.991521 0.985318 1.987416 16 0.988509 0.990038 1.989513 17 0.994341 0.991414 1.992397 18 0.993538 0.993946 1.993806 19 0.996367 0.994950 1.995424 20 0.996313 0.996342 1.996331 21 0.997724 0.997018 1.997254 22 0.997871 0.997798 1.997822 23 0.998596 0.998234 1.998355

GS 迭代法须迭代85次得到准确值X *={1, 1, 2};而SOR 方 法只须55次即得准确值. 由此可见,适当地选择松弛因子

ω,SOR 法具有明显的加速收敛效果. □

3.4

1. 向量和矩阵范数 (a) 向量范数

R n 空间的向量范数 || · || ,对任意

, ∈R n , 满足下列条件:

(1) ||x ||≥0; ||x ||=0⇔x =0 (正定性)

(2) ||αx ||=|α|⋅||x || (齐次性)

(3) ||x +y ||≤||x ||+||y || (三角不等式)

常见的向量范数有: (1) 列范数:

(2) 谱范数:(欧几里德范数或向量的长度, 模

)

(3) 行范数:

(4) p 范数:

:

||x ||∞=max(|x 2-x 1|,|y 2-y 1|) ; ||x ||1=|x 2-x 1|+|y 2-y 1| ;

||x ||2=(x 2-x 1) 2+(y 2-y 1) 2.

}依坐标收敛于向量x * 的充要条件是向量序列

{x (k ) }依范数收敛于向量x *, 即lim ||x (k ) -x *||=0.

k →∞

向量序列{x

(k )

(b) 矩阵范数

|| ,对任意 R m ⨯n 空间的向量范数 || ·

A , B ∈R m ⨯n , 满足下列条件:

(1) ||A ||≥0; ||A ||=0⇔A =0

(2) ||αA ||=|α|⋅||A ||(3) ||A +B ||≤||A ||+||B ||(4) || AB || ≤||A ||||B ||

||A ||∞=max ∑|a ij | (行和范数)

1≤i ≤n

n

j =1

||A ||1=max ∑|a ij | (列和范数)

1≤j ≤n

n

i =1

||A ||2=λm ax (A T A ) (谱范数)

若A 对称,则有

λm ax (A T A ) =λm ax (A 2) .

A ||2=ρ(A ) ,

矩阵A 的谱半径记为||

1≤i ≤n

ρ(A ) =m ax |λi |, 其中λi 为A 的特征根。

2. 迭代法基本定理

设有方程组X =BX +f , 对于任意初始向量X (0)及任意f , 迭代公式

X (k +1)=BX (k ) +f 收敛的充要条件是ρ(B )

对于任意初始向量X (0)及任意f , 迭代公式X (k +1)=BX (k ) +f ,

X (k ) -X *=BX (k -1) +f -(BX *+f )

=B (X (k -1) -X *)

=B 2(X (k -2) -X *) =B (X

k

(0)

-X )

k

*

由此可得, 迭代法收敛的充要条件是B 即,

→0(k →∞) .

ρ(B )

上述定理是线性方程组迭代解法收敛性分析的基本定理,然而由于

ρ(B ) 的计算往往比较困难,尽管有各种办法估计ρ(B ) 的上界,但往往

偏听偏大而不实用,由此导致定理的理论价值胜于实用价值,为满足实际判敛的需要,有如下定理

.

(迭代收敛的充分条件)

设有迭代公式X (k +1)=BX (k ) +f ,如果||B ||

3. 从方程组的系数矩阵A 判断迭代收敛性

实际中要求解的某些线性方程组,其系数矩阵往往具有一些特点,如系数矩阵为对称正定、对角元素占优等. 由这些方程组系数矩阵的特殊性,使得我们可以直接从方程组的系数矩阵A 出发来讨论迭代法的收敛性

.

设A =(a ij ) n ⨯n ∈R

n ⨯n

,满足

|a ii |≥∑|a ij |, i =1, 2, , n

j =1j ≠i

n

且至少有一个i 值,使得

|a ii |>∑|a ij |

j =1j ≠i

n

成立,则称A 为对角占优矩阵;若

|a ii |>∑|a ij |, i =1, 2, , n ,

j =1j ≠i

n

则称A 为严格对角占优矩阵

.

如果A =(a ij ) n ⨯n ∈R

n ⨯n

为严格对角占优矩阵,则对任意的初值

x (0), 解方程组AX =B 的Jacobi 法、Guess-Seidel 迭代法均收敛. HW: 3.1 3.2 3.3(上机实习)

§3.5

1

基本思想:

用高斯消去法求解线性方程组的基本思想是设法消去方程组的系数矩阵A 的主对角线下的元素,而将Ax =b 化为等价的上三角形方程组,然后再通过回代过程便可以获得方程组的解.

这种解线性方程组的方法,常称为高斯消去法(Gaussian Elimination)

.

2 例题分析:

⎧6x 1-2x 2⎪12x -8x ⎪12

利用高斯消去法求解方程组:⎨

⎪3x 1-13x 2⎪⎩-6x 1+4x 2

+2x 3+4x 4=12+6x 3+10x 4=34+9x 3+3x 4=27+x 3-18x 4=-38

2x 2

8x 2 13x 24x 2

+2x 3+4x 4=12+6x 3+10x 4=34+9x 3+3x 4=27+x 3-18x 4=-38

(1)

a i 1

r 1, i =2,3,4.得 利用r i -a 11

⎧6x 1-2x 2⎪ - 4x ⎪2⎨

⎪ ⎪⎩+2x 3+4x 4= 12+2x 3+ 2x 4= 10+8x 3+ x 4= 213x 3-14x 4=-26

(2)

a i (22)

利用r i -(2) r 2, i =3,4.得

a 22

⎧6x 1-2x 2+2x 3+4x 4= 12⎪ - 4x +2x + 2x = 10⎪234

(3) ⎨

⎪ 2x 3- 5x 4= -9⎪⎩ 13x 4=-21

a i (33)

利用r i -(3) r 3, i =4.得

a 33

⎧6x 1-2x 2+2x 3+4x 4= 12⎪ - 4x +2x + 2x = 10⎪234

(4) ⎨

⎪ 2x 3- 5x 4= -9⎪⎩ -3x 4=-3

显然,方程组(4)与(1)是等价的,其系数矩阵为上三角

状的,易于求解. 称以上过程为高斯消去法的消去过程. 通过方程组(4)的回代求解,可以得到准确解为

X *=[1, -3, -2,1]T .

这一过程为高斯消去法的回代过程. 2 高斯消去法算法的构造:

记方程组AX =b 为A (1)X =b (1), 其中,A (1)和b (1)的元素分别记为

(1)

a ij 、b i (1) , i 、j =1, 2, , n .

Step1:第一次消元 设a

(1) 11

1) (1)

倍,(i =2, „, n ), 目≠0,将增广矩阵的第i 行减去m i 1=a i (1/a 11

的是将增广矩阵的第一列内除每一个元素不变外,其余全部消为零,得到A (2)X =b (2), 即

(1)

⎡a 11⎢(1) a 21(1) (1) ⎢[A b ]=⎢ ⎢(1) ⎢⎣a n 1(1) ⎡a 11⎢⎢0⎢ ⎢⎢⎣0

(1) 1)

a 12... a 1(n (1) (1) a 22... a 2n

(1) (1) a n ... a 2nn (1) 1) a 12... a 1(n (2) (2) a 22... a 2n

(2) (2) a n ... a 2nn

b 1(1) ⎤

(1) ⎥b 2⎥⇒ ⎥⎥(1) b n ⎥⎦

b 1(1) ⎤(2) ⎥b 2⎥

=[A (2) b (2) ] ⎥⎥(2) b n ⎥⎦

其中

(2) (1)

⎧a ij =a ij -m i 1a 1(1j )

(i , j =2, ..., n ) ⎨(2) (1) (1)

⎩b i =b i -m i 1b 1

Step2:第k 次消元(2≤k ≤n -1)

(k )

设第k -1次消元已完成,且a kk

≠0,得到A (k ) X =b (k ) , 即

1)

a 1(n

(1) ⎡a 11⎢⎢⎢(k ) (k )

[A b ]=⎢

⎢⎢⎢⎣

(1) a 121)

a 1(k

(2) (2) (2) a 22 a 2 a k 2n

(k ) (k )

a kk a kn

(k ) (k ) a kn a nn

b 1(1) ⎤

(2) ⎥b 2⎥ ⎥ (k ) ⎥b k ⎥ ⎥(k ) ⎥b n ⎦

(k ) (k )

m ik =a ik /a kk

(i =k +1, ..., n ) ,

(k +1) (k ) (k ) ⎧a ij =a ij -m ik a kj ⎨(k +1) (k ) (k ) b =b -m b ⎩i i ik k

(i , j =k +1, ..., n )

如此反复,经过n -1次消元之后得到一个与原方程组等价的上三角形方程组.

(1)

⎡a 11⎢⎢⎢b (n ) ]=⎢

⎢⎢⎢⎣

(1)

a 12(2) a 22

1)

⎤⎡x 1⎤⎡b 1(1) ⎤a 1(n

(2) ⎥⎢⎥⎢b (2) ⎥... a 2x n ⎥⎢2⎥⎢2⎥. ⎥⎢. ⎥⎢. ⎥

=⎢⎥

⎥. . ⎢. ⎥

⎥⎢⎥⎢⎥. ⎥⎢. ⎥⎢. ⎥

(n ) (n )

⎥⎢x ⎢⎥b a nn ⎦⎦⎣n ⎦⎣n ⎥

...

[A (n )

Ste 只要a

回代

(n )

nn

≠0就可以回代求解

n

(n ) (n )

x n =b n /a nn

(i )

b -∑a ij x j (i ) i

x i =

3 高斯消去法[算法] Step1[消元]: 对k =1,2,…,n -1 ① ②

若a kk

(k )

a

j =i +1

(i ) ii

(i =n -1, ..., 1)

=0则停止计算

(k ) (k ) =a ik /a kk ;

对i =k +1,k+2,…,n

计算因子m ik

对j =k +1,k+2,…,n

(k +1) (k ) (k ) ⎧a ij =a ij -m ik a kj 计算⎨(k +1)

(k ) (k )

b =b -m b ⎩i i ik k

;

Step2[回代]: 对i =n , n -1,…,1

(i )

b -∑a ij x j

j =i +1

(i ) ii

(i ) i

n

x i =

a

(高斯消去法的条件)

(1)

若A 的所有顺序主子式均不为0,则高斯消元无需换行即可进行到底,且得到唯一解.

(2)

若消元过程中允许对增广矩阵进行行交换,则方程组Ax =b 可用消去法求解的充要条件是A 可逆.

3.6 高斯主元素消去法

1 主元素及其选取问题

Gauss 消去法第k 次消去是用第k 个方程

(k ) (k ) a kk x k + +a kn x n =b k (k )

来消去第k +1,…,n 个方程中的x k , 条件是a kk 关键元素,称为第k 次消去的主元(素). Gauss 消去法存在的问题是: (1) 顺序消元时一旦产生a kk (2) 此外,即使a kk

(k )

(k )

(k )

(k )

是实现第k 次消元的≠0. a kk

=0(这是经常可能的) ,消元过程则中断;

≠0但绝对值很小时,由于用它作除数,引起在消去过

程中出现数量级及舍入误差急剧增长的系数,而使最后的计算解严重地不可靠. 例:单精度求解方程组

⎧10-9x 1+x 2

+x 2⎩x 1

=1

=2

8 ⎧1

⎪⎪x 1=1-10-9=1. 00... 0100...

其准确解为 ⎨

8 ⎪

⎪⎩x 2=2-x 1=0. 99... 9899...

当利用Gauss 消元法时,

8

9

a 22=1-m 21⨯1=0. 0... 01⨯109-109= -10

b 2=2-m 21⨯1= -10

9

舍入误差

⎡10-911⎤⎡10-911⎤⎢⎥⇒⎢99⎥恶性传播 1120-10-10⎣⎦⎣⎦

⇒x

2=1, x 1=0×

基本思想:

主元素法是对Gauss 消去法的改进. 它全面或局部地选

取绝对值大的元素为主元素,仅对Gauss 消去法的步骤作某些技术性地修改,使之成为一种有效的方法. 从而保证和改善算法的数值稳定性. 2 完全主元素消去法

设方程组AX =b , 其中,A =(a ij ) n 为非奇异阵,X =(x 1, x 2,…,x n ) T , b =(b 1, b 2,…,b n ) T . 经过k -1次选主元消元后,得到下列等价方程组:

(1)

⎡a 11⎢⎢⎢(k ) (k )

[A b ]=⎢

⎢⎢⎢⎣

(1) a 12

1)

a 1(k

1)

a 1(n

(2) (2) (2) a 22 a 2 a k 2n

(k ) (k ) a kk a kn

(k ) a nk

(k )

a nn

b 1(1) ⎤(2) ⎥b 2⎥ ⎥ (k ) ⎥b k ⎥ ⎥(k ) ⎥b n ⎦

① 选主元过程

(k ) (k )

⎡a kk ⎤ a kn ⎢⎥

在矩阵 ⎢⎥中选取绝对值最大的元素为主元素,保证 (k ) (k ) ⎢⎥a a nn ⎦⎣nk

(k ) |a i (k , ) j |=max |a ij |≠0; 从而确定 i k , j k .

k

k

k ≤i , j ≤n

② 行变换和列交换

If i k ≠ k then 交换第 k 行与第i k 行;

j k ≠ k then 交换第 k 列与第j k 列;

值得注意的是,在全主元消去过程中,列交换已改

变了x 各分量的顺序,因此,必须在每次列交换的同时,记录调换后未知数的排列次序. ③ ④ ⑤

消元 回代求解

还原未知数的排列次序

2.1 全主元素Gauss 消去法[算法] Step1[消元]: 对k =1,2,…,n -1 ① ② ③

(k ) (k )

|a |=max |a 选主元 确定 i k , j k ∈{k , , n }, 满足i , j ij |≠0;

k

k

k ≤i , j ≤n

(k )

若i k , j k

a =0; 则停止计算,det A =0.

If i k ≠ k then 交换第 k 行与第i k 行; If j k ≠ k then 交换第 k 列与第j k 列;

④ 消元

对i =k +1,k+2,…,n

计算因子m ik

(k ) (k )

=a ik /a kk ;

对j =k +1,k+2,…,n

(k +1) (k ) (k ) ⎧a ij =a ij -m ik a kj 计算⎨(k +1)

(k ) (k )

b =b -m b ⎩i i ik k

;

Step2[回代]: ①

(n )

若nn

a =0则停止计算,det A =0.

② 对i =n , n -1,…,1

b y i =

(n )

i (n )

-∑a ij x j

n

a

j =i +1

(n ) ii

Step3[还原排列次序]: 对i =1,2,…,n -1

x * := y i

(3) 列主元素消去法

在计算机上实现主元素消去法意味着进行数的比较操作,全选主元素法需要相当多的计算时间,因此常采用局部选主元素的方法.

列主元素消去法依次按列选主元素,只须进行方程行 交换,不产生未知数次序的调换. ① 按列选主元过程 设方程组AX =b 的增广矩阵为

(1) ⎡a 11⎢(1) a 21(1) (1) ⎢ [A b ]=⎢ ⎢(1) ⎣a n 1

(1) a 12(1) a 22 (1) a n 2

1)

a 1(n

(1)

a 2n

(1)

a nn

b 1(1) ⎤⎥b 2(1) ⎥

⎥(1) ⎥b n ⎦

首先在A (1) 中第一列选取绝对值最大的元素为主元素,保证

|a i (1, ) 1|=max |a i (11) |≠0; 从而确定 i k .

k

1≤i ≤n

② 行变换

If i k ≠ 1 then 交换第 1 行与第i k 行;

重复上述过程,设已完成第k -1次按列选主元消元后,得到下列等价方程组:

(1) ⎡a 11⎢⎢⎢(k ) (k )

[A b ]=⎢

⎢⎢⎢⎣

(1) a 12

1)

a 1(k

1)

a 1(n

(2) (2) (2) a 22 a 2 a k 2n

(k ) (k ) a kk a kn

(k ) a nk

(k )

a nn

b 1(1) ⎤(2) ⎥b 2

⎥ ⎥ (k ) ⎥b k ⎥ ⎥(k ) ⎥b n ⎦

在方框内的诸元素中选取绝对值最大的元素为主元素,保证:

(k )

|a i (k , ) k |=max |a ij |≠0; 从而确定 i k .

k

k ≤i ≤n

If i k ≠ k then 交换第 k 行与第i k 行; 然后进行消元,如此进行,直至k =n -1为止. 3.1 列主元素Gauss 消去法[算法] Step1[消元]: 对k =1,2,…,n -1 ①

(k ) (k )

|a |=max |a 选主元 确定 i k ∈{k , , n }, 满足i , k ij |≠0;

k

k ≤i ≤n

(k )

② 若i k , k

a =0; 则停止计算,det A =0.

③ If i k ≠ k then 交换第 k 行与第i k 行; ④ 消元 对i =k +1,k+2,…,n

计算因子m ik

=a

(k )

ik

/a

(k ) kk ;

对j =k +1,k+2,…,n

(k +1) (k ) (k ) ⎧a ij =a ij -m ik a kj 计算⎨(k +1)

(k ) (k )

b =b -m b ⎩i i ik k

;

Step2[回代]: ①

(n )

a 若nn =0则停止计算,det A =0.

② 对i =n , n -1,…,1

b x i =

(n )

i (n )

-∑a ij x j

n

a

j =i +1

(n ) ii

(4) 例题分析:

⎧0. 001x 1+2. 000x 2+3. 000x 3=1. 000⎪

求解方程组:⎨-1. 000x 1+3. 712x 2+4. 623x 3=2. 000

⎪-2. 000x -1. 070x +5. 643x =3. 000⎩123

解之得:X *=(-0.479107 -0.033089 0.355552) T HW: 3.5

作业讲评3

⎧5x 1+2x 2+x 3=-12⎪

[3.1] 设有方程组⎨-x 1+4x 2+2x 3=20

⎪2x -3x +10x =3⎩123

(1) 考察用Jacobi'Method 、Gauss-Seidel'Method 解方程组的收敛性.

用Jacobi'Method 、Gauss-Seidel'Method 解方程,要求当|| x (k +1)-x (k ) ||∞

21⎤⎡5

⎢解:(1) 由于方程组系数矩阵A=-142⎥是一个严格对角占优矩⎢⎥⎢⎣2-310⎥⎦

阵,故用Jacobi'Method 、Gauss-Seidel'Method 进行迭代求解时算法均收敛.

(2) 用Jacobi'Method. 据此建立迭代公式:

(k ) (k )

⎧x 1(k +1) =-0. 4x 2-0. 2x 3-2. 4⎪(k +1) (k ) (k ) x =0. 25x -0. 5x +5 ⎨213

⎪x (k +1) =-0. 2x (k ) -0. 3x (k ) +0. 3

12⎩3

取迭代初值x 1

(0)

(0) (0)

=x 2=x 3=0, 其计算结果如表一.

Jacobi'Method 计算结果(表一)

迭代次数

0 1 2 3 4 5 6 7 8

x1 0 -2.4 -4.46 -4.556 -3.9914 -3.85794 -3.97123 -4.03031 -4.01386

x2 0 5 4.25 2.745 2.6275 2.9848 3.09225 3.02368 2.981464

x3 0 0.3 2.28 2.467 2.0347 1.88653 1.967028 2.02192 2.013165

10 11 12 13 14 15 16 17 18

-3.99542 -4.00024 -4.00122 -4.0002 -3.99973 -3.99989 -4.00005 -4.00004

-4

3.00259 3.003129 3.000009 2.9992 2.999826 3.000167 3.000081 2.999974 2.999974

1.99603 1.999862 2.000987 2.000247 1.9998 1.999894 2.000028 2.000033

2

利用Gauss-Seidel'Method ,据此建立迭代:

(k ) (k )

⎧x 1(k +1) =-0. 4x

2-0. 2x 3-2. 4⎪(k +1) (k +1) (k ) x =0. 25x -0. 5x +5 ⎨213

⎪x (k +1) =-0. 2x (k +1) -0. 3x (k +1) +0. 3

12⎩3

取迭代初值x 1

(0)

(0) (0) =x 2=x 3=0, 其计算结果如表二.

Gauss-Seidel'Method 计算结果(表二)

5 6 7 8

-4.00451 -3.99931 -3.99997 -4.00003

2.998116 3.000003 3.000075 2.999983

2.000338 1.999864 2.000017 2.000002

⎧a 11x 1+a 12x 2=b 1

[3.2] 设有方程组⎨

⎩a 21x 1+a 22x 2=b 2

(a 11, a 12≠0)

⎧(k ) 1(k -1) x =(b -a x ) 11122⎪a 11⎪

迭代公式为:⎨(k =1, 2, 3...)

1(k ) ⎪x 2=(b 2-a 21x 1(k -1) )

⎪a 22⎩

求证由上述迭代公式产生的向量序列{X (k ) }收敛的充要条件是

a 12a 21

γ=

a 11a 22

证明: 显然,上述迭代格式属于Jacobi 迭代格式,其迭代矩阵为X (k ) =BX (k -1) +f,

⎡⎢0

其中,B =⎢

⎢-a 21⎢⎣a 22

a 12

-⎥a 11

⎥,由迭代法基本定理得: 0⎥⎥⎦

2

X

(k )

a 12a 21

→X ⇔ρ(B )

a 11a 22

*

a 12a 21

γ

=

a 11a 22

[3.3] 用SOR 方法解下列方程组(取松弛因子ω

(k )

=1. 2), 要求 || x (k +1)-x -

||∞

-4

⎧2x 1+x 2=1, ⎨

⎩x 1-4x 2=5

.

解: SOR方法是Gauss-Seidel 法的一种改进(修正).

(k +1) (k ) ⎧=-0. 5x 2-0. 5⎪x 1

Gauss-Seidel'Method 迭代格式为:⎨,

(k +1) (k +1) ⎪x =0. 25x +1. 25⎩21

因此,SOR 法的迭代式为:

x

(k +1)

i

=(1-ω) x

(k ) i

+ω(k +1) i

=x i (k ) +ω(取迭代初值x 1

(0)

(k +1) i

-x i (k ) ) i =1, 2.

(0) =x 2=0, 其计算结果如表三.

SOR'Method 计算结果(表三)

次数 x1 x2 ||0 1 2 3 4 5 6

0 0.6 1.272 0.85824 1.0713408 0.9642927 1.0178566

X (k ) -X (k -1) ||∞

0 -1.32 -0.8544 -1.071648 -0.964268 -1.017859 -0.991071

1.32 0.672 0.41376 0.2131008 0.1070481 0.0535638

7 8 9 10 11 12 13 14 15 16

0.9910715 1.0044643 0.9977679 1.0011161 0.999442 1.000279 0.9998605 1.0000698 0.9999651 1.0000174

-1.004464 -0.997768 -1.001116 -0.999442 -1.000279 -0.99986 -1.00007 -0.999965 -1.000017 -0.999991

0.0267851 0.0133928 0.0066964 0.0033482 0.0016741 0.0008371 0.0004185 0.0002093 0.0001046 5.232E-05

3.7

1 矩阵A 的LU 分解

:

已给n 阶方阵A , 若能求得一个下三角方阵L 和一个上三角方阵

U , 使得A=LU, 则我们称方阵A 有LU 三角分解.

由高斯消去法,我们知道它是通过逐步消元过程,将方程组的系数矩阵A 转变为一个上三角矩阵,这实际上相当于用一系列初等矩阵左乘A . 2 高斯消去法的矩阵形式: Step1:第一次消元(a 11

(1)

≠0) :

(1) 1) a 12... a 1(n (1) (1) a 22... a 2n

(1)

⎡a 11⎢(1) a 21(1) (1) ⎢[A b ]=⎢ ⎢(1) ⎢⎣a n 1(1) ⎡a 11⎢⎢0⎢ ⎢⎢⎣0

(1) (1)

a n ... a 2nn (1) 1) a 12... a 1(n (2) (2) a 22... a 2n

(2) (2) a n ... a 2nn

b 1(1) ⎤

(1) ⎥b 2⎥⇒ ⎥⎥(1) b n ⎥⎦

b 1(1) ⎤(2) ⎥b 2⎥

=[A (2) b (2) ] ⎥⎥(2) b n ⎥⎦

即相当于:

0 0⎤⎡1

(1) (1) ⎢-m ⎥⎧m =a /a i 1i 111, 1 0

⎥ 其中, ⎨L 1=⎢21

⎩i =2, 3, , n . ⎥⎢

⎢⎥-m 0 1⎣n 1⎦

(1) (1) 1)

⎡a 11a 12... a 1(n ⎢(2) (2) 0a ... a 222n

L 1[A (1) b (1) ]=⎢

⎢ ⎢(2) (2) 0a ... a ⎢⎣n 2nn

.

b 1(1) ⎤

(2) ⎥b 2⎥

=[A (2) b (2) ] ⎥⎥(2) b n ⎥⎦

Step k :第k 次消元(a kk

L k [A

(k )

(k )

≠0) :

b (k ) ]=[A (k +1) b (k +1) ], 其中,

⎡1

⎢0⎢⎢ ⎢0L k =⎢

⎢0⎢⎢ ⎢ ⎢⎢⎣0

01

1-m k +1, k

m n , k

0 0

0000

(n -1)

0⎤ ⎥

0 0⎥

10 0⎥

010 ⎥00 ⎥

0 01⎥⎦

Step n -1:第n -1次消元(a n -1, n -1

≠0) :

(1) (1) 1)

⎡a 11a 12 a 1(n b 1(1) ⎤⎢(2) (2) (2) ⎥0a 22 a 2n b 2⎥(1) (1) ⎢L n -1L n -2 L 2L 1[A b ]==[A (n ) b (n ) ]⎢ ⎥⎢(n ) (n ) ⎥ 0a nn b n ⎦⎣011(n ) (n )

⇒[A (1) b (1) ]=L -11L -2 L n --21L n -[A b ]-1

⎧L =L L L ⎨(n ) U =A ⎩

其中

-1

1-12-1n -2

L

-1n -1

于是可以推出A =LU .

00 0⎤⎡1

⎢m ⎥1 0⎢21⎥

11

L =L -11L -2 L n --21L n -m 321 0⎥. -1=⎢m 31

⎢⎥ ⎢⎥⎢⎣m n 1m n 2 1⎥⎦

由上述讨论可知,高斯消去法实质上产生了一个将系数矩阵A 分解为上三角阵与下三角阵相乘的因式分解. 若A 的所有顺序主子式均不为0,则 A 的 LU 分解唯一(其中 L 为单位下三角阵). 设有方程组AX=b,并设A=LU, 于是 AX=LUX=b 其中,

UX=Y,

则 LY=b.

于是求解AX=b的问题等价于求解两个方程组UX=Y和LY=b. 具体的解法如下:

(1) 利用顺推过程解LY=b,其计算公式为:

y i =b i -∑l ij y j

j =1

i -1

(i =1, 2, , n ) .

(2) 利用回代过程解UX=Y,其计算公式为:

x i =(y i -∑u ij x j ) /u ii

j =i +1

n

(i =n , n -1, , 1) .

上述方法称为求解线性方程组的三角直接分解法. 这种分 解又称为Doolittle 分解法. 3 Doolittle 分解法[算法] Step1[分解]: ① 对i =1,2,…,n

u 1i =a 1i l i 1=a i 1/u 11

对r =2,3…,n

;

② 计算U 的第r 行, L 的第r 列元素

u ri =a ri -∑l rk u ki

k =1r -1k =1

r -1

(i =r , r +1, , n )

l ir =(a ir -∑l ik u kr ) /u rr

Step2[顺推过程]: 求解LY=b

(i =r , r +1, , n , 且r ≠n )

y i =b i -

i -1

j =1

l ij y j

(i =1, 2, , n )

Step3[回代过程]: 回代过程解UX=Y

x i =(y i -∑u ij x j ) /u ii

j =i +1

n

(i =n , n -1, , 1) .

4 算例

用Doolittle 分解法解方程组

⎡123⎤⎡x 1⎤⎡14⎤⎢252⎥⎢x ⎥=⎢18⎥ ⎢⎥⎢2⎥⎢⎥⎢⎣315⎥⎦⎢⎣x 3⎥⎦⎢⎣20⎥⎦

解: 用Doolittle 算法计算得:

3⎤⎡100⎤⎡12

A =⎢210⎥⎢01-4⎥=LU

⎢⎥⎢⎥⎢⎣3-51⎥⎦⎢⎣00-24⎥⎦

解得LY=(14,18,20)T , 得Y =(14,-10,-72)T

UX=(14,-10,-72)T , 得X =(1,2,3)T

3.8 追赶法

1 三对角方程组

具有如下形式的方程组:

⎡b 1c 1⎤⎡x 1⎤⎡d 1⎤⎢a b ⎥⎢x ⎥⎢d ⎥c 2⎢22⎥⎢2⎥⎢2⎥

⎢⎥⎢ ⎥=⎢ ⎥ ⎢⎥⎢⎥⎢⎥

a b c n -1n -1n -1⎥⎢⎢⎥⎢⎥

⎢a n b n ⎥⎣⎦⎢⎣x n ⎥⎦⎢⎣d n ⎥⎦

称为三对角方程组.

特点:其系数矩阵为一种带状的稀疏矩阵,非零元素集中分布在主对角线及相邻两条次对角线上,且系数矩阵为严格对角占优阵,即

⎧|b 1|>|c 1|

⎨|b i |>|a i |+|c i |⎪|b |>|a |⎩n n

a i ≠0, c i ≠0, i =2, 3, , n -1.

利用高斯消元法,经过n-1次消元后,可得等价的方程组:

⎡1u 1⎤⎡x 1⎤⎡q 1⎤⎢01u ⎥⎢x ⎥⎢q ⎥

2

⎢⎥⎢2⎥⎢2⎥⎢ ⎥⎢ ⎥=⎢ ⎥ ⎢⎥⎢x ⎥⎢q ⎥

1u n -1⎥⎢n -1⎥⎢n -1⎥⎢

⎢1⎥⎣⎦⎢⎣x n ⎥⎦⎢⎣q n ⎥⎦

其中,

u 1=c 1/b 1, q 1=d 1/b 1,

u i =c i /(b i -u i -1a i ) (i =2, 3, , n -1) q i =(d i -x i -1a i ) /(b i -u i -1a i ) (i =2, 3, , n )

利用回代依次求出u i

⇐追的过程

, q i , i =1, 2, , n ,于是,

⎧x n =q n ⎨

⎩x i =q i -u i x i +1(i =n -1, n -2,... , 1)

HW: 3.6 3.7 3.8(希望上机实习)

⇐赶的过程

3.9 其它应用

1 计算|A |

设A =(aij ) n : a) det(A )=det(A T ); b) 数a 乘A 的一行得:det c) A 的两行互换得:det

=a det(A ); =-det(A );

=det(A );

d) A 的一行乘以a 加到另一行得:det e) A 的两行成比例:det(A )=0;

f) det(AB )=det(A ) ·det(B ); 其中B=(b ij ) n

由以上定理可知,通过高斯消元法的计算可得到行列式 的值.

例1 用列主元素法求det(A ) 的值,其中

⎡11-3-2⎤A =⎢-23111⎥

⎢⎥

-22⎥⎢⎣1⎦

1) (2) (n -1) (n )

A |=a (11a 22 a n -1, n -1a n , n ,因此,

解:由矩阵A 的LU 分解过程,可知|

若用列主元素法求行列式的值,只须将每一步的主元素相乘即可,当然要注意行列式的值的符号改变. 其计算过程如下所示.

⎡ 11.0000 -3.0000 -2.0000⎤

⎢-23.0000 11.0000 1.0000⎥⇒(行交换)⎢⎥⎢⎣ 1.0000 -2.0000 2.0000⎥⎦⎡-23.0000 11.0000 1.0000⎤

⎢ 11.0000 -3.0000 -2.0000⎥⇒(消元) ⎢⎥⎢⎣ 1.0000 -2.0000 2.0000⎥⎦⎡-23.0000 11.0000 1.0000⎤

⎢ 0.0000 2.2609 -1.5217⎥⇒(消元) ⎢⎥⎢⎣ 0.0000 -1.5217 2.0435⎥⎦⎡-23.0000 11.0000 1.0000⎤

⎢ 0.0000 2.2609 -1.5217⎥⇒|A |=

⎢⎥⎢⎣ 0.0000 0.0000 1.0192⎥⎦

1 计算A -1

在某些应用中,如在统计学中,可能还需要计算矩阵A 的逆,并且将它明显地表示为A -1.

1.1 利用A 的LU 分解计算A -1 设A =(a ij ) n 为满秩矩阵,则

AX=I, (1) 这里I 为单位矩阵,显然X 为A 的可逆矩阵A -1. 将方程(1)改写为

A [X (1), X (2),…,X (n ) ]=[I (1), I (2),…,I (n ) ] (2) 其中,X (j ) , I (j ) 分别表示X 和I 的第j 列.

于是,方程(2)又可改写为n 个线性方程组的形式: AX (j ) =I(j ) , 1≤

j ≤n (3)

n 个方程组的系数矩阵相同,故可应用LU 分解法来进行计算,这样A -1=[X (1), X (2),…,X (n ) ].并且能够极大地节省计算工作量.

1.2 利用高斯消元法计算A -1

⎡-18-2⎢例如:对矩阵A =-649-10⎢⎢⎣-434-5

解: ⎤⎥,求A -1. ⎥⎥⎦

⎡-18-2100⎤

[A I ]=⎢-649-10010⎥⎢⎥⎢⎣-434-5001⎥⎦ ⎡ 1 0 0 95 -28 18⎤

=⎢ 0 1 0 10 -3 2⎥⎢⎥⎢⎣ 0 0 1 -8 2 -1⎥⎦

故 ⎡ 95 -28 18⎤A -1=⎢ 10 -3 2⎥ ⎢⎥⎢⎣-8 2 -1⎥⎦

3.10 误差分析

1 问题的提出

设方程组AX =b , 其中,A =(a ij ) n 为非奇异阵,X =(x 1, x 2,…,x n ) T , b =(b 1, b 2,…,b n ) T . 由于原始数据a ij , b i 往往是观测数据,难免带有误差,因此,我们下面讨论原始数据的微小变化对方程组的影响.

2 例题

11⎤⎡1⎛1⎫ ⎪⎢ 23 4⎥2 ⎪x 1⎫⎢⎥⎛⎪ 1⎪⎢ 11 1⎥ x 2⎪= ⎪ ⎢345⎥ ⎪3x 3⎭ ⎪⎢1⎥⎝111⎪ ⎢⎥ ⎪456⎝4⎭⎣⎦

的准确解为x *=(1, 0, 0) T , 当向量b 以较小的扰动时, 即

111b =(+ε, +ε, +ε) T , 这时方程组的准确解为234

x *=(1+492ε, -1860ε, 1500ε) T , 说明右端项的微小变化引起了解的很大扰动,其原因是由方程组本身的状态所决定的.

下面分别讨论右端项b i 的误差对解的影响以及系数矩阵元素a ij 的误差对解的影响.

2 右端项b i 的误差对解的影响

设 A 精确,b 有误差δb ,得到的解x +δx , 即A (x +δx ) =b +δb ⇒δx =A δb -1

⇒||δx ||≤||A -1||⋅||δb ||

而 ||b ||=||Ax ||≤||A ||⋅||x ||

||δx ||||δb ||-1≤||A ||⋅||A ||⋅, ||x ||||b ||

||δb ||-1上式说明右端项的相对误差在解中放大了||A ||⋅||A ||倍. ||b ||

3 系数矩阵元素a ij 的误差对解的影响

设b 精确,A 有误差δA ,得到的解为x +δx ,即

(A +δA )(x +δx ) =b

A (x +δx ) +δA (x +δx ) =b

⇒δx =-A -1δA (x +δx )

||δx ||⇒≤||A -1||⋅||δA ||||x +δx ||

||δA ||=||A ||⋅||A ||⋅||A ||-1

或者,

(A +δA ) x +(A +δA ) δx =b

⇒A (I +A -1δA ) δx =-δAx

(只要||δA ||充分小,使得||A -1δA ||≤||A -1||⋅||δA ||

⇒δx =-(I +A -1δA ) -1A -1δAx

||δx ||||A -1||⋅||δA ||⇒≤||x ||1-||A -1||⋅||δA ||

||δA || -1||A ||⋅||A ||⋅||A ||=-11-||A ||⋅||A ||⋅||A ||

||δA ||充分小,矩阵A ||δA ||的相对误差在解中可能放大了||A ||

||A ||⋅||A -1||倍

.

称cond(A )=||A ||⋅||A -1||为矩阵A 的条件数.

当cond(A )>>1时,则方程组是“病态”的;当cond(A ) 较小时,则方程组是“良态”的. 通常的条件数有:

(1) cond(A ) ∞=

(2) cond(A ) 2=||A ||∞⋅||A -1||∞ ||A ||2⋅||A -1||2=λ(A T A ) /λ(A T A ) max min

max |λ|特别地,若 A 对称,则cond (A ) 2=. min |λ|

3 例题

0. 99⎫⎛1已知A = ⎪⎝0. 990. 98⎭, 求A 的条件数.

⎧λ1=1. 980050504解: 由det(λI -A ) =0⇒⎨ , 于是

⎩λ2=-0. 000050504

λ1“病态”cond (A ) 2=≈39206>>1.说明由A 构成的系数矩阵方程组是λ2

的.

3.11 总结

[高斯消去法] 是解线性方程组直接方法的基础. 将线性方程组约化为等价的三角形方程组再求解是直接法的基本解法. 在约化过程中,引进选主元素的技巧是为了保证方法的数值稳定性所采取的必要措施. 如全选主元素消去法;列选主元素消去法等.

[直接三角分解法] 是高斯消去法的变形. 从代数上看,直接三角分解法和高斯消去法本质上是一致的. 但从实际应用效果来看是有差异的. 如用Doolittle 分解法解具有相同系数矩阵但右端向量不同的方程组AX=B=(b 1, b 2, „, b m ) 是相当便利的,每解一个方程组AX=bi 仅需增加n 2次乘除法运算.

迭代法是一种逐次逼近方法,注意到在使用迭代法时,X k +1=BX (k ) +f ,其迭代矩阵B 和迭代向量f 在计算过程中始终不变,迭代法具有循环的计算公式、方法简单. 此外,应注意收敛性与收敛速度问题. 收敛性是迭代法的前提,针对不同的问题,分析并采用适当的数值算法,如Guass-Seidel 方法、SOR 方法等.

对以上算法的分析,立足点是在计算机上实现. 因此,我们对于方法的掌握不仅在数学推导和数学公式上,而且应当深入思考方法的计算机实现过程,以加深对数值计算的认识和理解. □

第三章 线性方程组的解法

3.0 引 言

重要性:解线性代数方程组的有效方法在计算数学和科学计算中具有特

殊的地位和作用. 如弹性力学、电路分析、热传导和振动、以及社会科学及定量分析商业经济中的各种问题.

分类:线性方程组的解法可分为直接法和迭代法两种方法.

(a) 直接法:对于给定的方程组,在没有舍入误差的假设下,能在预定的运算次数内求得精确解. 最基本的直接法是Gauss 消去法,重要的直接法全都受到Gauss 消去法的启发. 计算代价高.

(b ) 迭代法:基于一定的递推格式,产生逼近方程组精确解的近似序列. 收敛性是其为迭代法的前提,此外,存在收敛速度与误差估计问题. 简单实用, 诱人.

(AX =b )

1

基本思想:

与解f (x )=0 的不动点迭代相类似, 将AX =b 改写为X =BX +f 的形式, 建立雅可比方法的迭代格式:X k +1=BX (k ) +f , 其中,B 称为迭代矩阵. 其计算精度可控, 特别适用于求解系数为大型稀疏矩阵(sparse matrices)的方程组. 2

问题:

(a) 如何建立迭代格式?

(b) 向量序列{X k }是否收敛以及收敛条件? 3 例题分析:

⎧10x 1-x 2-2x 3=7. 2

考虑解方程组⎨-x 1+10x 2-2x 3=8. 3 (1)

⎪-x -x +5x =4. 2⎩123

其准确解为X *={1, 1.2, 1.3}. 建立与式(1)相等价的形式:

⎧x 1=0. 1x 2+0. 2x 3+0. 72

⎨x 2=0. 1x 1+0. 2x 3+0. 83 (2) ⎪x =0. 1x +0. 2x +0. 84⎩312

据此建立迭代公式:

(k ) k

⎧x 1(k +1) =0. 1x 2+0. 2x 3+0. 72⎪(k +1) (k ) (k ) x =0. 1x +0. 2x +0. 83 (3) ⎨213

⎪x (k +1) =0. 1x (k ) +0. 2x (k ) +0. 84

12⎩3

取迭代初值x

(0)

1

=x

(0) 2

=x

(0) 3

=0, 迭代结果如下表.

x1 x2 x3

0 0 0 0 1 0.72 0.83 0.84 2 0.971 1.07 1.15 3 1.057 1.1571 1.2482 4 1.08535 1.18534 1.28282 5 1.095098 1.195099 1.294138 6 1.098338 1.198337 1.298039 7 1.099442 1.199442 1.299335 8 1.099811 1.199811 1.299777 9 1.099936 1.199936 1.299924 10 1.099979 1.199979 1.299975 11 1.099993 1.199993 1.299991 12 1.099998 1.199998 1.299997 13 1.099999 1.199999 1.299999 14 1.1 1.2 1.3 15 1.1 1.2 1.3

4 Jocobi 迭代公式:

设方程组AX=b, 通过分离变量的过程建立 Jocobi 迭代公式,即

i =1

∑a ij x j =b i , a ii ≠0(i =1, 2, , n )

n 1

x i =(b i -∑a ij x j ) (i =1, 2, , n )

j =1a ii

j ≠i

n

由此我们可以得到Jacobi 迭代公式:

x

(k +1)

i

n 1

=(b i -∑a ij x i k ) (i =1, 2, , n )

j =1a ii

j ≠i

[Jacobi迭代公式的算法] 1: 初始化. n , (a ij ), (b j ), (x 1) , M . 2: 执行k =1直到M 为止. ①

执行i =1直到n 为止.

u i ←(b i -∑a ij x j ) /a ii ;

j =1j ≠i

n

② ③

执行i =1直到n 为止.

x i ←u i ;

输出k , (x i ).

另外,我们也可以建立Jacobi 迭代公式的矩阵形式. 设方程组AX =b ,

其中,A =(a ij ) n 为非奇异阵,

X =(x 1, x 2,…,x n ) T , b =(b 1, b 2,…,b n ) T

将系数阵A 分解为: A =U +D +L ,

U 为上三角矩阵,D 为对角矩阵,L 为下三角矩阵.

AX =b 可改写为 (U +D +L ) X =b

⇔ X =D -1b -D -1(U +L ) X

由此可得矩阵形式的Jocobi 迭代公式: X k +1=BX (k ) +f

§3.2

注意到利用Jocobi 迭代公式计算x i

(k +1)

时,已经计算好

(k ) k ) x 1(k ) , x 2, , x i (-1的值,而Jocobi 迭代公式并不利用这些最新的近似值计

算,仍用x 1

(k )

(k ) k )

, x 2, , x i (-1. 这启发我们可以对其加以改进, 即在每个分量的

计算中尽量利用最新的迭代值,得到

x

(k +1)

i

i -1n 1(k +1)

=(b i -∑a ij x j -∑a ij x k ) (i =1, 2, , n ) j

j =1j =i +1a ii

上式称为Gauss-Seidel 迭代法. 其矩阵形式是

X =-(D +L ) -1UX +(D +L ) -1b , X k +1=BX (k ) +f .

迭代次数 x1 x2 x3 0 0 0 0 1 0.72 0.902 1.1644 2 1.04308 1.167188 1.282054 3 1.09313 1.195724 1.297771

4 1.099126 1.199467 1.299719 5 1.09989 1.199933 1.299965 6 1.099986 1.199992 1.299996 7 1.099998 1.199999 1.299999 8 1.1 1.2 1.3

.3

1

基本思想:

逐次超松弛迭代法(Successive Over Relaxation Method,简写为SOR) 可以看作带参数ω的高斯-塞德尔迭代法,是G-S 方法的一种修正或加速. 是求解大型稀疏矩阵方程组的有效方法之一. 2 SOR 算法的构造:

设方程组AX =b , 其中,A =(a ij ) n 为非奇异阵,X =(x 1, x 2,…,x n ) T , b =(b 1, b 2,…,b n ) T . 假设已算出x (k ) ,

(k +1) i

i -1n 1(k +1)

=(b i -∑a ij x j -∑a ij x k j ) (i =1, 2, , n ) (1)

j =1j =i +1a ii

相当于用高斯-塞德尔方法计算一个分量的公式. 若对某个参数ω,作(k +1)

与i

x i (k ) 加权的平均, 即

(k +1)

i

x

k +1i

=(1-ω) x

(k ) i

+ω=x i (k ) +ω((k +1) i

-x i (k ) )

(2)

其中,ω称为松弛因子.

用(1)式代入(2)式,就得到解方程组AX =b 的逐次超松弛迭代公式:

⎧x i (k +1) =x i (k ) +∆x i ⎪

i -1n ⎨ (3) ω(k +1) (k )

a ij x j -∑a ij x j ) (i =1, 2, , n ) ⎪∆x i =a (b i -∑j =1j =i

⎩ii

显然,当取ω=1时,式(3)就是高斯-塞德尔迭代公式. 3 例题分析:

SOR 方法解方程组

⎧4x 1-2x 2-x 3=0⎪

⎨-2x 1+4x 2-2x 3=-2 (1) ⎪-x -2x +3x =3⎩123

其准确解为X *={1, 1, 2}. 建立与式(1)相等价的形式:

⎪x 1=0. 5x 2+0. 25x 3⎪

⎨x 2=0. 5x 1+0. 5x 3-0. 5 (2) ⎪12

⎪x 3=x 1+x 2+1⎩33

据此建立迭代公式:

⎧(k +1)

(k ) k

=0. 5x 2+0. 25x 3⎪x 1

⎪(k +1) (k )

⎨x 2=0. 5x 1(k ) +0. 5x 3-0. 5 (3)

⎪1(k ) 2(k ) (k +1)

⎪x 3=x 1+x 2+1⎩33

(0) (0) (0)

利用SOR 算法,取迭代初值x 1=x 2=x 3=1,

ω=1.5,迭代结果如下表.

逐次超松弛迭代法

次数 x1 x2 x3 1 0.625000 0.062500 1.750000 2 0.390625 0.882813 1.468750 3 1.017578 0.516602 1.808594

4 0.556885 0.880981 1.710449 5 1.023712 0.743423 1.868103 6 0.746250 0.908419 1.838737 7 0.997715 0.860264 1.913894 8 0.864050 0.936742 1.908605 9 0.986259 0.922225 1.945523 10 0.928110 0.958649 1.947493 11 0.985242 0.955944 1.966198 12 0.961661 0.973818 1.969521 13 0.988103 0.974699 1.979289 14 0.979206 0.983746 1.982172 15 0.991521 0.985318 1.987416 16 0.988509 0.990038 1.989513 17 0.994341 0.991414 1.992397 18 0.993538 0.993946 1.993806 19 0.996367 0.994950 1.995424 20 0.996313 0.996342 1.996331 21 0.997724 0.997018 1.997254 22 0.997871 0.997798 1.997822 23 0.998596 0.998234 1.998355

GS 迭代法须迭代85次得到准确值X *={1, 1, 2};而SOR 方 法只须55次即得准确值. 由此可见,适当地选择松弛因子

ω,SOR 法具有明显的加速收敛效果. □

3.4

1. 向量和矩阵范数 (a) 向量范数

R n 空间的向量范数 || · || ,对任意

, ∈R n , 满足下列条件:

(1) ||x ||≥0; ||x ||=0⇔x =0 (正定性)

(2) ||αx ||=|α|⋅||x || (齐次性)

(3) ||x +y ||≤||x ||+||y || (三角不等式)

常见的向量范数有: (1) 列范数:

(2) 谱范数:(欧几里德范数或向量的长度, 模

)

(3) 行范数:

(4) p 范数:

:

||x ||∞=max(|x 2-x 1|,|y 2-y 1|) ; ||x ||1=|x 2-x 1|+|y 2-y 1| ;

||x ||2=(x 2-x 1) 2+(y 2-y 1) 2.

}依坐标收敛于向量x * 的充要条件是向量序列

{x (k ) }依范数收敛于向量x *, 即lim ||x (k ) -x *||=0.

k →∞

向量序列{x

(k )

(b) 矩阵范数

|| ,对任意 R m ⨯n 空间的向量范数 || ·

A , B ∈R m ⨯n , 满足下列条件:

(1) ||A ||≥0; ||A ||=0⇔A =0

(2) ||αA ||=|α|⋅||A ||(3) ||A +B ||≤||A ||+||B ||(4) || AB || ≤||A ||||B ||

||A ||∞=max ∑|a ij | (行和范数)

1≤i ≤n

n

j =1

||A ||1=max ∑|a ij | (列和范数)

1≤j ≤n

n

i =1

||A ||2=λm ax (A T A ) (谱范数)

若A 对称,则有

λm ax (A T A ) =λm ax (A 2) .

A ||2=ρ(A ) ,

矩阵A 的谱半径记为||

1≤i ≤n

ρ(A ) =m ax |λi |, 其中λi 为A 的特征根。

2. 迭代法基本定理

设有方程组X =BX +f , 对于任意初始向量X (0)及任意f , 迭代公式

X (k +1)=BX (k ) +f 收敛的充要条件是ρ(B )

对于任意初始向量X (0)及任意f , 迭代公式X (k +1)=BX (k ) +f ,

X (k ) -X *=BX (k -1) +f -(BX *+f )

=B (X (k -1) -X *)

=B 2(X (k -2) -X *) =B (X

k

(0)

-X )

k

*

由此可得, 迭代法收敛的充要条件是B 即,

→0(k →∞) .

ρ(B )

上述定理是线性方程组迭代解法收敛性分析的基本定理,然而由于

ρ(B ) 的计算往往比较困难,尽管有各种办法估计ρ(B ) 的上界,但往往

偏听偏大而不实用,由此导致定理的理论价值胜于实用价值,为满足实际判敛的需要,有如下定理

.

(迭代收敛的充分条件)

设有迭代公式X (k +1)=BX (k ) +f ,如果||B ||

3. 从方程组的系数矩阵A 判断迭代收敛性

实际中要求解的某些线性方程组,其系数矩阵往往具有一些特点,如系数矩阵为对称正定、对角元素占优等. 由这些方程组系数矩阵的特殊性,使得我们可以直接从方程组的系数矩阵A 出发来讨论迭代法的收敛性

.

设A =(a ij ) n ⨯n ∈R

n ⨯n

,满足

|a ii |≥∑|a ij |, i =1, 2, , n

j =1j ≠i

n

且至少有一个i 值,使得

|a ii |>∑|a ij |

j =1j ≠i

n

成立,则称A 为对角占优矩阵;若

|a ii |>∑|a ij |, i =1, 2, , n ,

j =1j ≠i

n

则称A 为严格对角占优矩阵

.

如果A =(a ij ) n ⨯n ∈R

n ⨯n

为严格对角占优矩阵,则对任意的初值

x (0), 解方程组AX =B 的Jacobi 法、Guess-Seidel 迭代法均收敛. HW: 3.1 3.2 3.3(上机实习)

§3.5

1

基本思想:

用高斯消去法求解线性方程组的基本思想是设法消去方程组的系数矩阵A 的主对角线下的元素,而将Ax =b 化为等价的上三角形方程组,然后再通过回代过程便可以获得方程组的解.

这种解线性方程组的方法,常称为高斯消去法(Gaussian Elimination)

.

2 例题分析:

⎧6x 1-2x 2⎪12x -8x ⎪12

利用高斯消去法求解方程组:⎨

⎪3x 1-13x 2⎪⎩-6x 1+4x 2

+2x 3+4x 4=12+6x 3+10x 4=34+9x 3+3x 4=27+x 3-18x 4=-38

2x 2

8x 2 13x 24x 2

+2x 3+4x 4=12+6x 3+10x 4=34+9x 3+3x 4=27+x 3-18x 4=-38

(1)

a i 1

r 1, i =2,3,4.得 利用r i -a 11

⎧6x 1-2x 2⎪ - 4x ⎪2⎨

⎪ ⎪⎩+2x 3+4x 4= 12+2x 3+ 2x 4= 10+8x 3+ x 4= 213x 3-14x 4=-26

(2)

a i (22)

利用r i -(2) r 2, i =3,4.得

a 22

⎧6x 1-2x 2+2x 3+4x 4= 12⎪ - 4x +2x + 2x = 10⎪234

(3) ⎨

⎪ 2x 3- 5x 4= -9⎪⎩ 13x 4=-21

a i (33)

利用r i -(3) r 3, i =4.得

a 33

⎧6x 1-2x 2+2x 3+4x 4= 12⎪ - 4x +2x + 2x = 10⎪234

(4) ⎨

⎪ 2x 3- 5x 4= -9⎪⎩ -3x 4=-3

显然,方程组(4)与(1)是等价的,其系数矩阵为上三角

状的,易于求解. 称以上过程为高斯消去法的消去过程. 通过方程组(4)的回代求解,可以得到准确解为

X *=[1, -3, -2,1]T .

这一过程为高斯消去法的回代过程. 2 高斯消去法算法的构造:

记方程组AX =b 为A (1)X =b (1), 其中,A (1)和b (1)的元素分别记为

(1)

a ij 、b i (1) , i 、j =1, 2, , n .

Step1:第一次消元 设a

(1) 11

1) (1)

倍,(i =2, „, n ), 目≠0,将增广矩阵的第i 行减去m i 1=a i (1/a 11

的是将增广矩阵的第一列内除每一个元素不变外,其余全部消为零,得到A (2)X =b (2), 即

(1)

⎡a 11⎢(1) a 21(1) (1) ⎢[A b ]=⎢ ⎢(1) ⎢⎣a n 1(1) ⎡a 11⎢⎢0⎢ ⎢⎢⎣0

(1) 1)

a 12... a 1(n (1) (1) a 22... a 2n

(1) (1) a n ... a 2nn (1) 1) a 12... a 1(n (2) (2) a 22... a 2n

(2) (2) a n ... a 2nn

b 1(1) ⎤

(1) ⎥b 2⎥⇒ ⎥⎥(1) b n ⎥⎦

b 1(1) ⎤(2) ⎥b 2⎥

=[A (2) b (2) ] ⎥⎥(2) b n ⎥⎦

其中

(2) (1)

⎧a ij =a ij -m i 1a 1(1j )

(i , j =2, ..., n ) ⎨(2) (1) (1)

⎩b i =b i -m i 1b 1

Step2:第k 次消元(2≤k ≤n -1)

(k )

设第k -1次消元已完成,且a kk

≠0,得到A (k ) X =b (k ) , 即

1)

a 1(n

(1) ⎡a 11⎢⎢⎢(k ) (k )

[A b ]=⎢

⎢⎢⎢⎣

(1) a 121)

a 1(k

(2) (2) (2) a 22 a 2 a k 2n

(k ) (k )

a kk a kn

(k ) (k ) a kn a nn

b 1(1) ⎤

(2) ⎥b 2⎥ ⎥ (k ) ⎥b k ⎥ ⎥(k ) ⎥b n ⎦

(k ) (k )

m ik =a ik /a kk

(i =k +1, ..., n ) ,

(k +1) (k ) (k ) ⎧a ij =a ij -m ik a kj ⎨(k +1) (k ) (k ) b =b -m b ⎩i i ik k

(i , j =k +1, ..., n )

如此反复,经过n -1次消元之后得到一个与原方程组等价的上三角形方程组.

(1)

⎡a 11⎢⎢⎢b (n ) ]=⎢

⎢⎢⎢⎣

(1)

a 12(2) a 22

1)

⎤⎡x 1⎤⎡b 1(1) ⎤a 1(n

(2) ⎥⎢⎥⎢b (2) ⎥... a 2x n ⎥⎢2⎥⎢2⎥. ⎥⎢. ⎥⎢. ⎥

=⎢⎥

⎥. . ⎢. ⎥

⎥⎢⎥⎢⎥. ⎥⎢. ⎥⎢. ⎥

(n ) (n )

⎥⎢x ⎢⎥b a nn ⎦⎦⎣n ⎦⎣n ⎥

...

[A (n )

Ste 只要a

回代

(n )

nn

≠0就可以回代求解

n

(n ) (n )

x n =b n /a nn

(i )

b -∑a ij x j (i ) i

x i =

3 高斯消去法[算法] Step1[消元]: 对k =1,2,…,n -1 ① ②

若a kk

(k )

a

j =i +1

(i ) ii

(i =n -1, ..., 1)

=0则停止计算

(k ) (k ) =a ik /a kk ;

对i =k +1,k+2,…,n

计算因子m ik

对j =k +1,k+2,…,n

(k +1) (k ) (k ) ⎧a ij =a ij -m ik a kj 计算⎨(k +1)

(k ) (k )

b =b -m b ⎩i i ik k

;

Step2[回代]: 对i =n , n -1,…,1

(i )

b -∑a ij x j

j =i +1

(i ) ii

(i ) i

n

x i =

a

(高斯消去法的条件)

(1)

若A 的所有顺序主子式均不为0,则高斯消元无需换行即可进行到底,且得到唯一解.

(2)

若消元过程中允许对增广矩阵进行行交换,则方程组Ax =b 可用消去法求解的充要条件是A 可逆.

3.6 高斯主元素消去法

1 主元素及其选取问题

Gauss 消去法第k 次消去是用第k 个方程

(k ) (k ) a kk x k + +a kn x n =b k (k )

来消去第k +1,…,n 个方程中的x k , 条件是a kk 关键元素,称为第k 次消去的主元(素). Gauss 消去法存在的问题是: (1) 顺序消元时一旦产生a kk (2) 此外,即使a kk

(k )

(k )

(k )

(k )

是实现第k 次消元的≠0. a kk

=0(这是经常可能的) ,消元过程则中断;

≠0但绝对值很小时,由于用它作除数,引起在消去过

程中出现数量级及舍入误差急剧增长的系数,而使最后的计算解严重地不可靠. 例:单精度求解方程组

⎧10-9x 1+x 2

+x 2⎩x 1

=1

=2

8 ⎧1

⎪⎪x 1=1-10-9=1. 00... 0100...

其准确解为 ⎨

8 ⎪

⎪⎩x 2=2-x 1=0. 99... 9899...

当利用Gauss 消元法时,

8

9

a 22=1-m 21⨯1=0. 0... 01⨯109-109= -10

b 2=2-m 21⨯1= -10

9

舍入误差

⎡10-911⎤⎡10-911⎤⎢⎥⇒⎢99⎥恶性传播 1120-10-10⎣⎦⎣⎦

⇒x

2=1, x 1=0×

基本思想:

主元素法是对Gauss 消去法的改进. 它全面或局部地选

取绝对值大的元素为主元素,仅对Gauss 消去法的步骤作某些技术性地修改,使之成为一种有效的方法. 从而保证和改善算法的数值稳定性. 2 完全主元素消去法

设方程组AX =b , 其中,A =(a ij ) n 为非奇异阵,X =(x 1, x 2,…,x n ) T , b =(b 1, b 2,…,b n ) T . 经过k -1次选主元消元后,得到下列等价方程组:

(1)

⎡a 11⎢⎢⎢(k ) (k )

[A b ]=⎢

⎢⎢⎢⎣

(1) a 12

1)

a 1(k

1)

a 1(n

(2) (2) (2) a 22 a 2 a k 2n

(k ) (k ) a kk a kn

(k ) a nk

(k )

a nn

b 1(1) ⎤(2) ⎥b 2⎥ ⎥ (k ) ⎥b k ⎥ ⎥(k ) ⎥b n ⎦

① 选主元过程

(k ) (k )

⎡a kk ⎤ a kn ⎢⎥

在矩阵 ⎢⎥中选取绝对值最大的元素为主元素,保证 (k ) (k ) ⎢⎥a a nn ⎦⎣nk

(k ) |a i (k , ) j |=max |a ij |≠0; 从而确定 i k , j k .

k

k

k ≤i , j ≤n

② 行变换和列交换

If i k ≠ k then 交换第 k 行与第i k 行;

j k ≠ k then 交换第 k 列与第j k 列;

值得注意的是,在全主元消去过程中,列交换已改

变了x 各分量的顺序,因此,必须在每次列交换的同时,记录调换后未知数的排列次序. ③ ④ ⑤

消元 回代求解

还原未知数的排列次序

2.1 全主元素Gauss 消去法[算法] Step1[消元]: 对k =1,2,…,n -1 ① ② ③

(k ) (k )

|a |=max |a 选主元 确定 i k , j k ∈{k , , n }, 满足i , j ij |≠0;

k

k

k ≤i , j ≤n

(k )

若i k , j k

a =0; 则停止计算,det A =0.

If i k ≠ k then 交换第 k 行与第i k 行; If j k ≠ k then 交换第 k 列与第j k 列;

④ 消元

对i =k +1,k+2,…,n

计算因子m ik

(k ) (k )

=a ik /a kk ;

对j =k +1,k+2,…,n

(k +1) (k ) (k ) ⎧a ij =a ij -m ik a kj 计算⎨(k +1)

(k ) (k )

b =b -m b ⎩i i ik k

;

Step2[回代]: ①

(n )

若nn

a =0则停止计算,det A =0.

② 对i =n , n -1,…,1

b y i =

(n )

i (n )

-∑a ij x j

n

a

j =i +1

(n ) ii

Step3[还原排列次序]: 对i =1,2,…,n -1

x * := y i

(3) 列主元素消去法

在计算机上实现主元素消去法意味着进行数的比较操作,全选主元素法需要相当多的计算时间,因此常采用局部选主元素的方法.

列主元素消去法依次按列选主元素,只须进行方程行 交换,不产生未知数次序的调换. ① 按列选主元过程 设方程组AX =b 的增广矩阵为

(1) ⎡a 11⎢(1) a 21(1) (1) ⎢ [A b ]=⎢ ⎢(1) ⎣a n 1

(1) a 12(1) a 22 (1) a n 2

1)

a 1(n

(1)

a 2n

(1)

a nn

b 1(1) ⎤⎥b 2(1) ⎥

⎥(1) ⎥b n ⎦

首先在A (1) 中第一列选取绝对值最大的元素为主元素,保证

|a i (1, ) 1|=max |a i (11) |≠0; 从而确定 i k .

k

1≤i ≤n

② 行变换

If i k ≠ 1 then 交换第 1 行与第i k 行;

重复上述过程,设已完成第k -1次按列选主元消元后,得到下列等价方程组:

(1) ⎡a 11⎢⎢⎢(k ) (k )

[A b ]=⎢

⎢⎢⎢⎣

(1) a 12

1)

a 1(k

1)

a 1(n

(2) (2) (2) a 22 a 2 a k 2n

(k ) (k ) a kk a kn

(k ) a nk

(k )

a nn

b 1(1) ⎤(2) ⎥b 2

⎥ ⎥ (k ) ⎥b k ⎥ ⎥(k ) ⎥b n ⎦

在方框内的诸元素中选取绝对值最大的元素为主元素,保证:

(k )

|a i (k , ) k |=max |a ij |≠0; 从而确定 i k .

k

k ≤i ≤n

If i k ≠ k then 交换第 k 行与第i k 行; 然后进行消元,如此进行,直至k =n -1为止. 3.1 列主元素Gauss 消去法[算法] Step1[消元]: 对k =1,2,…,n -1 ①

(k ) (k )

|a |=max |a 选主元 确定 i k ∈{k , , n }, 满足i , k ij |≠0;

k

k ≤i ≤n

(k )

② 若i k , k

a =0; 则停止计算,det A =0.

③ If i k ≠ k then 交换第 k 行与第i k 行; ④ 消元 对i =k +1,k+2,…,n

计算因子m ik

=a

(k )

ik

/a

(k ) kk ;

对j =k +1,k+2,…,n

(k +1) (k ) (k ) ⎧a ij =a ij -m ik a kj 计算⎨(k +1)

(k ) (k )

b =b -m b ⎩i i ik k

;

Step2[回代]: ①

(n )

a 若nn =0则停止计算,det A =0.

② 对i =n , n -1,…,1

b x i =

(n )

i (n )

-∑a ij x j

n

a

j =i +1

(n ) ii

(4) 例题分析:

⎧0. 001x 1+2. 000x 2+3. 000x 3=1. 000⎪

求解方程组:⎨-1. 000x 1+3. 712x 2+4. 623x 3=2. 000

⎪-2. 000x -1. 070x +5. 643x =3. 000⎩123

解之得:X *=(-0.479107 -0.033089 0.355552) T HW: 3.5

作业讲评3

⎧5x 1+2x 2+x 3=-12⎪

[3.1] 设有方程组⎨-x 1+4x 2+2x 3=20

⎪2x -3x +10x =3⎩123

(1) 考察用Jacobi'Method 、Gauss-Seidel'Method 解方程组的收敛性.

用Jacobi'Method 、Gauss-Seidel'Method 解方程,要求当|| x (k +1)-x (k ) ||∞

21⎤⎡5

⎢解:(1) 由于方程组系数矩阵A=-142⎥是一个严格对角占优矩⎢⎥⎢⎣2-310⎥⎦

阵,故用Jacobi'Method 、Gauss-Seidel'Method 进行迭代求解时算法均收敛.

(2) 用Jacobi'Method. 据此建立迭代公式:

(k ) (k )

⎧x 1(k +1) =-0. 4x 2-0. 2x 3-2. 4⎪(k +1) (k ) (k ) x =0. 25x -0. 5x +5 ⎨213

⎪x (k +1) =-0. 2x (k ) -0. 3x (k ) +0. 3

12⎩3

取迭代初值x 1

(0)

(0) (0)

=x 2=x 3=0, 其计算结果如表一.

Jacobi'Method 计算结果(表一)

迭代次数

0 1 2 3 4 5 6 7 8

x1 0 -2.4 -4.46 -4.556 -3.9914 -3.85794 -3.97123 -4.03031 -4.01386

x2 0 5 4.25 2.745 2.6275 2.9848 3.09225 3.02368 2.981464

x3 0 0.3 2.28 2.467 2.0347 1.88653 1.967028 2.02192 2.013165

10 11 12 13 14 15 16 17 18

-3.99542 -4.00024 -4.00122 -4.0002 -3.99973 -3.99989 -4.00005 -4.00004

-4

3.00259 3.003129 3.000009 2.9992 2.999826 3.000167 3.000081 2.999974 2.999974

1.99603 1.999862 2.000987 2.000247 1.9998 1.999894 2.000028 2.000033

2

利用Gauss-Seidel'Method ,据此建立迭代:

(k ) (k )

⎧x 1(k +1) =-0. 4x

2-0. 2x 3-2. 4⎪(k +1) (k +1) (k ) x =0. 25x -0. 5x +5 ⎨213

⎪x (k +1) =-0. 2x (k +1) -0. 3x (k +1) +0. 3

12⎩3

取迭代初值x 1

(0)

(0) (0) =x 2=x 3=0, 其计算结果如表二.

Gauss-Seidel'Method 计算结果(表二)

5 6 7 8

-4.00451 -3.99931 -3.99997 -4.00003

2.998116 3.000003 3.000075 2.999983

2.000338 1.999864 2.000017 2.000002

⎧a 11x 1+a 12x 2=b 1

[3.2] 设有方程组⎨

⎩a 21x 1+a 22x 2=b 2

(a 11, a 12≠0)

⎧(k ) 1(k -1) x =(b -a x ) 11122⎪a 11⎪

迭代公式为:⎨(k =1, 2, 3...)

1(k ) ⎪x 2=(b 2-a 21x 1(k -1) )

⎪a 22⎩

求证由上述迭代公式产生的向量序列{X (k ) }收敛的充要条件是

a 12a 21

γ=

a 11a 22

证明: 显然,上述迭代格式属于Jacobi 迭代格式,其迭代矩阵为X (k ) =BX (k -1) +f,

⎡⎢0

其中,B =⎢

⎢-a 21⎢⎣a 22

a 12

-⎥a 11

⎥,由迭代法基本定理得: 0⎥⎥⎦

2

X

(k )

a 12a 21

→X ⇔ρ(B )

a 11a 22

*

a 12a 21

γ

=

a 11a 22

[3.3] 用SOR 方法解下列方程组(取松弛因子ω

(k )

=1. 2), 要求 || x (k +1)-x -

||∞

-4

⎧2x 1+x 2=1, ⎨

⎩x 1-4x 2=5

.

解: SOR方法是Gauss-Seidel 法的一种改进(修正).

(k +1) (k ) ⎧=-0. 5x 2-0. 5⎪x 1

Gauss-Seidel'Method 迭代格式为:⎨,

(k +1) (k +1) ⎪x =0. 25x +1. 25⎩21

因此,SOR 法的迭代式为:

x

(k +1)

i

=(1-ω) x

(k ) i

+ω(k +1) i

=x i (k ) +ω(取迭代初值x 1

(0)

(k +1) i

-x i (k ) ) i =1, 2.

(0) =x 2=0, 其计算结果如表三.

SOR'Method 计算结果(表三)

次数 x1 x2 ||0 1 2 3 4 5 6

0 0.6 1.272 0.85824 1.0713408 0.9642927 1.0178566

X (k ) -X (k -1) ||∞

0 -1.32 -0.8544 -1.071648 -0.964268 -1.017859 -0.991071

1.32 0.672 0.41376 0.2131008 0.1070481 0.0535638

7 8 9 10 11 12 13 14 15 16

0.9910715 1.0044643 0.9977679 1.0011161 0.999442 1.000279 0.9998605 1.0000698 0.9999651 1.0000174

-1.004464 -0.997768 -1.001116 -0.999442 -1.000279 -0.99986 -1.00007 -0.999965 -1.000017 -0.999991

0.0267851 0.0133928 0.0066964 0.0033482 0.0016741 0.0008371 0.0004185 0.0002093 0.0001046 5.232E-05

3.7

1 矩阵A 的LU 分解

:

已给n 阶方阵A , 若能求得一个下三角方阵L 和一个上三角方阵

U , 使得A=LU, 则我们称方阵A 有LU 三角分解.

由高斯消去法,我们知道它是通过逐步消元过程,将方程组的系数矩阵A 转变为一个上三角矩阵,这实际上相当于用一系列初等矩阵左乘A . 2 高斯消去法的矩阵形式: Step1:第一次消元(a 11

(1)

≠0) :

(1) 1) a 12... a 1(n (1) (1) a 22... a 2n

(1)

⎡a 11⎢(1) a 21(1) (1) ⎢[A b ]=⎢ ⎢(1) ⎢⎣a n 1(1) ⎡a 11⎢⎢0⎢ ⎢⎢⎣0

(1) (1)

a n ... a 2nn (1) 1) a 12... a 1(n (2) (2) a 22... a 2n

(2) (2) a n ... a 2nn

b 1(1) ⎤

(1) ⎥b 2⎥⇒ ⎥⎥(1) b n ⎥⎦

b 1(1) ⎤(2) ⎥b 2⎥

=[A (2) b (2) ] ⎥⎥(2) b n ⎥⎦

即相当于:

0 0⎤⎡1

(1) (1) ⎢-m ⎥⎧m =a /a i 1i 111, 1 0

⎥ 其中, ⎨L 1=⎢21

⎩i =2, 3, , n . ⎥⎢

⎢⎥-m 0 1⎣n 1⎦

(1) (1) 1)

⎡a 11a 12... a 1(n ⎢(2) (2) 0a ... a 222n

L 1[A (1) b (1) ]=⎢

⎢ ⎢(2) (2) 0a ... a ⎢⎣n 2nn

.

b 1(1) ⎤

(2) ⎥b 2⎥

=[A (2) b (2) ] ⎥⎥(2) b n ⎥⎦

Step k :第k 次消元(a kk

L k [A

(k )

(k )

≠0) :

b (k ) ]=[A (k +1) b (k +1) ], 其中,

⎡1

⎢0⎢⎢ ⎢0L k =⎢

⎢0⎢⎢ ⎢ ⎢⎢⎣0

01

1-m k +1, k

m n , k

0 0

0000

(n -1)

0⎤ ⎥

0 0⎥

10 0⎥

010 ⎥00 ⎥

0 01⎥⎦

Step n -1:第n -1次消元(a n -1, n -1

≠0) :

(1) (1) 1)

⎡a 11a 12 a 1(n b 1(1) ⎤⎢(2) (2) (2) ⎥0a 22 a 2n b 2⎥(1) (1) ⎢L n -1L n -2 L 2L 1[A b ]==[A (n ) b (n ) ]⎢ ⎥⎢(n ) (n ) ⎥ 0a nn b n ⎦⎣011(n ) (n )

⇒[A (1) b (1) ]=L -11L -2 L n --21L n -[A b ]-1

⎧L =L L L ⎨(n ) U =A ⎩

其中

-1

1-12-1n -2

L

-1n -1

于是可以推出A =LU .

00 0⎤⎡1

⎢m ⎥1 0⎢21⎥

11

L =L -11L -2 L n --21L n -m 321 0⎥. -1=⎢m 31

⎢⎥ ⎢⎥⎢⎣m n 1m n 2 1⎥⎦

由上述讨论可知,高斯消去法实质上产生了一个将系数矩阵A 分解为上三角阵与下三角阵相乘的因式分解. 若A 的所有顺序主子式均不为0,则 A 的 LU 分解唯一(其中 L 为单位下三角阵). 设有方程组AX=b,并设A=LU, 于是 AX=LUX=b 其中,

UX=Y,

则 LY=b.

于是求解AX=b的问题等价于求解两个方程组UX=Y和LY=b. 具体的解法如下:

(1) 利用顺推过程解LY=b,其计算公式为:

y i =b i -∑l ij y j

j =1

i -1

(i =1, 2, , n ) .

(2) 利用回代过程解UX=Y,其计算公式为:

x i =(y i -∑u ij x j ) /u ii

j =i +1

n

(i =n , n -1, , 1) .

上述方法称为求解线性方程组的三角直接分解法. 这种分 解又称为Doolittle 分解法. 3 Doolittle 分解法[算法] Step1[分解]: ① 对i =1,2,…,n

u 1i =a 1i l i 1=a i 1/u 11

对r =2,3…,n

;

② 计算U 的第r 行, L 的第r 列元素

u ri =a ri -∑l rk u ki

k =1r -1k =1

r -1

(i =r , r +1, , n )

l ir =(a ir -∑l ik u kr ) /u rr

Step2[顺推过程]: 求解LY=b

(i =r , r +1, , n , 且r ≠n )

y i =b i -

i -1

j =1

l ij y j

(i =1, 2, , n )

Step3[回代过程]: 回代过程解UX=Y

x i =(y i -∑u ij x j ) /u ii

j =i +1

n

(i =n , n -1, , 1) .

4 算例

用Doolittle 分解法解方程组

⎡123⎤⎡x 1⎤⎡14⎤⎢252⎥⎢x ⎥=⎢18⎥ ⎢⎥⎢2⎥⎢⎥⎢⎣315⎥⎦⎢⎣x 3⎥⎦⎢⎣20⎥⎦

解: 用Doolittle 算法计算得:

3⎤⎡100⎤⎡12

A =⎢210⎥⎢01-4⎥=LU

⎢⎥⎢⎥⎢⎣3-51⎥⎦⎢⎣00-24⎥⎦

解得LY=(14,18,20)T , 得Y =(14,-10,-72)T

UX=(14,-10,-72)T , 得X =(1,2,3)T

3.8 追赶法

1 三对角方程组

具有如下形式的方程组:

⎡b 1c 1⎤⎡x 1⎤⎡d 1⎤⎢a b ⎥⎢x ⎥⎢d ⎥c 2⎢22⎥⎢2⎥⎢2⎥

⎢⎥⎢ ⎥=⎢ ⎥ ⎢⎥⎢⎥⎢⎥

a b c n -1n -1n -1⎥⎢⎢⎥⎢⎥

⎢a n b n ⎥⎣⎦⎢⎣x n ⎥⎦⎢⎣d n ⎥⎦

称为三对角方程组.

特点:其系数矩阵为一种带状的稀疏矩阵,非零元素集中分布在主对角线及相邻两条次对角线上,且系数矩阵为严格对角占优阵,即

⎧|b 1|>|c 1|

⎨|b i |>|a i |+|c i |⎪|b |>|a |⎩n n

a i ≠0, c i ≠0, i =2, 3, , n -1.

利用高斯消元法,经过n-1次消元后,可得等价的方程组:

⎡1u 1⎤⎡x 1⎤⎡q 1⎤⎢01u ⎥⎢x ⎥⎢q ⎥

2

⎢⎥⎢2⎥⎢2⎥⎢ ⎥⎢ ⎥=⎢ ⎥ ⎢⎥⎢x ⎥⎢q ⎥

1u n -1⎥⎢n -1⎥⎢n -1⎥⎢

⎢1⎥⎣⎦⎢⎣x n ⎥⎦⎢⎣q n ⎥⎦

其中,

u 1=c 1/b 1, q 1=d 1/b 1,

u i =c i /(b i -u i -1a i ) (i =2, 3, , n -1) q i =(d i -x i -1a i ) /(b i -u i -1a i ) (i =2, 3, , n )

利用回代依次求出u i

⇐追的过程

, q i , i =1, 2, , n ,于是,

⎧x n =q n ⎨

⎩x i =q i -u i x i +1(i =n -1, n -2,... , 1)

HW: 3.6 3.7 3.8(希望上机实习)

⇐赶的过程

3.9 其它应用

1 计算|A |

设A =(aij ) n : a) det(A )=det(A T ); b) 数a 乘A 的一行得:det c) A 的两行互换得:det

=a det(A ); =-det(A );

=det(A );

d) A 的一行乘以a 加到另一行得:det e) A 的两行成比例:det(A )=0;

f) det(AB )=det(A ) ·det(B ); 其中B=(b ij ) n

由以上定理可知,通过高斯消元法的计算可得到行列式 的值.

例1 用列主元素法求det(A ) 的值,其中

⎡11-3-2⎤A =⎢-23111⎥

⎢⎥

-22⎥⎢⎣1⎦

1) (2) (n -1) (n )

A |=a (11a 22 a n -1, n -1a n , n ,因此,

解:由矩阵A 的LU 分解过程,可知|

若用列主元素法求行列式的值,只须将每一步的主元素相乘即可,当然要注意行列式的值的符号改变. 其计算过程如下所示.

⎡ 11.0000 -3.0000 -2.0000⎤

⎢-23.0000 11.0000 1.0000⎥⇒(行交换)⎢⎥⎢⎣ 1.0000 -2.0000 2.0000⎥⎦⎡-23.0000 11.0000 1.0000⎤

⎢ 11.0000 -3.0000 -2.0000⎥⇒(消元) ⎢⎥⎢⎣ 1.0000 -2.0000 2.0000⎥⎦⎡-23.0000 11.0000 1.0000⎤

⎢ 0.0000 2.2609 -1.5217⎥⇒(消元) ⎢⎥⎢⎣ 0.0000 -1.5217 2.0435⎥⎦⎡-23.0000 11.0000 1.0000⎤

⎢ 0.0000 2.2609 -1.5217⎥⇒|A |=

⎢⎥⎢⎣ 0.0000 0.0000 1.0192⎥⎦

1 计算A -1

在某些应用中,如在统计学中,可能还需要计算矩阵A 的逆,并且将它明显地表示为A -1.

1.1 利用A 的LU 分解计算A -1 设A =(a ij ) n 为满秩矩阵,则

AX=I, (1) 这里I 为单位矩阵,显然X 为A 的可逆矩阵A -1. 将方程(1)改写为

A [X (1), X (2),…,X (n ) ]=[I (1), I (2),…,I (n ) ] (2) 其中,X (j ) , I (j ) 分别表示X 和I 的第j 列.

于是,方程(2)又可改写为n 个线性方程组的形式: AX (j ) =I(j ) , 1≤

j ≤n (3)

n 个方程组的系数矩阵相同,故可应用LU 分解法来进行计算,这样A -1=[X (1), X (2),…,X (n ) ].并且能够极大地节省计算工作量.

1.2 利用高斯消元法计算A -1

⎡-18-2⎢例如:对矩阵A =-649-10⎢⎢⎣-434-5

解: ⎤⎥,求A -1. ⎥⎥⎦

⎡-18-2100⎤

[A I ]=⎢-649-10010⎥⎢⎥⎢⎣-434-5001⎥⎦ ⎡ 1 0 0 95 -28 18⎤

=⎢ 0 1 0 10 -3 2⎥⎢⎥⎢⎣ 0 0 1 -8 2 -1⎥⎦

故 ⎡ 95 -28 18⎤A -1=⎢ 10 -3 2⎥ ⎢⎥⎢⎣-8 2 -1⎥⎦

3.10 误差分析

1 问题的提出

设方程组AX =b , 其中,A =(a ij ) n 为非奇异阵,X =(x 1, x 2,…,x n ) T , b =(b 1, b 2,…,b n ) T . 由于原始数据a ij , b i 往往是观测数据,难免带有误差,因此,我们下面讨论原始数据的微小变化对方程组的影响.

2 例题

11⎤⎡1⎛1⎫ ⎪⎢ 23 4⎥2 ⎪x 1⎫⎢⎥⎛⎪ 1⎪⎢ 11 1⎥ x 2⎪= ⎪ ⎢345⎥ ⎪3x 3⎭ ⎪⎢1⎥⎝111⎪ ⎢⎥ ⎪456⎝4⎭⎣⎦

的准确解为x *=(1, 0, 0) T , 当向量b 以较小的扰动时, 即

111b =(+ε, +ε, +ε) T , 这时方程组的准确解为234

x *=(1+492ε, -1860ε, 1500ε) T , 说明右端项的微小变化引起了解的很大扰动,其原因是由方程组本身的状态所决定的.

下面分别讨论右端项b i 的误差对解的影响以及系数矩阵元素a ij 的误差对解的影响.

2 右端项b i 的误差对解的影响

设 A 精确,b 有误差δb ,得到的解x +δx , 即A (x +δx ) =b +δb ⇒δx =A δb -1

⇒||δx ||≤||A -1||⋅||δb ||

而 ||b ||=||Ax ||≤||A ||⋅||x ||

||δx ||||δb ||-1≤||A ||⋅||A ||⋅, ||x ||||b ||

||δb ||-1上式说明右端项的相对误差在解中放大了||A ||⋅||A ||倍. ||b ||

3 系数矩阵元素a ij 的误差对解的影响

设b 精确,A 有误差δA ,得到的解为x +δx ,即

(A +δA )(x +δx ) =b

A (x +δx ) +δA (x +δx ) =b

⇒δx =-A -1δA (x +δx )

||δx ||⇒≤||A -1||⋅||δA ||||x +δx ||

||δA ||=||A ||⋅||A ||⋅||A ||-1

或者,

(A +δA ) x +(A +δA ) δx =b

⇒A (I +A -1δA ) δx =-δAx

(只要||δA ||充分小,使得||A -1δA ||≤||A -1||⋅||δA ||

⇒δx =-(I +A -1δA ) -1A -1δAx

||δx ||||A -1||⋅||δA ||⇒≤||x ||1-||A -1||⋅||δA ||

||δA || -1||A ||⋅||A ||⋅||A ||=-11-||A ||⋅||A ||⋅||A ||

||δA ||充分小,矩阵A ||δA ||的相对误差在解中可能放大了||A ||

||A ||⋅||A -1||倍

.

称cond(A )=||A ||⋅||A -1||为矩阵A 的条件数.

当cond(A )>>1时,则方程组是“病态”的;当cond(A ) 较小时,则方程组是“良态”的. 通常的条件数有:

(1) cond(A ) ∞=

(2) cond(A ) 2=||A ||∞⋅||A -1||∞ ||A ||2⋅||A -1||2=λ(A T A ) /λ(A T A ) max min

max |λ|特别地,若 A 对称,则cond (A ) 2=. min |λ|

3 例题

0. 99⎫⎛1已知A = ⎪⎝0. 990. 98⎭, 求A 的条件数.

⎧λ1=1. 980050504解: 由det(λI -A ) =0⇒⎨ , 于是

⎩λ2=-0. 000050504

λ1“病态”cond (A ) 2=≈39206>>1.说明由A 构成的系数矩阵方程组是λ2

的.

3.11 总结

[高斯消去法] 是解线性方程组直接方法的基础. 将线性方程组约化为等价的三角形方程组再求解是直接法的基本解法. 在约化过程中,引进选主元素的技巧是为了保证方法的数值稳定性所采取的必要措施. 如全选主元素消去法;列选主元素消去法等.

[直接三角分解法] 是高斯消去法的变形. 从代数上看,直接三角分解法和高斯消去法本质上是一致的. 但从实际应用效果来看是有差异的. 如用Doolittle 分解法解具有相同系数矩阵但右端向量不同的方程组AX=B=(b 1, b 2, „, b m ) 是相当便利的,每解一个方程组AX=bi 仅需增加n 2次乘除法运算.

迭代法是一种逐次逼近方法,注意到在使用迭代法时,X k +1=BX (k ) +f ,其迭代矩阵B 和迭代向量f 在计算过程中始终不变,迭代法具有循环的计算公式、方法简单. 此外,应注意收敛性与收敛速度问题. 收敛性是迭代法的前提,针对不同的问题,分析并采用适当的数值算法,如Guass-Seidel 方法、SOR 方法等.

对以上算法的分析,立足点是在计算机上实现. 因此,我们对于方法的掌握不仅在数学推导和数学公式上,而且应当深入思考方法的计算机实现过程,以加深对数值计算的认识和理解. □


相关内容

  • 关于常微分方程求解公式的注记
  • 2008年12月 阴山学刊 Dee.2008第22卷第4期 YINSHANACADEMICJOURNAL V01.22 No.4 关于常微分方程求解公式的注记 穆 勇 (长江大学信息与数学学院.湖北荆州434023) 摘 要:要对常微分方程中的一阶线性齐次方程以及一阶线性非齐次方程的求解公式进行一下 ...

  • 弹性力学及有限元
  • 弹性力学与有限单元法(报告) 姓名: 尚建波 学号: [1**********]4 班级:土木F1307 第一题(20分) 变分法中的δ符号与微积分中的d符号均表示微小变化,请问二者有何关 系?如何理解在理论上有了δ则不需要有d符号. 第二题(20分) 设y=y(x),证明:当y=y(x)满足固定边 ...

  • 线性微分方程组初值问题的求解公式
  • 第! ! 卷第" 期宝鸡文理学院学报#自然科学版$ #$, &-. /0'&120&345&''676&18. 9:0/; (! ! ) &(" %&' (! **! >0. ! **! 年+月 ? ? ? ? ? ? ...

  • 预处理共轭梯度法求解线性方程组
  • [摘 要]针对共轭梯度法求解线性方程组,提出一种预处理思想.基于次思想,首先给出预处理矩阵,然后求解预处理线性方程组,再使用共轭梯度法求解.最后通过几个数值试验,与直接使用共轭梯度法求解线性方程组相比较,本文的方法提高了收敛速度. [关键词]线性方程组,预处理,共轭梯度法 中图分类号:E911 文献 ...

  • 有限元法基本原理与应用
  • 有限元法基本原理与应用 班级 机械2081 姓名 方志平 指导老师 钟相强 摘要:有限元法的基础是变分原理和加权余量法,其基本求解思想是把计算域划分为有限 个互不重叠的单元,在每个单元内,选择一些合适的节点作为求解函数的插值点,将微分 方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成 ...

  • 常微分方程数值解法及其应用
  • 东北师范大学 硕士学位论文 常微分方程数值解法及其应用 姓名:李晓红 申请学位级别:硕士 专业:应用数学 指导教师:李佐锋 20080501 摘要 微分方程初值问题模型是数学建模竞赛中常见的一类数学模型.对于一些简单而典型的微分方程模型,譬如线性方程.某些特殊的一阶非线性方程等是可以设法求出其解析解 ...

  • 常微分方程数值解法的研究
  • 2011年第5期总第107期 高职教育 No.5. 2011Sum 107 常微分方程数值解法的研究 霍晓程 (株洲职业技术学院 湖南株洲 412001) 摘 要:本文简要介绍了常微分方程的初值问题和边值问题,同时对初值问题和边值问题常用的数值解法做了归纳介绍,并对未来的发展做了展望. 关键词:常微 ...

  • 非线性电路
  • 非线性电路学习报告 电路是由电气.电子器件按某种特定的目的而相互连接所形成的系统的总称.当电路中至少存在一个非线性电路元件时(例如非线性电阻.非线性电感元件等) ,其运动规律要由非线性微分方程或非线性算子来描述,我们称之为非线性电路或非线性系统. 一.非线性电路的特点: 1.非线性电路不满足叠加定理 ...

  • 极小曲面的有限元解法
  • 第29卷 第2期 2010年3月许昌学院学报 JOURNAL OF XUCHANG UN I V ERSI TY Vol . 29. No . 2 Mar . 2010 文章编号:1671-9824(2010) 02-0014-03 极小曲面的有限元解法 郑恩伟, 郭延涛 (许昌学院数学科学学院, ...